A digit string is an ordered list of numbers (integers from 0 to 9 in our case), such as 15397. This is not to be interpreted as fifteen thousand three hundred and ninety seven, but as one followed by five, three, nine, and seven.

\(1\times (5-3)+9\div 7\) is an expression of the digit string 15397. Between the digits we place some operation from addition, subtraction, multiplication and division. We are allowed to use an operation as many times as we like (including not at all), to use as many parentheses as we like, but we cannot reorder the digits and every two digits must have an operation between them.

A solution of the digit string is some expression that evaluates to zero. \(9-9\), \(1+2-3\), \(8-2\times 4\) and \(2-6\div (1+2)\) are solutions of 99, 123, 824 and 2612 respectively. Clearly any digit string which contains a zero can be solved using only multiplications. Eg. \(5\times0\times3\) is a solution of 503.\(1\times (5-3)+9\div 7\) is clearly not a solution to 15397, but can you find one? \((1+5+3-9)\times7\) Can you solve 194878? No. You can't. It is impossible.

Coming soon. Actually is not essential though, given the 'When have we found them all?' section, because if our program terminates, we have found them all, and so there are finitely many.

For each $L>0$, $U(L)$ denotes the number of digit strings of length $L$ that are not solvable. The idea is we calculate $U(1)$, then $U(2), U(3)$ and so on.

When $U(L)=0$

This means that any digit string of length $L$ has a solution. Take any digit string of length $>L$. It has a substring which has length $L$, for example its first $L$ digits. But this substring has a solution, and so the string is solvable. Thus no string of length $>L$ is solvable, so we are finished.

Thus if we, for each successive length, count how many strings of that length cannot be solved, we can stop searching once we reach a length for which the counter is zero. We haved proved above this will happen no later than length 20, but hopefully sooner! Also, there is another way we may stop...

When $U(L)=1$

This means that only one digit string of length $L$ cannot be solved. Take any digit string of length $>L$. It has (at least) two substrings of length $L$: the ones starting at the first digit and the second digit. There are two possibilities: both them are the unsolvable of length $L$, or otherwise. If otherwise, then at least one of them is not that unsolvable, so has a solution. Then the string would have a solvable substring, so is solvable itself. But if both of them are the unsolvable of length $L$, then their first digits match. Then the first two digits of the string are equal, so the string is solvable. Thus no string of length $>L$ is solvable and we are done.

We denote a string of length $L$ as being the first $L$ letters of the alphabet. Clearly $E(1)=1$ as $a$ is the only expression of the string $a$. There is no space to place an operation!

Here are the $E(2)=4$ of a string of length 2:

And here are the $E(3)=24$ of a string of length 3. Check that you agree.

$$\begin{align} 
&a+b+c &&a+b-c &&&a+b\times c &&&&a+b\div c \\ 
&a-b+c &&a-b-c &&&a-b\times c &&&&a-b\div c \\ 
&a\times b+c &&a\times b-c &&&a\times b\times c &&&&a\times b\div c \\ 
&a\div b+c &&a\div b-c &&&a\div b\times c &&&&a\div b\div c \\ 
&a\times (b+c) &&a\times (b-c) &&&(a+ b)\times c &&&&(a- b)\times c \\ 
&a\div (b+c) &&a\div (b-c) &&&(a+b)\div c &&&&(a- b)\div c
\end{align} $$

However we see there is a lot of redundancy in this listing. Everytime $+$ appears, there is another string that is identical except that $+$ is replaced by $-$, and vice verasa. The same is true for the pair $\times, \div$. Removing all of the expressions with $-$ or $\div$, we cut down the list by a factor of $2^{L-1}$, a factor of two for each of the $L-1$ operations. So, we can write the above lists more compactly as:

$$\begin{align}
&L=1 &&L=2 &&&L=3 \\
&a &&a+b &&&a+b+c \\
& &&a\times b &&&a+b\times c \\
& && &&&a\times b+ c \\
& && &&&(a+b)\times c \\
& && &&&a\times (b+ c) \\
& && &&&a\times b\times c 
\end{align}$$

