Answering the first question is easiest: as long as we haven’t yet emitted the digit 9, we can always emit the next digit – and this is the only valid digit to emit. We encode this logic in Haskell as

emit_digit :: State -> [(Char, State)]
emit_digit (State prev_digit op_count) =
  let
    emission =
      if prev_digit == 9 then
        mempty
      else
        pure (digit_to_char (prev_digit + 1))
  in
    map (, State (prev_digit + 1) op_count) emission

In other words, we figure out if an emission cannot happen (represented by mempty5⁵ In this context, mempty is actually just an empty list [], but using the standard mempty constant makes it more clear that the purpose is to emit nothing.) or emits a pure char based on the prev_digit. Since in Haskell we don’t mutate state directly, this function amends the emission with the new values of the state variables prev_digit and op_count.

The logic for when operators are allowed is only a tiny bit more complicated since it also needs to account for how many digits are emitted – we cannot have more operators than digits, or we generate a malformed expression. Other than that, it follows the same pattern.

emit_operator :: State -> [(Char, State)]
emit_operator (State prev_digit op_count) =
  let
    emission =
      if op_count >= prev_digit - 1 then
        mempty
      else
        "+-*/"
  in
    map (, State prev_digit (op_count + 1)) emission

In the case we have enough operators for the digits emitted, we will not emit any operators (that mempty again), otherwise any operator is fair game.

This function has the same signature and structure because both these functions conform to the same interface: given the current state (prev_digit and ops_count), return a list of possible emissions along with the updated state appropriate for each.

Now that we know what are valid emissions, we can write the function to tie the generation together. This is where we draw the rest of the owl. First we need to know when we are finished – we need a base case for the iteration.

finished (State prev_digit ops_count) =
  prev_digit == 9 && ops_count == 8

Then we can start combining emissions.

emissions = [emit_digit, emit_operator]

generate state =
  if finished state then
    return []
  else do
    (emitted, next_state) <- foldMap ($ state) emissions
    map (emitted :) (generate next_state)

The base case is fairly obvious: if we have nothing more to emit (neither digits nor operators) we return an empty list. (Noting that return in the Haskell case refers to something other than control flow.)

If there are more things to generate, we get all possible emissions by running foldMap over our two emission types (digits and operators) to first give them the current state and then extract a list of all possible emissions. Then the do syntax with the list monad will try all these variations “simultaneously”.6⁶ In terms of execution it happens sequentially, but my intuition for the list monad is that it allows us to run the same computation across multiple values as if they were one value. We continue generation from the next state, but prepend whatever emission we had from among our alternatives.

With this, all expressions is what we get when we start running this method with a state of (0,0), i.e.

expressions = generate (State 0 0)

To print those that evaluate to 100, we filter, loop, and print in our main function.

main =
  forM_ (filter (\e -> eval e == 100) expressions) $ \e -> do
    putStr e
    putStr " = "
    print (eval e)

With \(n\) digits, we know there will be \(n-1\) operators. Each of these operators are selected freely from an alphabet of four, meaning the number of operator strings (ignoring digits) we can see in the expressions is

\[4^{n-1}.\]

These operators must be distributed between and among the digits. The expression will always have \(2n-1\) graphemes, of which \(n-1\) are operators. We actually know the first two graphemes will always be digits because there cannot be operators there, so realistically, there are \(2n-3\) graphemes where we can pretend we can place \(n-1\) operators freely. The number of ways to do this is

\[\binom{2n-3}{n-1}.\]

The real number of ways to place the operators is lower, because we cannot have more operators than digits. If we remember the random walk argument from an earlier article, we can think of this as disallowing all paths that touch below the diagonal. How many paths do that? I don’t think I can explain it, but I have learned that when there will be N coin throws, of which H are heads, there are \(\binom{N}{H+c}\) paths that touch either +c or -c over the diagonal.7⁷ Thanks de Finetti for writing Theory of Probability and thanks spaced repetition for allowing me to learn even some of the parts of it I don’t fully understand! So in this case, we remove the \(\binom{N}{H+1}\) paths that touch just below the diagonal.

In total then, the number of expressions should be given by

\[4^{n-1} \left[\binom{2n-3}{n-1} - \binom{2n-3}{n}\right].\]

If we evalute this, we find out it matches the experimental figures above perfectly8⁸ Remarkably, because I didn’t actually expect it to. My combinatorics is not strong enough that I can have any sort of confidence in my guesses..

Empirically, we would have guessed 1–9 would generate somewhere around 95 million expressions. With combinatorics, we compute the exact number to be 93.7 million and change.

It seems from the above experimental table that the run time per expression increases as we increase the digit interval used, but we can guess it will land somewhere in the interval 8–12 ms per thousand expressions for 9 digits.

Then we multiply! The run time ought to be

\[9.37 \times 10^7 \times (10 \pm 2) \times 10^{-6}\]

seconds. The rpn calculator in Emacs can compute with error forms, so when we type this expression into it as

10p2
e-6*9.37e7*60/

we learn that it will take \(16 \pm 3\) minutes. Extrapolation from the experimental absolute runtimes indicates around 17 minutes. Either way, I won’t hold my breath.

Fortunately, Haskell evaluates everything lazily, so even though we wrote the code as though it would collect into huge lists, it does not – memory requirements stay low, and the code prints out results in a streaming fashion as it finds them. Here are some examples:

12345678-9+*--+** = 100.0
123456+78--//+9*+ = 100.0
123456/78-/9+**+* = 100.0
12345+6*//-78+9*+ = 100.0
12345+*/67**8+9** = 100.0

It did eventually finish with

12/3/4*567+8+9+** = 100.0

after 15 minutes and 49 seconds!