Notebook

Peter Norvig
5 January 2016
revised 18 May 2018

Four 4s, Five 5s, and Countdown to 2016¶

On January 1, 2016 Alex Bellos posed this New Year's puzzle:

Fill in the blanks so that this equation makes arithmetical sense:

10 ␣ 9 ␣ 8 ␣ 7 ␣ 6 ␣ 5 ␣ 4 ␣ 3 ␣ 2 ␣ 1 = 2016

You are allowed to use only the four basic arithmetical operations: +, -, ×, ÷. But brackets (parentheses) can be used wherever needed. So, for example, the solution could begin

(10 + 9) * (8 ... or 10 + (9 * 8) ...

Let's see if we can solve this puzzle, and some of the related ones from Alex's first and second post.

Four Operations, No Brackets¶

We'll start with a simpler version of the puzzle: a countdown with no brackets.

There are nine blanks, each of which can be filled by one of four operators, so there are 9⁴ = 262,144 possibilities, few enough that we can enumerate them all, using itertools.product to get sequences of operators, and then str.format to plug them into blanks, and then eval to evaluate the string:

In [1]:

from itertools import product
from functools import lru_cache # Used later
from collections import Counter, defaultdict # Used later

In [2]:

ops = next(product(('+', '-', '*', '/'), repeat=9))
ops

Out[2]:

('+', '+', '+', '+', '+', '+', '+', '+', '+')

In [3]:

'10{}9{}8{}7{}6{}5{}4{}3{}2{}1'.format(*ops)

Out[3]:

'10+9+8+7+6+5+4+3+2+1'

In [4]:

eval(_)

Out[4]:

We need to catch errors such as dividing by zero, so I'll define a wrapper function, evaluate, to do that, and I'll define simple_countdown to put the pieces together:

In [5]:

def evaluate(exp):
    "eval(exp), or None if there is an arithmetic error."
    try:
        return eval(exp)
    except ArithmeticError:
        return None

def simple_countdown(target, operators=('+', '-', '*', '/')):
    "All solutions to the countdown puzzle (with no brackets)."
    exps = ('10{}9{}8{}7{}6{}5{}4{}3{}2{}1'.format(*ops)
            for ops in product(operators, repeat=9))
    return [exp for exp in exps if evaluate(exp) == target]

In [6]:

simple_countdown(2016)

Out[6]:

[]

Too bad; we did all that work and didn't find a solution. What years can we find solutions for? I'll modify simple_countdown to take a collection of target years rather than a single one, and return a dict of the form {year: 'expression'} for each expression that evaluates to one of the target years. I'll also generalize it to allow any format string, not just the countdown string.

In [7]:

def simple_countdown(targets, operators=('+', '-', '*', '/'), 
                      fmt='10{}9{}8{}7{}6{}5{}4{}3{}2{}1'):
    "All solutions to the countdown puzzle (with no brackets)."
    exps = (fmt.format(*ops)
            for ops in product(operators, repeat=fmt.count('{}')))
    return {int(evaluate(exp)): exp for exp in exps if evaluate(exp) in targets}

simple_countdown(range(1900, 2100))

Out[7]:

{1981: '10*9*8+7*6*5*4*3/2+1',
 1979: '10*9*8+7*6*5*4*3/2-1',
 1980: '10*9*8+7*6*5*4*3/2/1',
 2019: '10*9*8*7/6/5*4*3+2+1',
 2017: '10*9*8*7/6/5*4*3+2-1',
 2018: '10*9*8*7/6/5*4*3+2/1',
 2015: '10*9*8*7/6/5*4*3-2+1',
 2013: '10*9*8*7/6/5*4*3-2-1',
 2014: '10*9*8*7/6/5*4*3-2/1'}

Interesting: in the 20th and 21st centuries, there are only two "golden eras" where the bracketless countdown equation works: the three year period centered on 1980, and the seven year period that is centered on 2016, but omits 2016.

Four Operations, With Brackets¶

Now we return to the original problem. Can I make use of what I have so far? Well, my second version of simple_countdown accepts different format strings, so I could give it all possible bracketed format strings and let it fill in all possible operator combinations. How long would that take? I happen to remember that the number of bracketings of n operands is the nth Catalan number, which I looked up and found to be 4862. A single call to simple_countdown took about 2 seconds, so we could do 4862 calls in about 160 minutes. That's not terrible, but (a) it would still take work to generate all the bracketings, and (b) I think I can find another way that is much faster.

I'll restate the problem as: given the sequence of numbers (10, 9, ... 1), create an expression whose value is 2016. But to get there I'll solve a more general problem: given any sequence of numbers, like (10, 9, 8), what {value: expression} pairs can I make with them? I'll define expressions(numbers) to return a dict of {value: expression} for all expressions (strings) whose numeric value is value, and can be made from numbers, for example:

expressions((10,)) ⇒ {10: '10'}
expressions((9, 8)) ⇒ {1: '(9-8)', 1.125: '(9/8)', 17: '(9+8)', 72: '(9*8)'}

I'll use the idea of dynamic programming: break the problem down into simpler subparts, compute an answer for each subpart, and remember intermediate results so we don't need to re-compute them later. How do we break the problem into parts? expressions((10, 9, 8)) should consist of all the ways of splitting (10, 9, 8) into two parts, finding all the expressions that can be made with each part, and combining pairs of expressions with any of the four operators:

expressions((10, 9, 8)) ⇒ {11: '(10+(9-8))', 27: '(10+(9+8))', 720: '(10*(9*8))', ...}

First the function splits:

In [8]:

c10 = (10, 9, 8, 7, 6, 5, 4, 3, 2, 1)

def splits(sequence):
    "Split sequence into two non-empty parts, in all ways."
    return [(sequence[:i], sequence[i:]) 
            for i in range(1, len(sequence))]

In [9]:

splits(c10)

Out[9]:

[((10,), (9, 8, 7, 6, 5, 4, 3, 2, 1)),
 ((10, 9), (8, 7, 6, 5, 4, 3, 2, 1)),
 ((10, 9, 8), (7, 6, 5, 4, 3, 2, 1)),
 ((10, 9, 8, 7), (6, 5, 4, 3, 2, 1)),
 ((10, 9, 8, 7, 6), (5, 4, 3, 2, 1)),
 ((10, 9, 8, 7, 6, 5), (4, 3, 2, 1)),
 ((10, 9, 8, 7, 6, 5, 4), (3, 2, 1)),
 ((10, 9, 8, 7, 6, 5, 4, 3), (2, 1)),
 ((10, 9, 8, 7, 6, 5, 4, 3, 2), (1,))]

Now the function expressions:

In [10]:

@lru_cache(None)
def expressions(numbers: tuple) -> dict:
    "Return {value: expr} for all expressions that can be made from numbers."
    if len(numbers) == 1: 
        return {numbers[0]: str(numbers[0])}
    else: 
        result = {}
        for (Lnums, Rnums) in splits(numbers):
            for (L, R) in product(expressions(Lnums), expressions(Rnums)):
                Lexp =  '(' + expressions(Lnums)[L]
                Rexp =        expressions(Rnums)[R] + ')'
                result[L * R] = Lexp + '*' + Rexp
                result[L - R] = Lexp + '-' + Rexp
                result[L + R] = Lexp + '+' + Rexp
                if R != 0: 
                    result[L / R] = Lexp + '/' + Rexp
        return result

For example, at one point in a call to expressions((10, 9, 8)) we would have:

 Lnums, Rnums = (10),   (9, 8)
 L,     R     = 10,     1
 Lexp,  Rexp  = '(10',  '(9-8))'

This would lead to us adding the following four entries to result:

 result[10] = '(10*(9-8))'
 result[9]  = '(10-(9-8))'
 result[11] = '(10+(9-8))'
 result[10] = '(10/(9-8))'

The decorator @lru_cache takes care of storing the intermediate results. Rather than catching errors, we just avoid division by 0.

Let's give it a try:

In [11]:

expressions((10,))

Out[11]:

{10: '10'}

In [12]:

expressions((9, 8))

Out[12]:

{72: '(9*8)', 1: '(9-8)', 17: '(9+8)', 1.125: '(9/8)'}

In [13]:

expressions((10, 9, 8))

Out[13]:

{720: '((10*9)*8)',
 -62: '(10-(9*8))',
 82: '((10*9)-8)',
 0.1388888888888889: '((10/9)/8)',
 10: '(10/(9-8))',
 9: '((10-9)+8)',
 11: '((10+9)-8)',
 170: '(10*(9+8))',
 -7: '((10-9)-8)',
 27: '((10+9)+8)',
 0.5882352941176471: '(10/(9+8))',
 11.25: '((10*9)/8)',
 8.875: '(10-(9/8))',
 11.125: '(10+(9/8))',
 8.88888888888889: '((10/9)*8)',
 98: '((10*9)+8)',
 8: '((10-9)*8)',
 0.125: '((10-9)/8)',
 152: '((10+9)*8)',
 2.375: '((10+9)/8)',
 -6.888888888888889: '((10/9)-8)',
 9.11111111111111: '((10/9)+8)'}

That looks reasonable. Let's solve the whole puzzle.

Countdown to 2016: A Solution¶

In [14]:

%time expressions(c10)[2016]

CPU times: user 29.9 s, sys: 833 ms, total: 30.7 s
Wall time: 30.8 s

Out[14]:

'(((((((((10*9)+8)*7)-6)-5)-4)*3)+2)+1)'

We have an answer! And in seconds, not hours, thanks to dynamic programming! Here are solutions for nearby years:

In [15]:

{y: expressions(c10)[y] for y in range(2010, 2026)}

Out[15]:

{2010: '((((((10+(((9*8)-7)*6))*5)+4)+3)+2)+1)',
 2011: '((((((10+9)*8)+7)*(6+(5*(4/3))))-2)-1)',
 2012: '((((((10+9)*8)+7)*(6+(5*(4/3))))-2)/1)',
 2013: '((((((10+9)*8)+7)*(6+(5*(4/3))))-2)+1)',
 2014: '(((((((10+((9*8)*7))-6)-5)*4)+3)-2)+1)',
 2015: '(((((((((10*9)+8)*7)-6)-5)-4)*3)+2)/1)',
 2016: '(((((((((10*9)+8)*7)-6)-5)-4)*3)+2)+1)',
 2017: '(((10*(((9*8)-((7*6)/(5+4)))*3))-2)-1)',
 2018: '(((10*(((9*8)-((7*6)/(5+4)))*3))-2)/1)',
 2019: '(((10*(((9*8)-((7*6)/(5+4)))*3))-2)+1)',
 2020: '(((((10+((9+((8+7)*6))*5))*4)+3)-2)-1)',
 2021: '((((((((10-9)+(8*7))*6)-5)*4)*3)/2)-1)',
 2022: '((((((((10-9)+(8*7))*6)-5)*4)*3)/2)/1)',
 2023: '(((((((10*9)*(8+7))/6)*(5+4))-3)+2)-1)',
 2024: '(((((((10*9)*(8+7))/6)*(5+4))-3)+2)/1)',
 2025: '(((((((10*9)*(8+7))/6)*(5+4))-3)+2)+1)'}

Counting Solutions¶

Alex Bellos had another challenge:

I was half hoping a computer scientist would let me know exactly how many solutions [to the Countdown problem] there are with only the four basic operations. Maybe someone will.

As it stands, my program can't answer that question, because I only keep one expression for each value.

Also, I'm not sure what it means to be a distinct solution. For example, are ((10+9)+8) and (10+(9+8)) different, or are they same, because they both are equivalent to (10+9+8)? Similarly, are ((3-2)-1) and (3-(2+1) different, or the same because they both are equivalent to (3 + -2 + -1)? I think the notion of "distinct solution" is just inherently ambiguous. My choice is to count each of these as distinct: every expression has exactly ten numbers, nine operators, and nine pairs of brackets, and if an expression differs in any character, it is different. But I won't argue with anyone who prefers a different definition of "distinct solution."

So how can I count expressions? One approach would be to go back to enumerating every equation (all 4862 × 4⁹ = 1.2 bilion of them) and checking which ones equal 2016. That would take about 40 hours with my Python program. A better approach is to mimic expressions, but to make a table of counts, rather than expression strings. I want:

counts[(10, 9, 8)][27] == 2

because there are 2 ways to make 27 with the numbers (10, 9, 8), namely, ((10+9)+8) and (10+(9+8)). And in general, if there are Lcount ways to make L using Lnums and Rcount ways to make R using Rnums then there are Lcount * Rcount ways of making L * R using Lnums followed by Rnums. So we have:

In [16]:

@lru_cache(None)
def counts(numbers: tuple) -> Counter:
    "Return a Counter of {value: count} for every value that can be made from numbers."
    if len(numbers) == 1: # Only one way to make an expression out of a single number
        return Counter(numbers)
    else: 
        result = Counter()
        for (Lnums, Rnums) in splits(numbers):
            for L in counts(Lnums):
                for R in counts(Rnums):
                    count = counts(Lnums)[L] * counts(Rnums)[R]
                    result[L + R] += count
                    result[L - R] += count
                    result[L * R] += count
                    if R != 0:
                        result[L / R] += count
        return result

In [17]:

counts((2, 2)) # corresponds to {0: '2-2', 1.0: '2/2', 4: '2+2' or '2*2'}

Out[17]:

Counter({4: 2, 0: 1, 1.0: 1})

In [18]:

counts((10, 9, 8))[27] # (10+9)+8 or 10+(9+8)

Out[18]:

Looks good to me. Now let's see if we can answer the question.

In [19]:

counts(c10)[2016]

Out[19]:

This says there are 30,066 distinct expressions for 2016.

But we're forgetting about round-off error.

Let's find all the values that are very near to 2016:

In [20]:

{y: counts(c10)[y] 
 for y in counts(c10)
 if abs(y - 2016) < 1e-10}

Out[20]:

{2016.0000000000002: 5792,
 2016.0: 30066,
 2015.9999999999995: 1930,
 2015.999999999997: 15,
 2016.0000000000005: 510,
 2016.0000000000018: 264,
 2016.0000000000023: 12,
 2015.9999999999998: 5868,
 2015.9999999999993: 14,
 2016.000000000002: 18,
 2015.999999999999: 10}

I suspect that all of these actually should be exactly 2016. To be absolutely sure, I could re-do the calculations using exact rational arithmetic, as provided by the fractions.Fraction data type. From experience I know that would be an order of magnitude slower, so instead I'll just add up the counts:

In [21]:

sum(_.values())

Out[21]:

I have more confidence in this answer, 44,499, than in 30,066, but I wouldn't accept it as definitive until it was independently verified and passed an extensive test suite. And of course, if you have a different definition of "distinct solution," you will get a different answer.

Four 4s¶

Alex Bellos continues with a related puzzle:

The most famous “fill in the gaps in the equation” puzzle is known as the four fours, because every equation is of the form

4 ␣ 4 ␣ 4 ␣ 4 = X.

In the classic form of the puzzle you must find a solution for X = 0 to 9 using just addition, subtraction, multiplication and division.

This puzzle goes back to a 1914 publication by the "most famous" mathematician/magician W. W. Rouse Ball. The solution is easy:

In [22]:

{i: expressions((4, 4, 4, 4))[i] for i in range(10)}

Out[22]:

{0: '(((4-4)-4)+4)',
 1: '(((4/4)-4)+4)',
 2: '((4/(4+4))*4)',
 3: '(((4+4)+4)/4)',
 4: '(((4-4)/4)+4)',
 5: '(((4*4)+4)/4)',
 6: '(((4+4)/4)+4)',
 7: '((4-(4/4))+4)',
 8: '(((4+4)/4)*4)',
 9: '(((4/4)+4)+4)'}

Note that I didn't do anything special to take advantage of the fact that in (4, 4, 4, 4) the digits are all the same. Happily, lru_cache does that for me automatically! If I split that into (4, 4) and (4, 4), when it comes time to do the second half of the split, the result will already be in the cache, and so won't be recomputed.

New Mathematical Operations¶

Bellos then writes:

If you want to show off, you can introduce new mathematical operations such as powers, square roots, concatenation and decimals, ... or use the factorial symbol, !.

Bellos is suggesting the following operations:

Powers: (4 ^ 4) is 4 to the fourth power, or 256.
Square roots: √4 is the square root of 4, or 2.
Concatenation: 44 is forty-four.
Decimals: 4.4 is four and four tenths.
Factorials 4! is 4 × 3 × 2 × 1, or 24.

There are some complications to deal with:

Irrationals: √2 is an irrational number; so we can't do exact rational arithmetic.
Imaginaries: √-1 is an imaginary number; do we want to deal with that?
Overflow: (10. ^ (9. ^ 8.)), as a float, gives an OverflowError.
Out of memory: (10 ^ (9 ^ (8 * 7))), as an int, gives an OutOfMemoryError (even if your memory uses every atom on Earth).
Infinite expressions: We could add an infinite number of square root and/or factorial signs to any expression: √√√√√√(4!!!!!!!!)...
Round-off error: (49*(1/49)) evaluates to 0.9999999999999999, not 1.

We could try to manage this with symbolic algebra, perhaps using SymPy, but that seems complicated, so instead I will make some arbitrary decisions:

Use floating point approximations for everything, rational or irrational.
Disallow irrationals.
Disallow overflow errors (a type of arithmetic error).
Put limits on what is allow for exponentiation.
Allow only two nested unary operations to avoid infinite expressions.

I'll define the function do to do an arithmetic computation, catch any errors, and try to correct round-off errors. The idea is that since my expressions start with integers, any result that is very close to an integer is probably actually that integer. So I'll correct (49*(1/49)) to be 1.0. Of course, an expression like (1+(10^-99)) is also very close to 1.0, but it should not be rounded off. Instead, I'll try to avoid such expressions by silently dropping any value that is outside the range of 10^-10 to 10¹⁰ (except of course I will accept an exact 0).

In [23]:

from operator import add, sub, mul, truediv
from math import sqrt, factorial

def do(op, *args): 
    "Return op(*args), trying to correct for roundoff error, or `None` on Error or too big/small."
    try:
        val = op(*args)
        return (val        if val == 0 else
                None       if isinstance(val, complex) else
                None       if not (1e-10 < abs(val) < 1e10) else
                round(val) if abs(val - round(val)) < 1e-12 else 
                val)
    except (ArithmeticError, ValueError):
        return None

In [24]:

assert do(truediv, 12, 4) == 3 and do(truediv, 12, 0) == None
assert do(mul, 49, do(truediv, 1, 49)) == 1
assert do(pow, 10, -99) == None

I'll take this opportunity to refactor expressions to have these subfunctions:

digit_expressions(result, numbers): returns a dict like {0.44: '.44', 4.4: '4.4', 44.0: '44'}
add_unary_expressions(result): add expressions like √44 and 4! to result and return result.
add_binary_expressions(result, numbers): add expressions like 4+√44 and 4^4 to result and return result.

In [25]:

@lru_cache(None)
def expressions(numbers: tuple) -> dict:
    "Return {value: expr} for all expressions that can be made from numbers."
    return add_unary_expressions(add_binary_expressions(digit_expressions(numbers), numbers))

def digit_expressions(digits: tuple) -> dict:
    "All ways of making numbers from these digits (in order), and a decimal point."
    D = ''.join(map(str, digits))
    exps = [(D[:i] + '.' + D[i:]).rstrip('.')
            for i in range(len(D) + 1)
            if not (D.startswith('0') and i > 1)]
    return {float(exp): exp for exp in exps}
      
def add_binary_expressions(result: dict, numbers: tuple) -> dict:
    "Add binary expressions by splitting numbers and combining with an op."
    for (Lnums, Rnums) in splits(numbers):
        for (L, R) in product(expressions(Lnums), expressions(Rnums)):
            Lexp = '(' + expressions(Lnums)[L]
            Rexp =       expressions(Rnums)[R] + ')'
            assign(result, do(truediv, L, R), Lexp + '/' + Rexp)
            assign(result, do(mul,     L, R), Lexp + '*'  + Rexp)
            assign(result, do(add,     L, R), Lexp + '+'  + Rexp)
            assign(result, do(sub,     L, R), Lexp + '-'  + Rexp)
            if -10 <= R <= 10 and (L > 0 or int(R) == R):
                assign(result, do(pow, L, R), Lexp + '^' + Rexp)
    return result
                
def add_unary_expressions(result: dict, nesting_level=2) -> dict:
    "Add unary expressions: -v, √v and v! to result"
    for _ in range(nesting_level):
        for v in tuple(result):
            exp = result[v]
            if -v not in result:
                assign(result, -v, '-' + exp)
            if 0 < v <= 100 and 120 * v == round(120 * v): 
                assign(result, sqrt(v), '√' + exp)
            if 3 <= v <= 6 and v == int(v):
                assign(result, factorial(v), exp + '!')
    return result

The function assign will silently drop expressions whose value is None, and if two expressions have the same value, assign keeps the one with the lower "weight"—a measure of the characters in the expression. This shows simpler expressions in favor of complex ones: for four 4s, the entry for 0 will be 44-44, not something like -√((-√4*--4)*(-√4--√4)).

In [26]:

def assign(result: dict, value: float, exp: str): 
    "Assign result[value] = exp, unless we already have a lighter exp or value is None."
    if (value is not None 
        and (value not in result or weight(exp) < weight(result[value]))): 
        result[value] = exp
        
PENALTIES = defaultdict(int, {'√':7, '!':5, '.':1, '^':3, '/':1, '-':1})
        
def weight(exp: str) -> int: return len(exp) + sum(PENALTIES[c] for c in exp)

I'll define a function to print a table of consecutive integers (starting at 0) that can be made by a sequence of numbers:

In [27]:

def show(numbers: tuple, upto=10000, clear=True):
    """Print expressions for integers from 0 up to limit or the first unmakeable integer."""
    if clear: expressions.cache_clear() # To free up memory
    print('{:,} entries for {}\n'.format(len(expressions(numbers)), numbers))
    for i in range(upto + 1): 
        if i not in expressions(numbers):
            return i
        print('{:4} = {}'.format(i, unbracket(expressions(numbers)[i])))
              
def unbracket(exp: str) -> str:
    "Strip outer parens from exp, if they are there"
    return (exp[1:-1] if exp.startswith('(') and exp.endswith(')') else exp)

In [28]:

%time show((4, 4, 4, 4))

705,959 entries for (4, 4, 4, 4)

   0 = 44-44
   1 = 44/44
   2 = 4*(4/(4+4))
   3 = (4+(4+4))/4
   4 = 4+(4*(4-4))
   5 = (4+(4*4))/4
   6 = 4+((4+4)/4)
   7 = (44/4)-4
   8 = 4+(4*(4/4))
   9 = 4+(4+(4/4))
  10 = 44/4.4
  11 = (4/4)+(4/.4)
  12 = (4+44)/4
  13 = 4!-(44/4)
  14 = (4+(.4*4))/.4
  15 = 4+(44/4)
  16 = .4*(44-4)
  17 = (4/4)+(4*4)
  18 = .4+(.4*44)
  19 = (4+(4-.4))/.4
  20 = 4*(4+(4/4))
  21 = (4+4.4)/.4
  22 = √4*(44/4)
  23 = ((4*4!)-4)/4
  24 = 4+(4+(4*4))
  25 = (4+(4*4!))/4
  26 = (4/.4)+(4*4)
  27 = 4+(4!-(4/4))
  28 = 44-(4*4)
  29 = 4+(4/(.4*.4))
  30 = (4+(4+4))/.4
  31 = 4!+((4+4!)/4)
  32 = (4*4)+(4*4)
  33 = 4!+((4-.4)/.4)
  34 = 44-(4/.4)
  35 = 4!+(44/4)
  36 = 44-(4+4)
  37 = ((.4+4!)/.4)-4!
  38 = 44-(4!/4)
  39 = ((4*4)-.4)/.4
  40 = 44-√(4*4)
  41 = (.4+(4*4))/.4
  42 = √4+(44-4)
  43 = 44-(4/4)
  44 = 4+(44-4)
  45 = 44+(4/4)
  46 = 4+(44-√4)
  47 = 4!+(4!-(4/4))
  48 = 4*(4+(4+4))
  49 = 44+(√4/.4)
  50 = (4+(4*4))/.4
  51 = (.4+(4!-4))/.4
  52 = 4+(4+44)
  53 = 4!+(4!+(√4/.4))
  54 = 44+(4/.4)
  55 = 44/(.4+.4)
  56 = 4*(4+(4/.4))
  57 = ((.4+4!)/.4)-4
  58 = ((4^4)-4!)/4
  59 = (4!/.4)-(4/4)
  60 = 44+(4*4)
  61 = (4/4)+(4!/.4)
  62 = (4*(4*4))-√4
  63 = ((4^4)-4)/4
  64 = (4+4)*(4+4)
  65 = (4+(4^4))/4
  66 = √4+(4*(4*4))
  67 = √4+((√4+4!)/.4)
  68 = 4+(4*(4*4))
  69 = (4+(4!-.4))/.4
  70 = (4!+(4^4))/4
  71 = (4!+4.4)/.4
  72 = 4+(4!+44)
CPU times: user 9.72 s, sys: 88.3 ms, total: 9.81 s
Wall time: 9.84 s

Out[28]:

We can also solve the "2016 with four fours" puzzle:

In [29]:

expressions((4, 4, 4, 4))[2016]

Out[29]:

'(4!*(4!+(4!/.4)))'

In a separate video, Alex Bellos shows how to form every integer from 0 to infinity using four 4s, if the square root and log functions are allowed. The solution comes from Paul Dirac (although Dirac originally developed it for the "four 2s" problem).

Donald Knuth has conjectured that with floor, square root and factorial, you can make any positive integer with just one 4.

Below are some popular variants:

Four 2s¶

In [30]:

show((2, 2, 2, 2))

109,747 entries for (2, 2, 2, 2)

   0 = 22-22
   1 = 22/22
   2 = 2+(2*(2-2))
   3 = (2+(2*2))/2
   4 = .2*(22-2)
   5 = 2+(2+(2/2))
   6 = 2*(2+(2/2))
   7 = 2+(2/(.2*2))
   8 = 2+(2+(2*2))
   9 = (22/2)-2
  10 = 22/2.2
  11 = (2/2)+(2/.2)
  12 = (2+22)/2
  13 = 2+(22/2)
  14 = 2+(2+(2/.2))
  15 = (2+(2/2))/.2
  16 = 2*(2*(2*2))
  17 = 22-(√.2^-2)
  18 = 22-(2*2)
  19 = (2+(2-.2))/.2
  20 = 22-√(2*2)
  21 = 22-(2/2)
  22 = 2+(22-2)
  23 = 22+(2/2)
  24 = 2*(2+(2/.2))
  25 = 2/(.2*(.2*2))
  26 = 2+(2+22)
  27 = 22+(√.2^-2)
  28 = 2+(2+(2*2)!)
  29 = 2+(2+(.2^-2))
  30 = (2+(2*2))/.2

Out[30]:

Four 9s¶

In [31]:

show((9, 9, 9, 9))

539,849 entries for (9, 9, 9, 9)

   0 = 99-99
   1 = 99/99
   2 = (99/9)-9
   3 = (9+(9+9))/9
   4 = 9-(9/(.9+.9))
   5 = √9+((9+9)/9)
   6 = ((9+(9+9))/9)!
   7 = 9-((9+9)/9)
   8 = ((9*9)-9)/9
   9 = 9+(9*(9-9))
  10 = 99/9.9
  11 = 9+((9+9)/9)
  12 = (9+99)/9
  13 = 9+(√9+(9/9))
  14 = 9+(9/(.9+.9))
  15 = 9+((9+9)/√9)
  16 = 9+((9/.9)-√9)
  17 = 9+(9-(9/9))
  18 = 99-(9*9)
  19 = 9+(9+(9/9))
  20 = 9+(99/9)
  21 = (9+9.9)/.9
  22 = 9+(√9+(9/.9))
  23 = √9+((9+9)/.9)
  24 = (99/√9)-9
  25 = 9+(√9!+(9/.9))
  26 = (9*√9)-(9/9)
  27 = 9+(9+√(9*9))
  28 = 9+(9+(9/.9))
  29 = 9+((9+9)/.9)
  30 = (9+(9+9))/.9
  31 = (.9+(9*√9))/.9
  32 = (99-√9)/√9
  33 = √9*(99/9)
  34 = (√9+99)/√9
  35 = (√9!+99)/√9
  36 = 9+(9+(9+9))
  37 = (9/.9)+(9*√9)
  38 = √9!*(√9+(√9/.9))
  39 = 9+(9*(√9/.9))
  40 = (9+(9*√9))/.9
  41 = (.9+(√9!*√9!))/.9
  42 = 9+(99/√9)
  43 = √9+(√9!*(√9!/.9))
  44 = (9*√9!)-(9/.9)
  45 = 9*(9/(.9+.9))
  46 = (9/.9)+(√9!*√9!)
  47 = √9*(9+(√9!/.9))
  48 = √9!*(9-(9/9))
  49 = 9+(√9!*(√9!/.9))
  50 = ((9*√9!)-9)/.9
  51 = 9*(9-(√9/.9))
  52 = √9!*(9-(√9/9))
  53 = (9*√9!)-(9/9)
  54 = 9*((9+9)/√9)
  55 = 99/(.9+.9)
  56 = √9!*(9+(√9/9))
  57 = √9*(9+(9/.9))
  58 = √9!*(9+(√9!/9))
  59 = ((9*√9!)-.9)/.9
  60 = √9*((9+9)/.9)
  61 = (.9+(9*√9!))/.9

Out[31]:

Four 5s¶

In [32]:

show((5, 5, 5, 5))

202,937 entries for (5, 5, 5, 5)

   0 = 55-55
   1 = 55/55
   2 = (5/5)+(5/5)
   3 = (5+(5+5))/5
   4 = ((5*5)-5)/5
   5 = 5+(5*(5-5))
   6 = (55/5)-5
   7 = 5+((5+5)/5)
   8 = 5.5+(.5*5)
   9 = 5+(5-(5/5))
  10 = 55/5.5
  11 = 5.5+5.5
  12 = (5+55)/5
  13 = .5+(.5*(5*5))
  14 = 5+((5-.5)/.5)
  15 = (5*5)-(5+5)
  16 = 5+(55/5)
  17 = 5+(5!/(5+5))
  18 = (5-.5)/(.5*.5)
  19 = (5+(5-.5))/.5
  20 = 5+(5+(5+5))
  21 = (5+5.5)/.5
  22 = 55/(.5*5)
  23 = .5+(5*(5-.5))
  24 = (5*5)-(5/5)
  25 = .5*(55-5)
  26 = (5/5)+(5*5)
  27 = (.5*55)-.5
  28 = .5+(.5*55)
  29 = (5!+(5*5))/5
  30 = 55-(5*5)
  31 = 55-(5!/5)
  32 = ((5+5)/5)^5
  33 = .5*(5!*.55)
  34 = 5+(5+(5!/5))
  35 = 5+(5+(5*5))
  36 = (5!+(.5*5!))/5
  37 = (.5^-5)+√(5*5)
  38 = ((5!/5)-5)/.5

Out[32]:

Five 5s¶

In [33]:

%time show((5, 5, 5, 5, 5), False)

7,624,387 entries for (5, 5, 5, 5, 5)

   0 = 5*(55-55)
CPU times: user 2min 19s, sys: 1.04 s, total: 2min 20s
Wall time: 2min 21s

Countdown to 2018¶

One more thing: On January 1 2018, Michael Littman posted this:

2+0+1×8, 2+0-1+8, (2+0-1)×8, |2-0-1-8|, -2-0+1×8, -(2+0+1-8), sqrt(|2+0-18|), 2+0+1^8, 20-18, 2^(0×18), 2×0×1×8... Happy New Year!

Can we replicate that countdown? For 2018 and for following years?

In [34]:

show((2,0,1,8), 10)

10,303 entries for (2, 0, 1, 8)

   0 = 2*(0*18)
   1 = 2^(0*18)
   2 = 20-18
   3 = 2.0+(1^8)
   4 = 20*(1-.8)
   5 = -2+((0-1)+8)
   6 = (2*(0-1))+8
   7 = ((2*0)-1)+8
   8 = (2*0)+(1*8)
   9 = (2*0)+(1+8)
  10 = 2.0+(1*8)

In [35]:

show((2,0,1,9), 10)

14,957 entries for (2, 0, 1, 9)

   0 = 2*(0*19)
   1 = 20-19
   2 = 2+(0*19)
   3 = 2+(0.1+.9)
   4 = (2*0)+(1+√9)
   5 = 20/(1+√9)
   6 = -2+((0-1)+9)
   7 = (2*(0-1))+9
   8 = ((2*0)-1)+9
   9 = (2*0)+(1*9)
  10 = 20-(1+9)

Now we run into a problem. For the year 2020, the two zeroes don't give us much to work with, and we can't make all the integers in the countdown. But I can turn a 0 into a 1 with cos(0); maybe that will be enough? Let's try it. I'll modify add_unary_expressions to assign values for sin and cos:

In [36]:

from math import sin, cos

def add_unary_expressions(result: dict, nesting_level=2) -> dict:
    "Add unary expressions: -v, √v and v! to result dict."
    for _ in range(nesting_level):
        for v in tuple(result):
            exp = result[v]
            if -v not in result:
                assign(result, -v, '-' + exp)
            if 0 < v <= 100 and 120 * v == round(120 * v): 
                assign(result, sqrt(v), '√' + exp)
            if 3 <= v <= 6 and v == int(v):
                assign(result, factorial(v), exp + '!')
            if abs(v) in (0, 45, 90, 180):
                assign(result, sin(v), 'sin(' + exp + ')')
                assign(result, cos(v), 'cos(' + exp + ')')
    return result

PENALTIES.update(s=1, i=1, n=1, c=1, o=1)

In [37]:

show((2, 0, 2, 0), 10, clear=True) # Clear the cache so we get new results for (2, 0) etc.

8,195 entries for (2, 0, 2, 0)

   0 = 202*0
   1 = 20/20
   2 = 2+(0*20)
   3 = 2+(.02^0)
   4 = .20*20
   5 = (2^0)/.20
   6 = 2*(cos(0)+2.0)
   7 = 2+(cos(0)/.20)
   8 = 2^(cos(0)+2.0)
   9 = (20/2)-cos(0)
  10 = 2/0.20

In [38]:

show((2, 0, 2, 1), 10)

38,875 entries for (2, 0, 2, 1)

   0 = 2*(0*21)
   1 = -20+21
   2 = 2+(0*21)
   3 = 2+((0*2)+1)
   4 = 2.0*(2*1)
   5 = 2.0+(2+1)
   6 = 2.0*(2+1)
   7 = 2+(0.2^-1)
   8 = 2.0^(2+1)
   9 = (20/2)-1
  10 = 20/(2*1)

What's Next?¶

One exercise would be adding even more operators, such as:

Floor and Ceiling: ⌊5.5⌋ = 5, ⌈5.5⌉ = 6
Nth root: 3√8 = 2
Percent: 5% = 5/100
Repeating decimal: .4... = .44444444... = 4/9
Double factorial: 9!! = 9 × 7 × 5 × 3 × 1 = 945
Gamma function: Γ(n) = (n − 1)!
Prime counting function: π(n) = number of primes ≲ n; π(5) = 3
Transcendental functions: besides sin and cos, there's log, tan, arcsin, ...

In the xkcd forum they got up to 298 for five fives, using the double factorial and π functions. What would you like to do?