You can get back the original lists by including the copies of these with various numbers of $+$ and $\times$ replaced with $-$ and $\div$. Using just $+$ and $\times$ then, how many expressions can we make? I find 22, see that you agree. Do not worry about the column titles, they are explained later when I show you the procedure.

$$\begin{align}
&\text{Split after 1}, + 
&&\text{Split after 1}, \times 
&&&\text{Split after 2} 
&&&& \text{Split after 3} 
\\
&a+b+c+d
&&a\times(b+c+d)
&&&a\times b+c+d
&&&&a\times b\times c+d
\\
&a+b+c\times d
&&a\times(b+c\times d)
&&&a\times b+c\times d
&&&&a\times(b+c)+d
\\
&a+b\times c+d
&&a\times(b\times c+d)
&&&(a+b)\times(c+d)
&&&&(a+b)\times c+d
\\
&a+b\times c\times d
&&a\times b\times c\times d
&&&(a+b)\times c\times d
&&&&(a+b+c)\times d
\\
&a+b\times(c+d)
&&a\times b\times (c+d)
&&&
&&&&(a+b\times c)\times d
\\
&a+(b+c)\times d
&&a\times(b+c)\times d
&&&
&&&&(a\times b+c)\times d
\end{align}$$

This means that $E(4)=2^3 \times 22 = 176$. Good thing I didn't write them all out...

When constructing expressions, we need to know what we need should wrap parentheses around. For example $a+b\times c$ is different to $(a+b)\times c$, but the same as $a+(b\times c)$.

An expression is called free if there is an addition or subtraction not within parentheses. Here are some free expressions:

If a string is not free, we say it is bound. Of the expressions discussed already, the following are bound.

The number of free and bound strings of length $L$ are denoted by $F(L)$ and $B(L)$ respectively. Looking at the previous lists of expressions, check that you agree $F(1)=0, F(2)=2, F(3)=12, F(4)=88$ and $B(1)=0, B(2)=2, B(3)=12, B(4)=88$.

At this point we might conjecture that $F(L)=B(L)$ for every $L>1$. This is proved later :)

coming soon. For now, you can read the python implementation.

We aim to implement the above strategy in Python to find the expressions of length $L$. To do this, we need to know them for length less than $L$. The only expression of length 1 is bound, and we will represent it using a blank symbol, such as a spcae or an underscore, to signify this is where the digit should be placed.

# initialise
blank = ' '
free_exprs = [[], []]
bound_exprs = [[], [blank]]

Then we use the procedure designed above to find the free and bound expressions of length $L$. These will need to be appended to 'free_exprs' and 'bound_exprs' respectively to help find the expressions of length $L+1$ and so on. Because it is called only a handful of times, it does not need to be hyper-optimised, so I write them in these 'extend' blocks for readability.

def find_next_exprs(
  free_exprs, bound_exprs, 
  ops_add, ops_mult):

    L = len(free_exprs)
    new_free_exprs = []
    new_bound_exprs = []

    # split after first digit

    new_free_exprs.extend(
        f"{blank}{op}{right}" 
        for op in ops_add
        for right in free_exprs[-1]
    )
    
    new_free_exprs.extend(
        f"{blank}{op}{right}" 
        for op in ops_add
        for right in bound_exprs[-1]
    )
    
    new_bound_exprs.extend(
        f"{blank}{op}({right})" 
        for op in ops_mult
        for right in free_exprs[-1]
    )
    
    new_bound_exprs.extend(
        f"{blank}{op}{right}" 
        for op in ops_mult
        for right in bound_exprs[-1]
    )

    for i in range(2, L):
      # split after i'th digit

        new_free_exprs.extend(
            f"{left}{op}{right}" 
            for left in bound_exprs[i]
            for op in ops_add
            for right in free_exprs[L-i]
        )
        
        new_free_exprs.extend(
            f"{left}{op}{right}" 
            for left in bound_exprs[i]
            for op in ops_add
            for right in bound_exprs[L-i]
        )

        new_bound_exprs.extend(
            f"({left}){op}({right})" 
            for left in free_exprs[i]
            for op in ops_mult
            for right in free_exprs[L-i]
        )

        new_bound_exprs.extend(
            f"({left}){op}{right}" 
            for left in free_exprs[i]
            for op in ops_mult
            for right in bound_exprs[L-i]
        )

    return new_free_exprs, new_bound_exprs

