#!/usr/bin/env python
# coding: utf-8

# <div style="float:right"><i>Peter Norvig<br>5 January 2016</i></div>
# 
# # Countdown to 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 &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: +, -, &times;, ÷. But brackets (parentheses) can be used wherever needed. So, for example, the solution could begin
# 
# > (10 - 9) &times; (8 ...
# 
# > or
# 
# > 10 + (9 &times; 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.  We'll start with a  simpler version of the puzzle.

# # Four Operators, No Brackets
# 
# Suppose for the moment we are not allowed to use brackets.  Then there are nine blanks,  each of which can be filled by one of four operators, so the total number of possible expressions is:

# In[1]:


4 ** 9


# 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()


# 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]:


get_ipython().run_line_magic('time', 'solve_no_brackets(range(1900, 2100))')


# 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.

# # Four Operators, With Brackets
# 
# 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](https://en.wikipedia.org/wiki/Catalan_number), and we can [look it up](http://www.wolframalpha.com/input/?i=9th+catalan+number) to find that there are 4862 diffferent bracketing. If we enumerated and evaluated all of them, 4862 times half a minute is [40 hours](https://www.google.com/webhp?sourceid=chrome-instant&ion=1&espv=2&ie=UTF-8#q=4862+*+1%2F2+minute+in+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](https://en.wikipedia.org/wiki/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)


# 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))


# In[11]:


EXPS[(10, 9, 8)][27]


# That looks reasonable. Let's solve the whole puzzle.
# 
# # Countdown to 2016: A Solution

# In[12]:


get_ipython().run_line_magic('time', 'expressions(c10)[2016]')


# 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


# # Counting Solutions
# 
# 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 &times; 4<sup>9</sup> = 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 &times; 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))


# In[16]:


COUNTS[(10, 9, 8)]


# In[17]:


COUNTS[(10, 9, 8)][27]


# Looks good to me. Now let's repeat the computation, this time keeping track of `COUNTS`:

# In[18]:


clear()

get_ipython().run_line_magic('time', 'expressions(c10)[2016]')


# # The Answer (?)
# 
# Now we can read off the answer:

# In[19]:


COUNTS[c10][2016]


# 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:
# 
# - How many distinct values can be made?

# In[20]:


len(expressions(c10))


# - What's the smallest positive integer than cannot be made?

# 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)


# - What is the expected theoretical number of expressions? It is [Catalan(9)](http://www.wolframalpha.com/input/?i=9th+catalan+number) &times; 4<sup>9</sup>:

# In[22]:


4862 * 4 ** 9


# - How many expressions were actually entered in the `COUNTS` table?

# In[23]:


sum(COUNTS[c10].values())


# So, about 1% of the 1.27 billion theoretically possible expressions do not appear in `COUNTS`; these are the `ZeroDivisionErrors`.
# 
# # Dealing with Round-off Errors
# 
# Why don't I believe the answer of  30,066 expressions? Because floating point division can have round-off errors.  For example:

# In[24]:


2015 + 1/3 + 1/3 + 1/3


# In[25]:


2015 + 1/3 + 1/3 + 1/3 == 2016


# 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}


# 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


# # The Answer (!)

# In[28]:


clear()

get_ipython().run_line_magic('time', 'expressions(c10)[2016]')


# In[29]:


COUNTS[c10][2016]


# 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.

# # Dealing with the Exponentiation Operator
# 
# Now let's turn to another of Alex's puzzles: making 2016 and other target values from a string of four or five `4`s,  with exponentiation allowed. Exponentiation is tricky for five reasons:
# 
# - **Division by zero**: `(0 ** -1)` is the same as `(1 / 0)`, and gives a `ZeroDivisionError`.
# - **Irrationals**: `(2 ** (1 / 2))` is an irrational number; so we can't do exact rational arithmetic.
# - **Imaginaries**: `(-1 ** (1 / 2))` is an imaginary number, but Python gives a `ValueError`.
# - **Overflow**: `(10. ** (9. ** 8.))`, as a `float`, gives a `OverflowError`.
# - **Finite 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 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](http://www.sympy.org/en/index.html), 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]


# 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))


# 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)


# We didn't get there. 
# 
# # More Operations
# 
# With some research, I [see](http://www.infobarrel.com/Five_Fives_Problem_Recreational_Mathematics) that others who got up to 55 with five 5s used these three concepts:
# 
# - **digit concatenation**: `55`
# - **decimal point**: `.5`
# - **unary operations**: 
#   - **factorial**: `5! = 120` 
#   - **square root**:  &radic; `5`
# 
# 
# We'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()

get_ipython().run_line_magic('time', 'makeable(ff)')


# 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.

# # Even More Operations 
# 
# At the risk of making the computation take even longer, I'm going to add two more operations:
# 
# - **Floor**: &lfloor;*x*&rfloor; is the largest integer less than or equal to *x* (in other words, rounding down).
# - **Ceiling**: &lceil;*x*&rceil; is the smallest integer greater than or equal to *x* (in other words, rounding up).
# 
# 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
# 
# - ⌊√5⌋ = 2
# - ⌈5.5⌉! = 720
# 
# But still not allowing three consecutive applications, such as
# 
# - ⌈√5⌉! = 6
# 
# 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()

get_ipython().run_line_magic('time', 'makeable(ff)')


# The output got truncated by Jupyter Notebook after 1000 entries. Let's see what the first unmakeable integer actually is:

# In[40]:


unmakeable(ff)


# I think we've conquered the five 5s problem.
# 
# # Conclusion
# 
# We were able to solve the puzzles, and we now have some tools to explore new puzzles. What can you discover?