#!/usr/bin/env python # coding: utf-8 # Peter Norvig
5 January 2016
revised 18 May 2018 # # Four 4s, Five 5s, and Countdown to 2016 # # On January 1, 2016 Alex Bellos [posed](http://www.theguardian.com/science/2016/jan/04/can-you-solve-it-complete-the-equation-10-9-8-7-6-5-4-3-2-1-2016) 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](http://www.theguardian.com/science/2016/jan/04/can-you-solve-it-complete-the-equation-10-9-8-7-6-5-4-3-2-1-2016) and [second](http://www.theguardian.com/science/2016/jan/04/did-you-solve-it-complete-the-equation-10-9-8-7-6-5-4-3-2-1-2016) 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 # In[3]: '10{}9{}8{}7{}6{}5{}4{}3{}2{}1'.format(*ops) # In[4]: eval(_) # 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) # 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)) # 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](https://oeis.org/A000108), 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](https://en.wikipedia.org/wiki/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) # 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,)) # In[12]: expressions((9, 8)) # In[13]: expressions((10, 9, 8)) # That looks reasonable. Let's solve the whole puzzle. # # # Countdown to 2016: A Solution # In[14]: get_ipython().run_line_magic('time', 'expressions(c10)[2016]') # 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)} # # 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'} # In[18]: counts((10, 9, 8))[27] # (10+9)+8 or 10+(9+8) # Looks good to me. Now let's see if we can answer the question. # In[19]: counts(c10)[2016] # 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} # 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()) # 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](https://en.wikipedia.org/wiki/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](https://archive.org/details/mathematicalrecr00ball) 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)} # 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)))`](http://www.wolframalpha.com/input/?i=10+%5E+9+%5E+56), 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](http://www.sympy.org/en/index.html), 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]: get_ipython().run_line_magic('time', 'show((4, 4, 4, 4))') # We can also solve the "2016 with four fours" puzzle: # In[29]: expressions((4, 4, 4, 4))[2016] # In a [separate video](https://www.youtube.com/embed/Noo4lN-vSvw), 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](https://nebusresearch.wordpress.com/2014/04/18/how-dirac-made-every-number/) for the "four 2s" problem). # # Donald Knuth has [conjectured](https://www.tandfonline.com/doi/abs/10.1080/0025570X.1964.11975546) 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)) # # Four 9s # In[31]: show((9, 9, 9, 9)) # # Four 5s # In[32]: show((5, 5, 5, 5)) # # Five 5s # In[33]: get_ipython().run_line_magic('time', 'show((5, 5, 5, 5, 5), False)') # # # Countdown to 2018 # # One more thing: On January 1 2018, [Michael Littman](http://cs.brown.edu/~mlittman/) 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) # In[35]: show((2,0,1,9), 10) # 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. # In[38]: show((2, 0, 2, 1), 10) # # 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](http://forums.xkcd.com/viewtopic.php?f=14&t=116813&start=280) they got up to 298 for five fives, using the double factorial and π functions. What would you like to do?