def print_exprs(ops_add, ops_mult, L_max):
    free_exprs = [[], []]
    bound_exprs = [[], [blank]]
    for L in range(2, L_max):
        new_free_exprs, new_bound_exprs = find_next_exprs(
            free_exprs, bound_exprs, ops_add, ops_mult)
        free_exprs.append(new_free_exprs)
        bound_exprs.append(new_bound_exprs)
    
        a = len(new_free_exprs)
        b = len(new_bound_exprs)
        print()
        print(f"Free: F({L})={a}. Bound: B({L})={b}. Total: E({L})={a+b}")
        print(new_free_exprs)
        print(new_bound_exprs)

>>> print_exprs('+-', '*/', 4)

Free: F(2)=2. Bound: B(2)=2. Total: E(2)=4
['a+a', 'a-a']
['a*a', 'a/a']

Free: F(3)=12. Bound: B(3)=12. Total: E(3)=24
['a+a+a', 'a+a-a', 'a-a+a', 'a-a-a', 
'a+a*a', 'a+a/a', 'a-a*a', 'a-a/a', 
'a*a+a', 'a*a-a', 'a/a+a', 'a/a-a']
['a*(a+a)', 'a*(a-a)', 'a/(a+a)', 'a/(a-a)', 
'a*a*a', 'a*a/a', 'a/a*a', 'a/a/a', 
'(a+a)*a', '(a+a)/a', '(a-a)*a', '(a-a)/a']

Free: F(2)=1. Bound: B(2)=1. Total: E(2)=2
['a+a']                                                                           
['a*a']

Free: F(3)=3. Bound: B(3)=3. Total: E(3)=6
['a+a+a', 'a+a*a', 'a*a+a']
['a*(a+a)', 'a*a*a', '(a+a)*a']

Free: F(4)=11. Bound: B(4)=11. Total: E(4)=22
['a+a+a+a', 'a+a+a*a', 'a+a*a+a', 
'a+a*(a+a)', 'a+a*a*a', 'a+(a+a)*a', 'a*a+a+a', 
'a*a+a*a', 'a*(a+a)+a', 'a*a*a+a', '(a+a)*a+a']
['a*(a+a+a)', 'a*(a+a*a)', 'a*(a*a+a)', 
'a*a*(a+a)', 'a*a*a*a', 'a*(a+a)*a', '(a+a)*(a+a)', 
'(a+a)*a*a', '(a+a+a)*a', '(a+a*a)*a', '(a*a+a)*a']
        

To get an idea, let us run the following code. It is nothing special, just has some pretty print statements thrown in.

def print_num_of_exprs(ops_add, ops_mult, L_max):

    print(f"{'Length':<{6}}"
          f"{'Free':>{10}}"      
          f"{'Bound':>{10}}"  
          f"{'Total':>{10}}"
          )
    print("-" * 36)
    print(f"{1:>6}"
          f"{0:>10}" 
          f"{1:>10}" 
          f"{1:>10}"
          )

    free_exprs = [[], []]
    bound_exprs = [[], [blank]]

    for L in range(2, L_max):
        new_free_exprs, new_bound_exprs = find_next_exprs(
            free_exprs, bound_exprs, ops_add, ops_mult)
        free_exprs.append(new_free_exprs)
        bound_exprs.append(new_bound_exprs)
    
        a = len(new_free_exprs)
        b = len(new_bound_exprs)
        print(f"{L:>6}"
          f"{a:>10}" 
          f"{b:>10}" 
          f"{a+b:>10}"
          )

