Alex Bellos posed this New Year's puzzle:
Fill in the blanks so that this equation makes arithmetical sense:
10 &square 9 &square 8 &square 7 &square 6 &square 5 &square 4 &square 3 &square 2 &square 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. We'll start with a simpler version of the puzzle.
In [1]:
4 ** 9
Out[1]:
That's small enough that we can just enumerate all the expressions and see if any evaluates to 2016.
We will generate each expression as a string and then evaluate it with eval. But we need to catch errors such as dividing by zero, so we'll define a wrapper function, evaluate:
In [2]:
from __future__ import division
def evaluate(exp):
"eval exp, or return None if there is an arithmetic error."
try:
return eval(exp)
except ArithmeticError:
return None
We'll try each of the four operators in each of the nine blanks, evaluate each resulting expression, and collect the ones that evaluate to 2016:
In [3]:
def solve_no_brackets(ops=('+', '-', '*', '/')):
"All solutions to the countdown puzzle (with no brackets)."
exps = ('10'+a+'9'+b+'8'+c+'7'+d+'6'+e+'5'+f+'4'+g+'3'+h+'2'+i+'1'
for a in ops for b in ops for c in ops
for d in ops for e in ops for f in ops
for g in ops for h in ops for i in ops)
return [exp for exp in exps if evaluate(exp) == 2016]
In [4]:
solve_no_brackets()
Out[4]:
Too bad; we did all that work and didn't find a solution. What years can we find solutions for? Let's modify solve_no_brackets to take a collection of target years, and return a dict of the form {year: "expression"} for each expression that evaluates to one of the target years:
In [5]:
def solve_no_brackets(targets, ops=('+', '-', '*', '/')):
"Solutions to the no-bracket countdown puzzle matching one of targets."
exps = ('10'+a+'9'+b+'8'+c+'7'+d+'6'+e+'5'+f+'4'+g+'3'+h+'2'+i+'1'
for a in ops for b in ops for c in ops
for d in ops for e in ops for f in ops
for g in ops for h in ops for i in ops)
return {int(evaluate(exp)): exp
for exp in exps if evaluate(exp) in targets}
Let's try it:
In [6]:
%time solve_no_brackets(range(1900, 2100))
Out[6]:
Interesting: in the 20th and 21st centuries, there are only two "golden eras" where the countdown equation works: the three year period centered on 1980, and the seven year period that is centered on 2016, but omits 2016. Note it takes about half a minute to do the computation.
Now let's return to the original puzzle, with the brackets. How many ways are there of bracketing an expression with 9 binary operators? I happen to remember that this is given by the Catalan numbers, and we can look it up to find that there are 4862 diffferent bracketing. If we enumerated and evaluated all of them, 4862 times half a minute is 40 hours. I'm impatient, so I'd like a faster approach. One possibility would be to switch to a faster compiled language like C++ or Go. But I'll concentrate on a better algorithm and stick with Python.
I'll use the idea of dynamic programming: break the problem down into simpler subparts, and compute an answer for each subpart, and store intermediate results in a table so we don't need to re-compute them when we need them again.
How do we break the problem into parts? In general, any expression must consist of an operator with two operands (which might in turn be complex subexpressions). For example, a complete countdown expression might be of the form
(10 ... 8) + (7 ... 1)
where (10 ... 8) means some expression that starts with 10 and ends with 8. Of course we need not use '+' as the operator, and we need not split after the 8; we could use any of the four operators and split anywhere. Let's start by defining c10 as the tuple of integers forming the countdown from 10 to 1, and the function splits to split a tuple in all ways:
In [7]:
c10 = (10, 9, 8, 7, 6, 5, 4, 3, 2, 1)
def splits(items):
"Split sequence of items into two non-empty parts, in all ways."
return [(items[:i], items[i:])
for i in range(1, len(items))]
In [8]:
splits(c10)
Out[8]:
Now what I would like to do is build up a table that says, for every subsequence of the numbers, what are the expressions we can make with those numbers, and what do they evaluate to? We'll call the table EXPS. For example, with the subsequence (10, 9, 8), we would have:
EXPS[(10, 9, 8)] = {
27: '((10+9)+8)',
8: '((10-9)*8)',
-7: '(10-(9+8))',
...}
We'll do the same for every other subsequence, for example:
EXPS[(7, 6, 5, 4, 3, 2, 1)] = {
1: '((((7/((6/5)-4))+3)*2)*1)',
2: '((((7/((6/5)-4))+3)*2)+1)',
3.5: '((((7/((6/5)-4))+3)+2)+1)',
4: '((7-((6/5)*((4/3)+2)))+1)',
...}
Once we have the tables for these two subsequences, we can put them together to get the table for the complete countdown(10) by considering all ways of taking a value from the first table, then one of the four operators, then a value from the second table. For example, taking the first entry from each table, and the operator '+', we would have:
EXPS[(10, 9, 8, 7, 6, 5, 4, 3, 2, 1)] = {
28: '(((10+9)+8)+((((7/((6/5)-4))+3)*2)*1))',
...}
I can implement EXPS as a defaultdict of dicts, and define expressions(numbers) to fill in EXPS entries for numbers and all sub-sequences of numbers. Within fill_tables, note that Lnums and Rnums are sequences of numbers, such as (10, 9, 8) and (7, 6). L and R are values, such as 27 and 1. And Lexp and Rexp are strings, such as "((10+9)+8)" and "(7-6)". Rather than catching division-by-zero errors, we just avoid the division when the denominator is 0.
In [9]:
from collections import defaultdict, Counter
import itertools
EXPS = defaultdict(dict) # e.g., EXPS[(10, 9, 8)][27] == '((10+9)+8)'
def expressions(numbers):
"Fill EXPS table for numbers, and all sub-sequences of numbers. Return EXPS[numbers]"
if numbers in EXPS: # Already did the work
pass
elif len(numbers) == 1: # Only one way to make an expression out of a single number
expr(numbers, numbers[0], str(numbers[0]))
else: # Split in all ways; fill tables for left and right; combine tables in all ways
for (Lnums, Rnums) in splits(numbers):
for (L, R) in itertools.product(expressions(Lnums), expressions(Rnums)):
Lexp, Rexp = '(' + EXPS[Lnums][L], EXPS[Rnums][R] + ')'
if R != 0:
expr(numbers, L / R, Lexp + '/' + Rexp)
expr(numbers, L * R, Lexp + '*' + Rexp)
expr(numbers, L - R, Lexp + '-' + Rexp)
expr(numbers, L + R, Lexp + '+' + Rexp)
return EXPS[numbers]
def expr(numbers, value, exp):
"Record exp as an expression with the given value, covering the sequence of numbers."
EXPS[numbers][value] = exp
Let's give it a try:
In [10]:
expressions((10, 9, 8))
Out[10]:
In [11]:
EXPS[(10, 9, 8)][27]
Out[11]:
In [12]:
%time expressions(c10)[2016]
Out[12]:
We have an answer! And in a lot less than 40 hours, thanks to dynamic programming!
Removing unnecessry brackets, this equation is equivalent to:
In [13]:
(10 + 9 * 8 * 7 - 6 - 5) * 4 + 3 + 2 - 1
Out[13]:
Alex Bellos had another challenge:
I was half hoping a computer scientist would let me know exactly how many solutions 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, and each of these questions could reasonably be answered either way. 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 has 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 × 49 = 1.2 bilion of them) and checking which ones equal 2016. That would take about 40 hours with my Python program, but I could get it under 40 minutes in a more efficient language.
Another approach is to count subexpressions as the table is filled in. We won't enumerate all the expressions, just count them. I'll introduce a second table, COUNTS, such that
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)).
How do we compute the counts? By looking at every split and operator choice that can make the value, and summing up (over all of these) the product of the counts for the two sides. For example, there are 2 ways to make 27 with (10 ... 8), and it turns out there are 3526 ways to make 1 with (7 ... 1). So there are 2 × 3526 = 7052 ways to make 28 with this split by adding 27 and 1.
I'll make expr handle COUNTS as well as EXPS. And I'll define clear to clear out the cache of COUNTS and EXPS, so that we can fill the tables with our new, improved entries.
In [14]:
COUNTS = defaultdict(Counter) # e.g., COUNTS[(10, 9, 8)][27] == 2
def expressions(numbers):
"Fill EXPS table for numbers, and all sub-sequences of numbers. Return EXPS[numbers]"
if numbers in EXPS: # Already did the work
pass
elif len(numbers) == 1: # Only one way to make an expression out of a single number
expr(numbers, numbers[0], str(numbers[0]), 1)
else: # Split in all ways; fill tables for left and right; combine tables in all ways
for (Lnums, Rnums) in splits(numbers):
for (L, R) in itertools.product(expressions(Lnums), expressions(Rnums)):
Lexp, Rexp = '(' + EXPS[Lnums][L], EXPS[Rnums][R] + ')'
count = COUNTS[Lnums][L] * COUNTS[Rnums][R]
if R != 0:
expr(numbers, L / R, Lexp + '/' + Rexp, count)
expr(numbers, L * R, Lexp + '*' + Rexp, count)
expr(numbers, L - R, Lexp + '-' + Rexp, count)
expr(numbers, L + R, Lexp + '+' + Rexp, count)
return EXPS[numbers]
def expr(numbers, val, exp, count):
"Fill EXPS[numbers][val] with exp, and increment COUNTS."
EXPS[numbers][val] = exp
COUNTS[numbers][val] += count
def clear(): EXPS.clear(); COUNTS.clear()
In [15]:
clear()
expressions((10, 9, 8))
Out[15]:
In [16]:
COUNTS[(10, 9, 8)]
Out[16]:
In [17]:
COUNTS[(10, 9, 8)][27]
Out[17]:
Looks good to me. Now let's repeat the computation, this time keeping track of COUNTS:
In [18]:
clear()
%time expressions(c10)[2016]
Out[18]:
In [19]:
COUNTS[c10][2016]
Out[19]:
This says there are 30,066 distinct expressions for 2016.
But I don't believe it.
We'll see why in a bit, but first some trivia questions to answer:
In [20]:
len(expressions(c10))
Out[20]:
In [21]:
import itertools
def unmakeable(numbers):
"Smallest positive integer than can't be made by numbers."
return next(i for i in itertools.count(0) if i not in expressions(numbers))
unmakeable(c10)
Out[21]:
In [22]:
4862 * 4 ** 9
Out[22]:
COUNTS table?
In [23]:
sum(COUNTS[c10].values())
Out[23]:
In [24]:
2015 + 1/3 + 1/3 + 1/3
Out[24]:
In [25]:
2015 + 1/3 + 1/3 + 1/3 == 2016
Out[25]:
So there might be perfectly good solutions that are hiding in the EXPS table under 2015.9999999999998 (or some similar number) when they should be exactly 2016. Let's find all the values that are very near to 2016:
In [26]:
{val for val in expressions(c10)
if abs(val - 2016) < 0.0001}
Out[26]:
I suspect that all of these actually should be exactly 2016.
To find out for sure, let's re-do all the calculations using exact rational arithmetic, as provided by the fractions.Fraction data type. From experience I know that this will be an order of magnitude slower, so we're looking at maybe a 20 minute computation.
I'll replace the computation L / R with divide(L, R), which calls Fraction. To mitigate the expense of this computation, I make two optimizations. First, divide replaces a whole fraction, such as Fraction(6, 2), with an int, such as 3. Second, I modify expr to not fill in the EXPS[c10] table, except for EXPS[c10][2016]. The rationale is that we don't need the other entries, and by not storing them, we save a lot of memory, and thus save garbage collection time.
In [27]:
from fractions import Fraction
def expressions(numbers):
"Fill EXPS table for numbers, and all sub-sequences of numbers. Return EXPS[numbers]"
if numbers in EXPS: # Already did the work
pass
elif len(numbers) == 1: # Only one way to make an expression out of a single number
expr(numbers, numbers[0], str(numbers[0]), 1)
else: # Split in all ways; fill tables for left and right; combine tables in all ways
for (Lnums, Rnums) in splits(numbers):
for (L, R) in itertools.product(expressions(Lnums), expressions(Rnums)):
Lexp, Rexp = '(' + EXPS[Lnums][L], EXPS[Rnums][R] + ')'
count = COUNTS[Lnums][L] * COUNTS[Rnums][R]
if R != 0:
expr(numbers, divide(L, R), Lexp + '/' + Rexp, count)
expr(numbers, L * R, Lexp + '*' + Rexp, count)
expr(numbers, L - R, Lexp + '-' + Rexp, count)
expr(numbers, L + R, Lexp + '+' + Rexp, count)
return EXPS[numbers]
def divide(L, R):
"Exact rational division of L/R."
f = Fraction(L, R)
return f.numerator if f.denominator == 1 else f
def expr(numbers, value, exp, count):
"Fill EXPS[numbers][val] with exp, and increment COUNTS."
if numbers == c10 and value != 2016:
return
EXPS[numbers][value] = exp
COUNTS[numbers][value] += count
In [28]:
clear()
%time expressions(c10)[2016]
Out[28]:
In [29]:
COUNTS[c10][2016]
Out[29]:
That did indeed take about ten times longer, but we now have an answer that I have more confidence in (but I wouldn't accept it as definitive until it was independently verified or I had an extensive tst suite). And of course, if you have a different definition of "distinct solution," you will get a different answer.
Now let's turn to another of Alex's puzzles: making 2016 and other target values from a string of four or five 4s, with exponentiation allowed. Exponentiation is tricky for five reasons:
(0 ** -1) is the same as (1 / 0), and gives a ZeroDivisionError.(2 ** (1 / 2)) is an irrational number; so we can't do exact rational arithmetic.(-1 ** (1 / 2)) is an imaginary number, but Python gives a ValueError.(10. ** (9. ** 8.)), as a float, gives a OverflowError.(10 ** (9 ** (8 * 7))), as an int, gives an OutOfMemoryError (even if your memory was expanded to use every atom in the universe).How do we deal with this? We can't do exact rational arithmetic. We could try to do exact algebra, perhaps using SymPy, but that seems difficult
and computationally expensive, so instead I will abandon the goal of exact computation, and do everything in the domain of floats (reluctantly accepting that there will be some round-off errors). We'll coerce numbers to floats when we first put them in the table, and all subsequent operations will be with floats. I define a new function, expr2, to call expr, catching arithmetic errors. Since we are making some rather arbitrary decisions about what expressions are allowed (e.g. imaginary numbers are not), I'll give up on trying to maintain COUNTS.
In [30]:
from operator import add, sub, mul, truediv
def expressions(numbers):
"Fill EXPS table for numbers, and all sub-sequences of numbers. Return EXPS[numbers]"
if numbers in EXPS: # Already did the work
pass
elif len(numbers) == 1: # Only one way to make an expression out of a single number
expr(numbers, float(numbers[0]), str(numbers[0]))
else: # Split in all ways; fill tables for left and right; combine tables in all ways
for (Lnums, Rnums) in splits(numbers):
for (L, R) in itertools.product(expressions(Lnums), expressions(Rnums)):
Lexp, Rexp = '(' + EXPS[Lnums][L], EXPS[Rnums][R] + ')'
expr2(numbers, L, pow, R, Lexp + '**' + Rexp)
expr2(numbers, L, truediv, R, Lexp + '/' + Rexp)
expr2(numbers, L, mul, R, Lexp + '*' + Rexp)
expr2(numbers, L, add, R, Lexp + '+' + Rexp)
expr2(numbers, L, sub, R, Lexp + '-' + Rexp)
return EXPS[numbers]
def expr2(numbers, L, op, R, exp):
"Fill table entries for op(L, R), catching errors."
try:
expr(numbers, op(L, R), exp)
except (ArithmeticError, ValueError):
pass
def expr(numbers, value, exp): EXPS[numbers][value] = exp
Now we can solve the "2016 with five fours" puzzle:
In [31]:
clear()
expressions((4, 4, 4, 4, 4))[2016]
Out[31]:
I'll define a function to create a table of makeable integers:
In [32]:
def makeable(numbers):
"A table of {i: expression} for all integers i from 0 up to first unmakeable."
targets = range(unmakeable(numbers))
return {i: expressions(numbers)[i]
for i in targets}
We'll use this to see if we can solve the "0 to 9 with four fours" puzzle:
In [33]:
makeable((4, 4, 4, 4))
Out[33]:
Yes: we can get 0 to 9 (but not 10).
Now I'll see what integers we can make with five fives. Legend has it that you can get all the way up to 55:
In [34]:
ff = (5, 5, 5, 5, 5)
makeable(ff)
Out[34]:
We didn't get there.
With some research, I see that others who got up to 55 with five 5s used these three concepts:
55.55! = 120 5We'll refactor expressions to call these three new subfunctions:
digit_expressions: For every subsequence of numbers, we'll smush the digits together, and then make a table entry for those resulting digits as an integer, and with a decimal point in each possible position.binary_expressions: The code that previously was the main body of expressions.unary_expressions: Apply the unary operators to every entry in the table. (Because it applies to entries already in the table, make sure to call unary_expressions last.)We'll still do all computation in the domain of floats.
In [35]:
from math import sqrt, factorial
def expressions(numbers):
"Fill EXPS table for numbers, and all sub-sequences of numbers. Return EXPS[numbers]"
if numbers not in EXPS:
digit_expressions(numbers)
binary_expressions(numbers)
unary_expressions(numbers)
return EXPS[numbers]
def digit_expressions(numbers):
"Fill tables with expressions made from the digits of numbers, and a decimal point."
exp = ''.join(str(n) for n in numbers)
expr(numbers, float(exp), exp)
for d in range(len(exp)):
decimal = exp[:d] + '.' + exp[d:]
expr(numbers, float(decimal), decimal)
def binary_expressions(numbers):
"Fill tables with all expressions formed by splitting numbers and combining with an op."
for (Lnums, Rnums) in splits(numbers):
for (L, R) in itertools.product(expressions(Lnums), expressions(Rnums)):
Lexp, Rexp = '(' + EXPS[Lnums][L], EXPS[Rnums][R] + ')'
if 1 <= R <= 10 and (L > 0 or int(R) == R):
expr2(numbers, L, pow, R, Lexp + '**' + Rexp)
expr2(numbers, L, truediv, R, Lexp + '/' + Rexp)
expr2(numbers, L, mul, R, Lexp + '*' + Rexp)
expr2(numbers, L, add, R, Lexp + '+' + Rexp)
expr2(numbers, L, sub, R, Lexp + '-' + Rexp)
def unary_expressions(numbers):
for v in list(EXPS[numbers]):
exp = EXPS[numbers][v]
if v > 0:
expr(numbers, sqrt(v), '√' + exp)
if 2 <= v <= 6 and v == int(v):
expr(numbers, factorial(v), exp + '!')
Now that we have more variety in the types of expressions formed, I want to choose a "good" expression to represent each value. I'll modify expr so that when there are multiple expressions for a value, it chooses the one with the least "weight," where I define weight as the length of the string, plus a penalty of 1 for every square root sign (just because square root feels "heavier" than the other operations):
In [36]:
def expr(numbers, value, exp):
if value not in EXPS[numbers] or weight(exp) < weight(EXPS[numbers][value]):
EXPS[numbers][value] = exp
def weight(exp): return len(exp) + exp.count('√')
We'll try again:
In [37]:
clear()
%time makeable(ff)
Out[37]:
Wow! We almost tripled the 55 goal! I have to say, I would never have come up with the solution for 81 on my own. It works because (√(5 - .5) * √0.5) = √(4.5 * 0.5) = (√2.25) = 1.5.
At the risk of making the computation take even longer, I'm going to add two more operations:
These operations are useful because they produce integers, and our targets are the integers (from 0 up to whatever).
In addition, I'll allow two consecutive applications of unary operators, thus allowing expressions such as
But still not allowing three consecutive applications, such as
To compensate for these new expressions, I'll be pickier about when I allow the square root function to apply.
In [38]:
from math import floor, ceil
def unary_expressions(numbers):
for i in range(2):
for v in list(EXPS[numbers]):
exp = EXPS[numbers][v]
if 0 < v <= 100 and 4*v == round(4*v):
expr(numbers, sqrt(v), '√' + exp)
if 2 <= v <= 6 and v == round(v):
expr(numbers, factorial(v), exp + '!')
if v != round(v):
uexp = unbracket(exp)
expr(numbers, floor(v), '⌊' + uexp + '⌋')
expr(numbers, ceil(v), '⌈' + uexp + '⌉')
def unbracket(exp):
"Remove outer brackets from exp if they are there."
if exp.startswith('(') and exp.endswith(')'):
return exp[1:-1]
else:
return exp
Let's try:
In [39]:
clear()
%time makeable(ff)
Out[39]:
The output got truncated by Jupyter Notebook after 1000 entries. Let's see what the first unmakeable integer actually is:
In [40]:
unmakeable(ff)
Out[40]: