Entropic Thoughts

Failing at Combinatorics with Haskell

Failing at Combinatorics with Haskell

I stumbled across Jezen Thomas article on Solving a Maths Riddle with Bad Haskell recently. The challenge is described as

With the digits 123456789, make them add up to 100. They must stay in the same order but you can use addition, subtraction, multiplication, division, and parentheses. All digits must be used exactly once.

and Thomas, after showing his code and reasoning, concludes that

Running this produces 14 possible answers, and that’s without enumerating all those possible answers that would result from changing the order of operations with brackets.

I was interested in what some solutions might look like when parentheses are used, but I didn’t want to fiddle with parentheses.

Postfix notation to the rescue

There’s a neat trick to avoid fiddling with parentheses when dealing with arithmetic expressions: turn them into postfix notation. One of his examples,

1 * 2 * 3 - 4 * 5 + 6 * 7 + 8 * 9

could in postfix look like1 I say “could” because there are actually several ways to translate expressions like this into postfix, because an expression like 1 2 + 3 + can, thanks to the associative nature of addition, also be written 1 2 3 + +, but we’ll ignore this inconvenient fact to begin with.

1 2 * 3 * 4 5 * - 6 7 * + 8 9 * +

The reason this lets us avoid parentheses is that we specify the order of operations through the order of numbers and operators, rather than with parentheses.

This leads to another bonus of postfix notation: it’s remarkably easy to evaluate. We don’t have to turn to a generic eval function2 Or visit a shunting yard. but can instead write a tiny stack machine:

eval :: [Char] -> Float
eval expr =
  let
    eval' expr stack =
      case (expr, stack) of
        ([], [x]) -> x
        ('+':cs, a:b:ss) -> eval' cs ((b+a):ss)
        ('-':cs, a:b:ss) -> eval' cs ((b-a):ss)
        ('*':cs, a:b:ss) -> eval' cs ((b*a):ss)
        ('/':cs, a:b:ss) -> eval' cs ((b/a):ss)
        (c:cs, ss) -> eval' cs ((read (c:[])):ss)
  in
    eval' expr []

This is not particularly pretty or robust code3 When I originally wrote this I had exhaustive patterns and it returned an Either String Float to make testing and debugging invalid expressions easier, but it’s not needed for the final program which only ever produces valid expressions., but it takes a string representing a valid arithmetic expression in postfix notation, and computes its value.

Generating postfix expressions

We can generate postfix expressions one character at a time. Essentially, at any given point, we have to ask two questions:

  1. Which digits – if any – can we emit?
  2. Which operators – if any – can we emit? and

For both of these questions, we need to know two things: which was the last digit emitted, and how many operators have we emitted so far? The state of our generator will keep track of those two variables.4 Not in a particularly type-safe way, but meh. We’re doing this for fun, after all.

data State = State Int Int

Emitting digits

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.

Emitting operators

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.

Combining emissions

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)

Combinatorial explosion

If we run this for the digits 1–5 only, we get variations of (2+3)×5×4, which is the only way to make those digits into 100:

123+45*** = 100.0
123+4*5** = 100.0
123+4**5* = 100.0
123+*45** = 100.0
123+*4*5* = 100.0
12*3+45** = 100.0
12*3+4*5* = 100.0

It gets to this result in 0.02 seconds on my laptop, and needs to evaluate 3584 expressions, given by the length of the expressions list. We can build a table indicating these metrics as we increase the interval of digits it needs to use.

Digits Expressions Runtime (s)
1–5 3584 0.02
1–6 43008 0.23
1–7 540672 3.13
1–8 7028736 51.22

How long will it take to run for the interval 1–9, which the original question asks about? We could recognise that the number of expressions grows by roughly 13× between each increase, but that’s boring. We can also use combinatorics to figure out an exact value for the number of expressions.

Combinatorics to the rescue

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.

Estimating runtime

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!

Appendix A: Full code

The entire code I wrote didn’t use as pretty variable and function names as what I’ve shown so far, but here it is, in case I made a silly error when transcribing above.

I’ve left the traceShow import there just to be able to mention it. A classic beginner’s problem with Haskell is being unable to write debug prints to trace execution. That’s what traceShow does. It stands for a recognition that some side-effects are harmless during development, and the developer should have access to them. Debug prints are an example.

import Control.Monad (forM_)
import Data.Char (ord, chr)
import Debug.Trace (traceShow)

--------------------------------
-- Constants and tiny helpers
--------------------------------
max_digit = 9
max_ops = max_digit - 1

dig2chr i = chr (ord '0' + i)


--------------------------------
-- State object
--------------------------------
data State = State Int Int

emit_dig (State d o) =
  map (, State (d+1) o) $
    if d == max_digit then mempty else pure (dig2chr (d+1))

emit_op (State d o) =
  map (, State d (o+1)) $
    if o >= d - 1 then mempty else "+-*/"

finished (State d o) =
  d == max_digit && o == max_ops


--------------------------------
-- Combining emissions
--------------------------------
emissions = [emit_dig, emit_op]

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


--------------------------------
-- Evaluation
--------------------------------
eval :: [Char] -> Float
eval expr =
  let
    eval' expr stack =
      case (expr, stack) of
        ([], [x]) -> x
        ('+':cs, a:b:ss) -> eval' cs ((b+a):ss)
        ('-':cs, a:b:ss) -> eval' cs ((b-a):ss)
        ('*':cs, a:b:ss) -> eval' cs ((b*a):ss)
        ('/':cs, a:b:ss) -> eval' cs ((b/a):ss)
        (c:cs, ss) -> eval' cs ((read (c:[])):ss)
  in
    eval' expr []


--------------------------------
-- Filtering and printing
--------------------------------
main =
  let
    expressions = generate (State 0 0)
  in do
    forM_ (filter (\e -> eval e == 100) expressions) $ \e -> do
      putStr e
      putStr " = "
      print (eval e)

It’s a little longer than Jezen Thomas’ code, and not that much better, but it’s what I got. And I realise I really miss writing Haskell. If only it was easier to build and deploy to foreign systems.