Running this, we find there are... dangerously many. Quite a few, even. If we need to use the expressions of length 10, it will take forever, there are 100 billion of them! And God help us if we need the longer lengths, I don't even have the memory to store them!

>>> print_num_of_exprs('+-', '*/'', 12)

Length      Free     Bound     Total
------------------------------------
     2         2         2         4
     3        12        12        24
     4        88        88       176
     5       720       720      1440
     6      6304      6304     12608
     7     57792     57792    115584
     8    547712    547712   1095424
     9   5323008   5323008  10646016
    10  52761088  52761088 105522176
MemoryError
            

Searching for this last column on the OEIS suggests that $E(L) =$ A156017$(L-1)$. I later prove this to be the case, from which it follows that

in the sense that $f(n)\sim g(n) \iff \lim\limits_{n\to\infty} \frac{f(n)}{g(n)} = 1$. From this it is immediate that

which is to say that each time we increase the length by one, the number of expressions grows by about a factor of $6+4\sqrt{2} \approx 11.65685$. This grows even faster than the number of digit strings of length $L$! (assuming we allow less than 12 digits.)

Here is the rough blueprint of our method. There are a few things to note. The expressions can be made using the procedure of Chapter 0. We do not even consider digit strings that contain a zero, as they are trivially solvable. For each steo in the while loop, we increment the length by one and are counting the number of arrays of this that do not have a solution. As we proved earlier if this count is either 0 or 1 we can stop.

import itertools
        
base = 10
ops_mult='+-'
ops_mult='*/'
digits = list(range(1, base))

length = 1
count = base-1
total = count 

while count > 1:
    length += 1
    
    # find expressions with new length

    count = 0
    for array in itertools.product(digits, repeat=length):
        if has_no_solution(array):
            count += 1
    total += count

Problem 1 — Division by Zero

Imagine we are trying to check if $7123$ is solvable, and we try the expression $7\div(1+2-3)$. Our poor python script would not be happy, that's for sure.

This problem can easily be avoided, though. If the operation '$a\div b$' raises this error, it is because $b=0$. Clearly then $a\times b=0$. So, if we were to first try $7\times (1+2-3)$, this would give zero. We would stop trying expressions, and so would never try $7\div (1+2-3)$.

So, imagine if we sort the expressions by the number of divisions. Then, we would never try $a\div b$ without first verifying $a\times b$ was non-zero, which means $b$ is non-zero and $a\div b$ is not problematic. The cost to sorting the strings like this is not relatively large, as it needs to be done only once, when the list of expressions is generated.

Problem 2 — Floating Point Errors

The string 13946 is particularly annoying to work with, but serves as a good example. It can be solved by $1-3\div9-4\div6$. (In fact, this is the only solution. It not only needs a division to be solved, it actually needs 2!) You can see this is a solution, right? Well, Python can't...

We must try to avoid these floating point errors somehow. Consider working with pairs of integers $(a,b)$, for $b>0$ and $\text{gcd}(a,b)=1$. In this language, $(a,b)$ is understood as $a/b$. Then $(a,1)$ would be the integer $a$, and $(0,b)$ would be 0. These are the rationals $\mathbb{Q}$, which are closed under the four operations were are using. The operations can then be written purely in terms of integer arithmetic:

(At the end, remove common factors from the numerator and denominator.) Of course, I'm not actually bothered to code all of that up myself. There is a Python module called 'fractions' which does this for us. You can read about it here. '$\text{Fraction}(a,b)$' is an object that behaves just like $(a,b)$ above. We never need to worry about floating point errors again!

Two small comments on this. Firstly, yes, the expressions now take longer to compute. An addition now requires 3 multiplications, an addition, and a gcd calculation. But we need the accuracy this provides. Luckily, we already agreed to sort the expressions by how many divisions they have. So, we can use integer arithmetic on the first group, those with no divisions. Then, use the fractions approach only for the rest. This should salvage back some time.

Secondly, the floating point error problem could be avoided altogether by carefully redefining the game is such a way that at every step in the evaluation of the expression, all values must be integers. For example, you could consider 13946 to be unsolvable, because in the steops of calculating $1-3\div9-4\div6$ you have some non-integers $1/3$ and $2/3$. This is easy to do when actually playing in person. It is not so easy, and certainly would be slow, to ask the program to ensure we have an integer every time we divide. You are welcome to try that yourself :)

Using everything we have talked about, let's put this into action. First, when we generate the new expressions, we need to split them up into those with divisions, and those without. Also, we must order the ones with divisions by how many divisions they have. This can be done with the following function.

def order_exprs_by_divs(exprs, length):
    """
    Given a list of expressions, and their common length, 
    returns two arrays. First, those without division. 
    Second, those with division, sorted by how many divisions are used.
    """
    exprs_wo_div = []
    temp = [[] for _ in range(length)]
    
    for expr in exprs:
        div_count = expr.count('/')
        if div_count == 0:
            exprs_wo_div.append(expr)
        else:
            temp[div_count].append(expr)

    exprs_w_div = sum(temp, [])
    return exprs_wo_div, exprs_w_div

We will need to evaluate these expressions over each array of values. To prepare for this, we replace the blank character in a smart way. What I have chosen to do is replace the i'th appearance of the balnk character with 'x[i]', starting with i=0. Then, we can evaluate the string by just passing in an array of values called x. For example, using a as the blank character, the expression 'a*(a+a)' would then become 'x[0]*(x[1]+x[2])'. Using the array of values x=[1,2,3], the expression would give 1*(2+3)=5. Precompiling the expressions in advance like this also is just way faster than doing each time.

def precompile_exprs(exprs, indices):
    """
    Precompiles the expressions.
    Replaces blank with 'x[{}]',  then '{}' with the integers 0,1,... from indices.
    Must use eval(expr, {}, {x:values}) to evaluate.
    """
    return [compile(expr.replace(blank, 'x[{}]').format(*indices), 'string', 'eval') 
            for expr in exprs]

Next, we need to actually evaluate them when given a digit string, and stop if we hit a zero. We do integer arithmetic on the ones without divisions, and rational arithemtic on the rest. That way, if we find a zero without using a division, we save ourselves from needing to use rational arithmetic.

from fractions import Fraction
def has_no_solution(array, exprs_wo_div, exprs_w_div):
    """
    False if any expression evaluates to 0, True otherwise.
    Checks expressions in order of increasing 
    number of divisions
    """
    namespace = {'x': array}
    for expr in exprs_wo_div:
        if eval(expr, {}, namespace) == 0:
            return False
        
    namespace = {'x': [Fraction(v) for v in array]}
    for expr in exprs_w_div:
        if eval(expr, {}, namespace) == 0:
            return False
        
    return True

With these, out blueprint now looks like this. All that is needed is some print statements to show the count at the end of each loop, but I'll let you design your own.

base = 10
ops_add = '+-'
ops_mult = '*/'


# initial set up
start_message(base)
start_time = time.time()
digits = list(range(1, base))
length = 1
count = base-1
total = count
count_per_length = [0,count]


while count > 1:
    length += 1
    F, B = find_next_exprs(free_exprs, bound_exprs, ops_add, ops_mult)
    free_exprs.append(F)
    bound_exprs.append(B)

    exprs_wo_div, exprs_w_div = order_exprs_by_divs(F+B, length)
    exprs_wo_div = precompile_exprs(exprs_wo_div, list(range(length)))
    exprs_w_div = precompile_exprs(exprs_w_div, list(range(length)))

    count = 0
    for array in itertools.product(digits, repeat=length):
        if has_no_solution(array, exprs_wo_div, exprs_w_div):
            count += 1
    
    count_per_length.append(count)
    total += count
print(count_per_length)