fifteen_factorial = 15*14*13*12*11*10*9*8*7*6*5*4*3*2*1
print(fifteen_factorial)

"""
Explanation: Suggestions for lab exercises.
Variables and assignment
Exercise 1
Remember that $n! = n \times (n - 1) \times \dots \times 2 \times 1$. Compute $15!$, assigning the result to a sensible variable name.
Solution
End of explanation
"""


import math
print(math.factorial(15))
print("Result correct?", math.factorial(15) == fifteen_factorial)

"""
Explanation: Exercise 2
Using the math module, check your result for $15$ factorial. You should explore the help for the math library and its functions, using eg tab-completion, the spyder inspector, or online sources.
Solution
End of explanation
"""


print(math.factorial(5), math.sqrt(2*math.pi)*5**(5+0.5)*math.exp(-5))
print(math.factorial(10), math.sqrt(2*math.pi)*10**(10+0.5)*math.exp(-10))
print(math.factorial(15), math.sqrt(2*math.pi)*15**(15+0.5)*math.exp(-15))
print(math.factorial(20), math.sqrt(2*math.pi)*20**(20+0.5)*math.exp(-20))
print("Absolute differences:")
print(math.factorial(5) - math.sqrt(2*math.pi)*5**(5+0.5)*math.exp(-5))
print(math.factorial(10) - math.sqrt(2*math.pi)*10**(10+0.5)*math.exp(-10))
print(math.factorial(15) - math.sqrt(2*math.pi)*15**(15+0.5)*math.exp(-15))
print(math.factorial(20) - math.sqrt(2*math.pi)*20**(20+0.5)*math.exp(-20))
print("Relative differences:")
print((math.factorial(5) - math.sqrt(2*math.pi)*5**(5+0.5)*math.exp(-5)) / math.factorial(5))
print((math.factorial(10) - math.sqrt(2*math.pi)*10**(10+0.5)*math.exp(-10)) / math.factorial(10))
print((math.factorial(15) - math.sqrt(2*math.pi)*15**(15+0.5)*math.exp(-15)) / math.factorial(15))
print((math.factorial(20) - math.sqrt(2*math.pi)*20**(20+0.5)*math.exp(-20)) / math.factorial(20))

"""
Explanation: Exercise 3
Stirling's approximation gives that, for large enough $n$, 
\begin{equation}
  n! \simeq \sqrt{2 \pi} n^{n + 1/2} e^{-n}.
\end{equation}
Using functions and constants from the math library, compare the results of $n!$ and Stirling's approximation for $n = 5, 10, 15, 20$. In what sense does the approximation improve?
Solution
End of explanation
"""


def cuboid_volume(a, b, c):
    """
    Compute the volume of a cuboid with edge lengths a, b, c.
    Volume is abc. Only makes sense if all are non-negative.
    
    Parameters
    ----------
    
    a : float
        Edge length 1
    b : float
        Edge length 2
    c : float
        Edge length 3
        
    Returns
    -------
    
    volume : float
        The volume a*b*c
    """
    
    if (a < 0.0) or (b < 0.0) or (c < 0.0):
        print("Negative edge length makes no sense!")
        return 0
    
    return a*b*c

print(cuboid_volume(1,1,1))
print(cuboid_volume(1,2,3.5))
print(cuboid_volume(0,1,1))
print(cuboid_volume(2,-1,1))

"""
Explanation: We see that the relative error decreases, whilst the absolute error grows (significantly).
Basic functions
Exercise 1
Write a function to calculate the volume of a cuboid with edge lengths $a, b, c$. Test your code on sample values such as

$a=1, b=1, c=1$ (result should be $1$);
$a=1, b=2, c=3.5$ (result should be $7.0$);
$a=0, b=1, c=1$ (result should be $0$);
$a=2, b=-1, c=1$ (what do you think the result should be?).

Solution
End of explanation
"""


def fall_time(H):
    """
    Give the time in seconds for an object to fall to the ground
    from H metres.
    
    Parameters
    ----------
    
    H : float
        Starting height (metres)
        
    Returns
    -------
    
    T : float
        Fall time (seconds)
    """
    
    from math import sqrt
    from scipy.constants import g
    
    if (H < 0):
        print("Negative height makes no sense!")
        return 0
    
    return sqrt(2.0*H/g)

print(fall_time(1))
print(fall_time(10))
print(fall_time(0))
print(fall_time(-1))

"""
Explanation: In later cases, after having covered exceptions, I would suggest raising a NotImplementedError for negative edge lengths.
Exercise 2
Write a function to compute the time (in seconds) taken for an object to fall from a height $H$ (in metres) to the ground, using the formula
\begin{equation}
  h(t) = \frac{1}{2} g t^2.
\end{equation}
Use the value of the acceleration due to gravity $g$ from scipy.constants.g. Test your code on sample values such as

$H = 1$m (result should be $\approx 0.452$s);
$H = 10$m (result should be $\approx 1.428$s);
$H = 0$m (result should be $0$s);
$H = -1$m (what do you think the result should be?).

Solution
End of explanation
"""


def triangle_area(a, b, c):
    """
    Compute the area of a triangle with edge lengths a, b, c.
    Area is sqrt(s (s-a) (s-b) (s-c)). 
    s is (a+b+c)/2.
    Only makes sense if all are non-negative.
    
    Parameters
    ----------
    
    a : float
        Edge length 1
    b : float
        Edge length 2
    c : float
        Edge length 3
        
    Returns
    -------
    
    area : float
        The triangle area.
    """
    
    from math import sqrt
    
    if (a < 0.0) or (b < 0.0) or (c < 0.0):
        print("Negative edge length makes no sense!")
        return 0
    
    s = 0.5 * (a + b + c)
    return sqrt(s * (s-a) * (s-b) * (s-c))

print(triangle_area(1,1,1)) # Equilateral; answer sqrt(3)/4 ~ 0.433
print(triangle_area(3,4,5)) # Right triangle; answer 6
print(triangle_area(1,1,0)) # Not a triangle; answer 0
print(triangle_area(-1,1,1)) # Not a triangle; exception or 0.

"""
Explanation: Exercise 3
Write a function that computes the area of a triangle with edge lengths $a, b, c$. You may use the formula
\begin{equation}
  A = \sqrt{s (s - a) (s - b) (s - c)}, \qquad s = \frac{a + b + c}{2}.
\end{equation}
Construct your own test cases to cover a range of possibilities.
End of explanation
"""


from math import sqrt

x = 1.0
y = 1.0 + 1e-14 * sqrt(3.0)
print("The calculation gives {}".format(1e14*(y-x)))
print("The result should be {}".format(sqrt(3.0)))

"""
Explanation: Floating point numbers
Exercise 1
Computers cannot, in principle, represent real numbers perfectly. This can lead to problems of accuracy. For example, if
\begin{equation}
  x = 1, \qquad y = 1 + 10^{-14} \sqrt{3}
\end{equation}
then it should be true that
\begin{equation}
  10^{14} (y - x) = \sqrt{3}.
\end{equation}
Check how accurately this equation holds in Python and see what this implies about the accuracy of subtracting two numbers that are close together.
Solution
End of explanation
"""


a = 1e-3
b = 1e3
c = a
formula1_n3_plus = (-b + sqrt(b**2 - 4.0*a*c))/(2.0*a)
formula1_n3_minus = (-b - sqrt(b**2 - 4.0*a*c))/(2.0*a)
formula2_n3_plus = (2.0*c)/(-b + sqrt(b**2 - 4.0*a*c))
formula2_n3_minus = (2.0*c)/(-b - sqrt(b**2 - 4.0*a*c))
print("For n=3, first formula, solutions are {} and {}.".format(formula1_n3_plus, 
                                                                formula1_n3_minus))
print("For n=3, second formula, solutions are {} and {}.".format(formula2_n3_plus, 
                                                                 formula2_n3_minus))

a = 1e-4
b = 1e4
c = a
formula1_n4_plus = (-b + sqrt(b**2 - 4.0*a*c))/(2.0*a)
formula1_n4_minus = (-b - sqrt(b**2 - 4.0*a*c))/(2.0*a)
formula2_n4_plus = (2.0*c)/(-b + sqrt(b**2 - 4.0*a*c))
formula2_n4_minus = (2.0*c)/(-b - sqrt(b**2 - 4.0*a*c))
print("For n=4, first formula, solutions are {} and {}.".format(formula1_n4_plus, 
                                                                formula1_n4_minus))
print("For n=4, second formula, solutions are {} and {}.".format(formula2_n4_plus, 
                                                                 formula2_n4_minus))

"""
Explanation: We see that the first three digits are correct. This isn't too surprising: we expect 16 digits of accuracy for a floating point number, but $x$ and $y$ are identical for the first 14 digits.
Exercise 2
The standard quadratic formula gives the solutions to
\begin{equation}
  a x^2 + b x + c = 0
\end{equation}
as
\begin{equation}
  x = \frac{-b \pm \sqrt{b^2 - 4 a c}}{2 a}.
\end{equation}
Show that, if $a = 10^{-n} = c$ and $b = 10^n$ then
\begin{equation}
  x = \frac{10^{2 n}}{2} \left( -1 \pm \sqrt{1 - 10^{-4n}} \right).
\end{equation}
Using the expansion (from Taylor's theorem)
\begin{equation}
  \sqrt{1 - 10^{-4 n}} \simeq 1 - \frac{10^{-4 n}}{2} + \dots, \qquad n \gg 1,
\end{equation}
show that
\begin{equation}
  x \simeq -10^{2 n} + \frac{10^{-2 n}}{4} \quad \text{and} \quad -\frac{10^{-2n}}{4}, \qquad n \gg 1.
\end{equation}
Solution
This is pen-and-paper work; each step should be re-arranging.
Exercise 3
By multiplying and dividing by $-b \mp \sqrt{b^2 - 4 a c}$, check that we can also write the solutions to the quadratic equation as
\begin{equation}
  x = \frac{2 c}{-b \mp \sqrt{b^2 - 4 a c}}.
\end{equation}
Solution
Using the difference of two squares we get
\begin{equation}
  x = \frac{b^2 - \left( b^2 - 4 a c \right)}{2a \left( -b \mp \sqrt{b^2 - 4 a c} \right)}
\end{equation}
which re-arranges to give the required solution.
Exercise 4
Using Python, calculate both solutions to the quadratic equation
\begin{equation}
  10^{-n} x^2 + 10^n x + 10^{-n} = 0
\end{equation}
for $n = 3$ and $n = 4$ using both formulas. What do you see? How has floating point accuracy caused problems here?
Solution
End of explanation
"""


def g(f, X, delta):
    """
    Approximate the derivative of a given function at a point.
    
    Parameters
    ----------
    
    f : function
        Function to be differentiated
    X : real
        Point at which the derivative is evaluated
    delta : real
        Step length
        
    Returns
    -------
    
    g : real
        Approximation to the derivative
    """
    
    return (f(X+delta) - f(X)) / delta

"""
Explanation: There is a difference in the fifth significant figure in both solutions in the first case, which gets to the third (arguably the second) significant figure in the second case. Comparing to the limiting solutions above, we see that the larger root is definitely more accurately captured with the first formula than the second (as the result should be bigger than $10^{-2n}$).
In the second case we have divided by a very small number to get the big number, which loses accuracy.
Exercise 5
The standard definition of the derivative of a function is
\begin{equation}
  \left. \frac{\text{d} f}{\text{d} x} \right|{x=X} = \lim{\delta \to 0} \frac{f(X + \delta) - f(X)}{\delta}.
\end{equation}
We can approximate this by computing the result for a finite value of $\delta$:
\begin{equation}
  g(x, \delta) = \frac{f(x + \delta) - f(x)}{\delta}.
\end{equation}
Write a function that takes as inputs a function of one variable, $f(x)$, a location $X$, and a step length $\delta$, and returns the approximation to the derivative given by $g$.
Solution
End of explanation
"""


from math import exp
for n in range(1, 8):
    print("For n={}, the approx derivative is {}.".format(n, g(exp, 0.0, 10**(-2.0*n))))

"""
Explanation: Exercise 6
The function $f_1(x) = e^x$ has derivative with the exact value $1$ at $x=0$. Compute the approximate derivative using your function above, for $\delta = 10^{-2 n}$ with $n = 1, \dots, 7$. You should see the results initially improve, then get worse. Why is this?
Solution
End of explanation
"""


def isprime(n):
    """
    Checks to see if an integer is prime.
    
    Parameters
    ----------
    
    n : integer
        Number to check
        
    Returns
    -------
    
    isprime : Boolean
        If n is prime
    """
    
    # No number less than 2 can be prime
    if n < 2:
        return False
    
    # We only need to check for divisors up to sqrt(n)
    for m in range(2, int(n**0.5)+1):
        if n%m == 0:
            return False
    
    # If we've got this far, there are no divisors.
    return True

for n in range(50):
    if isprime(n):
        print("Function says that {} is prime.".format(n))

"""
Explanation: We have a combination of floating point inaccuracies: in the numerator we have two terms that are nearly equal, leading to a very small number. We then divide two very small numbers. This is inherently inaccurate.
This does not mean that you can't calculate derivatives to high accuracy, but alternative approaches are definitely recommended.
Prime numbers
Exercise 1
Write a function that tests if a number is prime. Test it by writing out all prime numbers less than 50.
Solution
This is a "simple" solution, but not efficient.
End of explanation
"""


n = 2
while (not isprime(n)) or (isprime(2**n-1)):
    n += 1
print("The first n such that 2^n-1 is not prime is {}.".format(n))

"""
Explanation: Exercise 2
500 years ago some believed that the number $2^n - 1$ was prime for all primes $n$. Use your function to find the first prime $n$ for which this is not true.
Solution
We could do this many ways. This "elegant" solution says:

Start from the smallest possible $n$ (2).
Check if $n$ is prime. If not, add one to $n$.
If $n$ is prime, check if $2^n-1$ is prime. If it is, add one to $n$.
If both those logical checks fail, we have found the $n$ we want.
End of explanation
"""


for n in range(2, 41):
    if isprime(n) and isprime(2**n-1):
        print("n={} is such that 2^n-1 is prime.".format(n))

"""
Explanation: Exercise 3
The Mersenne primes are those that have the form $2^n-1$, where $n$ is prime. Use your previous solutions to generate all the $n < 40$ that give Mersenne primes.
Solution
End of explanation
"""


def prime_factors(n):
    """
    Generate all the prime factors of n.
    
    Parameters
    ----------
    
    n : integer
        Number to be checked
        
    Returns
    -------
    
    factors : dict
        Prime factors (keys) and multiplicities (values)
    """
    
    factors = {}
    
    m = 2
    while m <= n:
        if n%m == 0:
            factors[m] = 1
            n //= m
            while n%m == 0:
                factors[m] += 1
                n //= m
        m += 1
        
    return factors

for n in range(17, 21):
    print("Prime factors of {} are {}.".format(n, prime_factors(n).keys()))
print("Multiplicities of prime factors of 48 are {}.".format(prime_factors(48).values()))

"""
Explanation: Exercise 4
Write a function to compute all prime factors of an integer $n$, including their multiplicities. Test it by printing the prime factors (without multiplicities) of $n = 17, \dots, 20$ and the multiplicities (without factors) of $n = 48$.
Note
One effective solution is to return a dictionary, where the keys are the factors and the values are the multiplicities.
Solution
This solution uses the trick of immediately dividing $n$ by any divisor: this means we never have to check the divisor for being prime.
End of explanation
"""


def divisors(n):
    """
    Generate all integer divisors of n.
    
    Parameters
    ----------
    
    n : integer
        Number to be checked
        
    Returns
    -------
    
    divs : list
        All integer divisors, including 1.
    """
    
    divs = [1]
    m = 2
    while m <= n/2:
        if n%m == 0:
            divs.append(m)
        m += 1
        
    return divs

for n in range(16, 21):
    print("The divisors of {} are {}.".format(n, divisors(n)))

"""
Explanation: Exercise 5
Write a function to generate all the integer divisors, including 1, but not including $n$ itself, of an integer $n$. Test it on $n = 16, \dots, 20$.
Note
You could use the prime factorization from the previous exercise, or you could do it directly.
Solution
Here we will do it directly.
End of explanation
"""


def isperfect(n):
    """
    Check if a number is perfect.
    
    Parameters
    ----------
    
    n : integer
        Number to check
        
    Returns
    -------
    
    isperfect : Boolean
        Whether it is perfect or not.
    """
    
    divs = divisors(n)
    sum_divs = 0
    for d in divs:
        sum_divs += d
    
    return n == sum_divs

for n in range(2,10000):
    if (isperfect(n)):
        factors = prime_factors(n)
        print("{} is perfect.\n"
              "Divisors are {}.\n"
              "Prime factors {} (multiplicities {}).".format(
            n, divisors(n), factors.keys(), factors.values()))

"""
Explanation: Exercise 6
A perfect number $n$ is one where the divisors sum to $n$. For example, 6 has divisors 1, 2, and 3, which sum to 6. Use your previous solution to find all perfect numbers $n < 10,000$ (there are only four!).
Solution
We can do this much more efficiently than the code below using packages such as numpy, but this is a "bare python" solution.
End of explanation
"""


%timeit isperfect(2**(3-1)*(2**3-1))

%timeit isperfect(2**(5-1)*(2**5-1))

%timeit isperfect(2**(7-1)*(2**7-1))

%timeit isperfect(2**(13-1)*(2**13-1))

"""
Explanation: Exercise 7
Using your previous functions, check that all perfect numbers $n < 10,000$ can be written as $2^{k-1} \times (2^k - 1)$, where $2^k-1$ is a Mersenne prime.
Solution
In fact we did this above already:

$6 = 2^{2-1} \times (2^2 - 1)$. 2 is the first number on our Mersenne list.
$28 = 2^{3-1} \times (2^3 - 1)$. 3 is the second number on our Mersenne list.
$496 = 2^{5-1} \times (2^5 - 1)$. 5 is the third number on our Mersenne list.
$8128 = 2^{7-1} \times (2^7 - 1)$. 7 is the fourth number on our Mersenne list.

Exercise 8 (bonus)
Investigate the timeit function in python or IPython. Use this to measure how long your function takes to check that, if $k$ on the Mersenne list then $n = 2^{k-1} \times (2^k - 1)$ is a perfect number, using your functions. Stop increasing $k$ when the time takes too long!
Note
You could waste considerable time on this, and on optimizing the functions above to work efficiently. It is not worth it, other than to show how rapidly the computation time can grow!
Solution
End of explanation
"""


def logistic(x0, r, N = 1000):
    sequence = [x0]
    xn = x0
    for n in range(N):
        xnew = r*xn*(1.0-xn)
        sequence.append(xnew)
        xn = xnew
    return sequence

"""
Explanation: It's worth thinking about the operation counts of the various functions implemented here. The implementations are inefficient, but even in the best case you see how the number of operations (and hence computing time required) rapidly increases.
Logistic map
Partly taken from Newman's book, p 120.
The logistic map builds a sequence of numbers ${ x_n }$ using the relation
\begin{equation}
  x_{n+1} = r x_n \left( 1 - x_n \right),
\end{equation}
where $0 \le x_0 \le 1$.
Exercise 1
Write a program that calculates the first $N$ members of the sequence, given as input $x_0$ and $r$ (and, of course, $N$).
Solution
End of explanation
"""


import numpy
from matplotlib import pyplot
%matplotlib inline

x0 = 0.5
N = 2000
sequence1 = logistic(x0, 1.5, N)
sequence2 = logistic(x0, 3.5, N)
pyplot.plot(sequence1[-100:], 'b-', label = r'$r=1.5$')
pyplot.plot(sequence2[-100:], 'k-', label = r'$r=3.5$')
pyplot.xlabel(r'$n$')
pyplot.ylabel(r'$x$')
pyplot.show()

"""
Explanation: Exercise 2
Fix $x_0=0.5$. Calculate the first 2,000 members of the sequence for $r=1.5$ and $r=3.5$ Plot the last 100 members of the sequence in both cases.
What does this suggest about the long-term behaviour of the sequence?
Solution
End of explanation
"""


import numpy
from matplotlib import pyplot
%matplotlib inline

r_values = numpy.linspace(1.0, 4.0, 401)
x0 = 0.5
N = 2000
for r in r_values:
    sequence = logistic(x0, r, N)
    pyplot.plot(r*numpy.ones_like(sequence[1000:]), sequence[1000:], 'k.')
pyplot.xlabel(r'$r$')
pyplot.ylabel(r'$x$')
pyplot.show()

"""
Explanation: This suggests that, for $r=1.5$, the sequence has settled down to a fixed point. In the $r=3.5$ case it seems to be moving between four points repeatedly.
Exercise 3
Fix $x_0 = 0.5$. For each value of $r$ between $1$ and $4$, in steps of $0.01$, calculate the first 2,000 members of the sequence. Plot the last 1,000 members of the sequence on a plot where the $x$-axis is the value of $r$ and the $y$-axis is the values in the sequence. Do not plot lines - just plot markers (e.g., use the 'k.' plotting style).
Solution
End of explanation
"""


def in_Mandelbrot(c, n_iterations = 100):
    z0 = 0.0 + 0j
    in_set = True
    n = 0
    zn = z0
    while in_set and (n < n_iterations):
        n += 1
        znew = zn**2 + c
        in_set = abs(znew) < 2.0
        zn = znew
    return in_set

"""
Explanation: Exercise 4
For iterative maps such as the logistic map, one of three things can occur:

The sequence settles down to a fixed point.
The sequence rotates through a finite number of values. This is called a limit cycle.
The sequence generates an infinite number of values. This is called deterministic chaos.

Using just your plot, or new plots from this data, work out approximate values of $r$ for which there is a transition from fixed points to limit cycles, from limit cycles of a given number of values to more values, and the transition to chaos.
Solution
The first transition is at $r \approx 3$, the next at $r \approx 3.45$, the next at $r \approx 3.55$. The transition to chaos appears to happen before $r=4$, but it's not obvious exactly where.
Mandelbrot
The Mandelbrot set is also generated from a sequence, ${ z_n }$, using the relation
\begin{equation}
  z_{n+1} = z_n^2 + c, \qquad z_0 = 0.
\end{equation}
The members of the sequence, and the constant $c$, are all complex. The point in the complex plane at $c$ is in the Mandelbrot set only if the $|z_n| < 2$ for all members of the sequence. In reality, checking the first 100 iterations is sufficient.
Note: the python notation for a complex number $x + \text{i} y$ is x + yj: that is, j is used to indicate $\sqrt{-1}$. If you know the values of x and y then x + yj constructs a complex number; if they are stored in variables you can use complex(x, y).
Exercise 1
Write a function that checks if the point $c$ is in the Mandelbrot set.
Solution
End of explanation
"""


c_values = [0.0, 2+2j, 2-2j, -2+2j, -2-2j]
for c in c_values:
    print("Is {} in the Mandelbrot set? {}.".format(c, in_Mandelbrot(c)))

"""
Explanation: Exercise 2
Check the points $c=0$ and $c=\pm 2 \pm 2 \text{i}$ and ensure they do what you expect. (What should you expect?)
Solution
End of explanation
"""


import numpy

def grid_Mandelbrot(N):
    x = numpy.linspace(-2.0, 2.0, N)
    X, Y = numpy.meshgrid(x, x)
    C = X + 1j*Y
    grid = numpy.zeros((N, N), int)
    for nx in range(N):
        for ny in range(N):
            grid[nx, ny] = int(in_Mandelbrot(C[nx, ny]))
    return grid

"""
Explanation: Exercise 3
Write a function that, given $N$

generates an $N \times N$ grid spanning $c = x + \text{i} y$, for $-2 \le x \le 2$ and $-2 \le y \le 2$;
returns an $N\times N$ array containing one if the associated grid point is in the Mandelbrot set, and zero otherwise.

Solution
End of explanation
"""


from matplotlib import pyplot
%matplotlib inline

pyplot.imshow(grid_Mandelbrot(100))

"""
Explanation: Exercise 4
Using the function imshow from matplotlib, plot the resulting array for a $100 \times 100$ array to make sure you see the expected shape.
Solution
End of explanation
"""


from math import log

def log_Mandelbrot(c, n_iterations = 100):
    z0 = 0.0 + 0j
    in_set = True
    n = 0
    zn = z0
    while in_set and (n < n_iterations):
        n += 1
        znew = zn**2 + c
        in_set = abs(znew) < 2.0
        zn = znew
    return log(n)

def log_grid_Mandelbrot(N):
    x = numpy.linspace(-2.0, 2.0, N)
    X, Y = numpy.meshgrid(x, x)
    C = X + 1j*Y
    grid = numpy.zeros((N, N), int)
    for nx in range(N):
        for ny in range(N):
            grid[nx, ny] = log_Mandelbrot(C[nx, ny])
    return grid

from matplotlib import pyplot
%matplotlib inline

pyplot.imshow(log_grid_Mandelbrot(100))

"""
Explanation: Exercise 5
Modify your functions so that, instead of returning whether a point is inside the set or not, it returns the logarithm of the number of iterations it takes. Plot the result using imshow again.
Solution
End of explanation
"""


pyplot.imshow(log_grid_Mandelbrot(1000)[600:800,400:600])

"""
Explanation: Exercise 6
Try some higher resolution plots, and try plotting only a section to see the structure. Note this is not a good way to get high accuracy close up images!
Solution
This is a simple example:
End of explanation
"""


class Eqint(object):
    
    def __init__(self, sequence):
        self.sequence = sequence
        
    def __repr__(self):
        return str(len(self.sequence))
    
    def __eq__(self, other):
        return len(self.sequence)==len(other.sequence)

"""
Explanation: Equivalence classes
An equivalence class is a relation that groups objects in a set into related subsets. For example, if we think of the integers modulo $7$, then $1$ is in the same equivalence class as $8$ (and $15$, and $22$, and so on), and $3$ is in the same equivalence class as $10$. We use the tilde $3 \sim 10$ to denote two objects within the same equivalence class.
Here, we are going to define the positive integers programmatically from equivalent sequences.
Exercise 1
Define a python class Eqint. This should be

Initialized by a sequence;
Store the sequence;
Define its representation (via the __repr__ function) to be the integer length of the sequence;
Redefine equality (via the __eq__ function) so that two eqints are equal if their sequences have same length.

Solution
End of explanation
"""


zero = Eqint([])
one_list = Eqint([1])
one_tuple = Eqint((1,))
one_string = Eqint('1')

print("Is zero equivalent to one? {}, {}, {}".format(zero == one_list, 
                                                     zero == one_tuple,
                                                     zero == one_string))
print("Is one equivalent to one? {}, {}, {}.".format(one_list == one_tuple,
                                                     one_list == one_string,
                                                     one_tuple == one_string))
print(zero)
print(one_list)
print(one_tuple)
print(one_string)

"""
Explanation: Exercise 2
Define a zero object from the empty list, and three one objects, from a single object list, tuple, and string. For example
python
one_list = Eqint([1])
one_tuple = Eqint((1,))
one_string = Eqint('1')
Check that none of the one objects equal the zero object, but all equal the other one objects. Print each object to check that the representation gives the integer length.
Solution
End of explanation
"""


class Eqint(object):
    
    def __init__(self, sequence):
        self.sequence = sequence
        
    def __repr__(self):
        return str(len(self.sequence))
    
    def __eq__(self, other):
        return len(self.sequence)==len(other.sequence)
    
    def __add__(a, b):
        return Eqint(tuple(a.sequence) + tuple(b.sequence))

"""
Explanation: Exercise 3
Redefine the class by including an __add__ method that combines the two sequences. That is, if a and b are Eqints then a+b should return an Eqint defined from combining a and bs sequences.
Note
Adding two different types of sequences (eg, a list to a tuple) does not work, so it is better to either iterate over the sequences, or to convert to a uniform type before adding.
Solution
End of explanation
"""


zero = Eqint([])
one_list = Eqint([1])
one_tuple = Eqint((1,))
one_string = Eqint('1')

sum_eqint = zero + one_list + one_tuple + one_string
print("The sum is {}.".format(sum_eqint))
print("The internal sequence is {}.".format(sum_eqint.sequence))

"""
Explanation: Exercise 4
Check your addition function by adding together all your previous Eqint objects (which will need re-defining, as the class has been redefined). Print the resulting object to check you get 3, and also print its internal sequence.
Solution
End of explanation
"""


positive_integers = []
zero = Eqint([])
positive_integers.append(zero)

N = 10
for n in range(1,N+1):
    positive_integers.append(Eqint(list(positive_integers)))
    
print("The 'final' Eqint is {}".format(positive_integers[-1]))
print("Its sequence is {}".format(positive_integers[-1].sequence))
print("That is, it contains all Eqints with length less than 10.")

"""
Explanation: Exercise 5
We will sketch a construction of the positive integers from nothing.

Define an empty list positive_integers.
Define an Eqint called zero from the empty list. Append it to positive_integers.
Define an Eqint called next_integer from the Eqint defined by a copy of positive_integers (ie, use Eqint(list(positive_integers)). Append it to positive_integers.
Repeat step 3 as often as needed.

Use this procedure to define the Eqint equivalent to $10$. Print it, and its internal sequence, to check.
Solution
End of explanation
"""


def normal_form(numerator, denominator):
    from fractions import gcd
    
    factor = gcd(numerator, denominator)
    return numerator//factor, denominator//factor

print(normal_form(3, 2))
print(normal_form(15, 3))
print(normal_form(20, 42))

"""
Explanation: Rational numbers
Instead of working with floating point numbers, which are not "exact", we could work with the rational numbers $\mathbb{Q}$. A rational number $q \in \mathbb{Q}$ is defined by the numerator $n$ and denominator $d$ as $q = \frac{n}{d}$, where $n$ and $d$ are coprime (ie, have no common divisor other than $1$).
Exercise 1
Find a python function that finds the greatest common divisor (gcd) of two numbers. Use this to write a function normal_form that takes a numerator and divisor and returns the coprime $n$ and $d$. Test this function on $q = \frac{3}{2}$, $q = \frac{15}{3}$, and $q = \frac{20}{42}$.
Solution
End of explanation
"""


class Rational(object):
    """
    A rational number.
    """
    
    def __init__(self, numerator, denominator):
        
        n, d = normal_form(numerator, denominator)
        
        self.numerator = n
        self.denominator = d
        
        return None
    
    def __repr__(self):
        
        max_length = max(len(str(self.numerator)), len(str(self.denominator)))
        
        if self.denominator == 1:
            frac = str(self.numerator)
        else:
            numerator = str(self.numerator)+'\n'
            bar = max_length*'-'+'\n'
            denominator = str(self.denominator)
            frac = numerator+bar+denominator
        
        return frac

q1 = Rational(3, 2)
print(q1)
q2 = Rational(15, 3)
print(q2)
q3 = Rational(20, 42)
print(q3)

"""
Explanation: Exercise 2
Define a class Rational that uses the normal_form function to store the rational number in the appropriate form. Define a __repr__ function that prints a string that looks like $\frac{n}{d}$ (hint: use len(str(number)) to find the number of digits of an integer). Test it on the cases above.
Solution
End of explanation
"""


class Rational(object):
    """
    A rational number.
    """
    
    def __init__(self, numerator, denominator):
        
        n, d = normal_form(numerator, denominator)
        
        self.numerator = n
        self.denominator = d
        
        return None
    
    def __add__(a, b):
        
        numerator = a.numerator * b.denominator + b.numerator * a.denominator
        denominator = a.denominator * b.denominator
        return Rational(numerator, denominator)
    
    def __repr__(self):
        
        max_length = max(len(str(self.numerator)), len(str(self.denominator)))
        
        if self.denominator == 1:
            frac = str(self.numerator)
        else:
            numerator = str(self.numerator)+'\n'
            bar = max_length*'-'+'\n'
            denominator = str(self.denominator)
            frac = numerator+bar+denominator
        
        return frac

print(Rational(1,2) + Rational(1,3) + Rational(1,6))

"""
Explanation: Exercise 3
Overload the __add__ function so that you can add two rational numbers. Test it on $\frac{1}{2} + \frac{1}{3} + \frac{1}{6} = 1$.
Solution
End of explanation
"""


class Rational(object):
    """
    A rational number.
    """
    
    def __init__(self, numerator, denominator):
        
        n, d = normal_form(numerator, denominator)
        
        self.numerator = n
        self.denominator = d
        
        return None
    
    def __add__(a, b):
        
        numerator = a.numerator * b.denominator + b.numerator * a.denominator
        denominator = a.denominator * b.denominator
        return Rational(numerator, denominator)
    
    def __mul__(a, b):
        
        numerator = a.numerator * b.numerator
        denominator = a.denominator * b.denominator
        return Rational(numerator, denominator)
    
    def __repr__(self):
        
        max_length = max(len(str(self.numerator)), len(str(self.denominator)))
        
        if self.denominator == 1:
            frac = str(self.numerator)
        else:
            numerator = str(self.numerator)+'\n'
            bar = max_length*'-'+'\n'
            denominator = str(self.denominator)
            frac = numerator+bar+denominator
        
        return frac

print(Rational(1,3)*Rational(15,2)*Rational(2,5))

"""
Explanation: Exercise 4
Overload the __mul__ function so that you can multiply two rational numbers. Test it on $\frac{1}{3} \times \frac{15}{2} \times \frac{2}{5} = 1$.
Solution
End of explanation
"""


class Rational(object):
    """
    A rational number.
    """
    
    def __init__(self, numerator, denominator):
        
        n, d = normal_form(numerator, denominator)
        
        self.numerator = n
        self.denominator = d
        
        return None
    
    def __add__(a, b):
        
        numerator = a.numerator * b.denominator + b.numerator * a.denominator
        denominator = a.denominator * b.denominator
        return Rational(numerator, denominator)
    
    def __mul__(a, b):
        
        numerator = a.numerator * b.numerator
        denominator = a.denominator * b.denominator
        return Rational(numerator, denominator)
    
    def __rmul__(self, other):
        
        numerator = self.numerator * other
        return Rational(numerator, self.denominator)
    
    def __sub__(a, b):
        
        return a + (-1)*b
    
    def __repr__(self):
        
        max_length = max(len(str(self.numerator)), len(str(self.denominator)))
        
        if self.denominator == 1:
            frac = str(self.numerator)
        else:
            numerator = str(self.numerator)+'\n'
            bar = max_length*'-'+'\n'
            denominator = str(self.denominator)
            frac = numerator+bar+denominator
        
        return frac

half = Rational(1,2)
print(2*half)
print(half+(-1)*half)
print(half-half)

"""
Explanation: Exercise 5
Overload the __rmul__ function so that you can multiply a rational by an integer. Check that $\frac{1}{2} \times 2 = 1$ and $\frac{1}{2} + (-1) \times \frac{1}{2} = 0$. Also overload the __sub__ function (using previous functions!) so that you can subtract rational numbers and check that $\frac{1}{2} - \frac{1}{2} = 0$.
Solution
End of explanation
"""


class Rational(object):
    """
    A rational number.
    """
    
    def __init__(self, numerator, denominator):
        
        n, d = normal_form(numerator, denominator)
        
        self.numerator = n
        self.denominator = d
        
        return None
    
    def __add__(a, b):
        
        numerator = a.numerator * b.denominator + b.numerator * a.denominator
        denominator = a.denominator * b.denominator
        return Rational(numerator, denominator)
    
    def __mul__(a, b):
        
        numerator = a.numerator * b.numerator
        denominator = a.denominator * b.denominator
        return Rational(numerator, denominator)
    
    def __rmul__(self, other):
        
        numerator = self.numerator * other
        return Rational(numerator, self.denominator)
    
    def __sub__(a, b):
        
        return a + (-1)*b
    
    def __float__(a):
        
        return float(a.numerator) / float(a.denominator)
    
    def __repr__(self):
        
        max_length = max(len(str(self.numerator)), len(str(self.denominator)))
        
        if self.denominator == 1:
            frac = str(self.numerator)
        else:
            numerator = str(self.numerator)+'\n'
            bar = max_length*'-'+'\n'
            denominator = str(self.denominator)
            frac = numerator+bar+denominator
        
        return frac

print(float(Rational(1,2)))
print(float(Rational(1,3)))
print(float(Rational(1,11)))

"""
Explanation: Exercise 6
Overload the __float__ function so that float(q) returns the floating point approximation to the rational number q. Test this on $\frac{1}{2}, \frac{1}{3}$, and $\frac{1}{11}$.
Solution
End of explanation
"""


class Rational(object):
    """
    A rational number.
    """
    
    def __init__(self, numerator, denominator):
        
        n, d = normal_form(numerator, denominator)
        
        self.numerator = n
        self.denominator = d
        
        return None
    
    def __add__(a, b):
        
        numerator = a.numerator * b.denominator + b.numerator * a.denominator
        denominator = a.denominator * b.denominator
        return Rational(numerator, denominator)
    
    def __mul__(a, b):
        
        numerator = a.numerator * b.numerator
        denominator = a.denominator * b.denominator
        return Rational(numerator, denominator)
    
    def __rmul__(self, other):
        
        numerator = self.numerator * other
        return Rational(numerator, self.denominator)
    
    def __sub__(a, b):
        
        return a + (-1)*b
    
    def __float__(a):
        
        return float(a.numerator) / float(a.denominator)
    
    def __lt__(a, b):
        
        return a.numerator * b.denominator < a.denominator * b.numerator
    
    def __repr__(self):
        
        max_length = max(len(str(self.numerator)), len(str(self.denominator)))
        
        if self.denominator == 1:
            frac = str(self.numerator)
        else:
            numerator = '\n'+str(self.numerator)+'\n'
            bar = max_length*'-'+'\n'
            denominator = str(self.denominator)
            frac = numerator+bar+denominator
        
        return frac

q_list = [Rational(n//2, n) for n in range(2, 12)]
print(sorted(q_list))

"""
Explanation: Exercise 7
Overload the __lt__ function to compare two rational numbers. Create a list of rational numbers where the denominator is $n = 2, \dots, 11$ and the numerator is the floored integer $n/2$, ie n//2. Use the sorted function on that list (which relies on the __lt__ function).
Solution
End of explanation
"""


def wallis_rational(N):
    """
    The partial product approximation to pi using the first N terms of Wallis' formula.
    
    Parameters
    ----------
    
    N : int
        Number of terms in product
        
    Returns
    -------
    
    partial : Rational
        A rational number approximation to pi
    """
    
    partial = Rational(2,1)
    for n in range(1, N+1):
        partial = partial * Rational((2*n)**2, (2*n-1)*(2*n+1))
    return partial

pi_list = [wallis_rational(n) for n in range(1, 21)]
print(pi_list)
print(sorted(pi_list))

import numpy
print(numpy.pi-numpy.array(list(map(float, pi_list))))

"""
Explanation: Exercise 8
The Wallis formula for $\pi$ is
\begin{equation}
  \pi = 2 \prod_{n=1}^{\infty} \frac{ (2 n)^2 }{(2 n - 1) (2 n + 1)}.
\end{equation}
We can define a partial product $\pi_N$ as
\begin{equation}
  \pi_N = 2 \prod_{n=1}^{N} \frac{ (2 n)^2 }{(2 n - 1) (2 n + 1)},
\end{equation}
each of which are rational numbers.
Construct a list of the first 20 rational number approximations to $\pi$ and print them out. Print the sorted list to show that the approximations are always increasing. Then convert them to floating point numbers, construct a numpy array, and subtract this array from $\pi$ to see how accurate they are.
Solution
End of explanation
"""


lhs = 27**5 + 84**5 + 110**5 + 133**5
rhs = 144**5

print("Does the LHS {} equal the RHS {}? {}".format(lhs, rhs, lhs==rhs))

"""
Explanation: The shortest published Mathematical paper
A candidate for the shortest mathematical paper ever shows the following result:
\begin{equation}
  27^5 + 84^5 + 110^5 + 133^5 = 144^5.
\end{equation}
This is interesting as

This is a counterexample to a conjecture by Euler ... that at least $n$ $n$th powers are required to sum to an $n$th power, $n > 2$.

Exercise 1
Using python, check the equation above is true.
Solution
End of explanation
"""


import numpy
import itertools

"""
Explanation: Exercise 2
The more interesting statement in the paper is that
\begin{equation}
  27^5 + 84^5 + 110^5 + 133^5 = 144^5.
\end{equation}

[is] the smallest instance in which four fifth powers sum to a fifth power.

Interpreting "the smallest instance" to mean the solution where the right hand side term (the largest integer) is the smallest, we want to use python to check this statement.
You may find the combinations function from the itertools package useful.
End of explanation
"""


input_list = numpy.arange(1, 7)
combinations = list(itertools.combinations(input_list, 4))
print(combinations)

"""
Explanation: The combinations function returns all the combinations (ignoring order) of r elements from a given list. For example, take a list of length 6, [1, 2, 3, 4, 5, 6] and compute all the combinations of length 4:
End of explanation
"""


n_combinations = 144*143*142*141/24
print("Number of combinations of 4 objects from 144 is {}".format(n_combinations))

"""
Explanation: We can already see that the number of terms to consider is large.
Note that we have used the list function to explicitly get a list of the combinations. The combinations function returns a generator, which can be used in a loop as if it were a list, without storing all elements of the list.
How fast does the number of combinations grow? The standard formula says that for a list of length $n$ there are
\begin{equation}
  \begin{pmatrix} n \ k \end{pmatrix} = \frac{n!}{k! (n-k)!}
\end{equation}
combinations of length $k$. For $k=4$ as needed here we will have $n (n-1) (n-2) (n-3) / 24$ combinations. For $n=144$ we therefore have
End of explanation
"""


from matplotlib import pyplot
%matplotlib inline

n = numpy.arange(5, 51)
N = numpy.zeros_like(n)
for i, n_c in enumerate(n):
    combinations = list(itertools.combinations(numpy.arange(1,n_c+1), 4))
    N[i] = len(combinations)

pyplot.figure(figsize=(12,6))
pyplot.loglog(n, N, linestyle='None', marker='x', color='k', label='Combinations')
pyplot.loglog(n, n**4, color='b', label=r'$n^4$')
pyplot.xlabel(r'$n$')
pyplot.ylabel(r'$N$')
pyplot.legend(loc='upper left')
pyplot.show()

"""
Explanation: Exercise 2a
Show, by getting python to compute the number of combinations $N = \begin{pmatrix} n \ 4 \end{pmatrix}$ that $N$ grows roughly as $n^4$. To do this, plot the number of combinations and $n^4$ on a log-log scale. Restrict to $n \le 50$.
Solution
End of explanation
"""


nmax=145
range_to_power = numpy.arange(1, nmax)**5
lhs_combinations = list(itertools.combinations(range_to_power, 4))

"""
Explanation: With 17 million combinations to work with, we'll need to be a little careful how we compute.
One thing we could try is to loop through each possible "smallest instance" (the term on the right hand side) in increasing order. We then check all possible combinations of left hand sides.
This is computationally very expensive as we repeat a lot of calculations. We repeatedly recalculate combinations (a bad idea). We repeatedly recalculate the powers of the same number.
Instead, let us try creating the list of all combinations of powers once.
Exercise 2b

Construct a numpy array containing all integers in $1, \dots, 144$ to the fifth power. 
Construct a list of all combinations of four elements from this array.
Construct a list of sums of all these combinations.
Loop over one list and check if the entry appears in the other list (ie, use the in keyword).

Solution
End of explanation
"""


lhs_sums = []
for lhs_terms in lhs_combinations:
    lhs_sums.append(numpy.sum(numpy.array(lhs_terms)))

"""
Explanation: Then calculate the sums:
End of explanation
"""


for i, lhs in enumerate(lhs_sums):
    if lhs in range_to_power:
        rhs_primitive = int(lhs**(0.2))
        lhs_primitive = (numpy.array(lhs_combinations[i])**(0.2)).astype(int)
        print("The LHS terms are {}.".format(lhs_primitive))
        print("The RHS term is {}.".format(rhs_primitive))

"""
Explanation: Finally, loop through the sums and check to see if it matches any possible term on the RHS:
End of explanation
"""


def dvdt(v, t, sigma, rho, beta):
    """
    Define the Lorenz system.
    
    Parameters
    ----------
    
    v : list
        State vector
    t : float
        Time
    sigma : float
        Parameter
    rho : float
        Parameter
    beta : float
        Parameter
    
    Returns
    -------
    
    dvdt : list
        RHS defining the Lorenz system
    """
    
    x, y, z = v
    
    return [sigma*(y-x), x*(rho-z)-y, x*y-beta*z]

"""
Explanation: Lorenz attractor
The Lorenz system is a set of ordinary differential equations which can be written
\begin{equation}
  \frac{\text{d} \vec{v}}{\text{d} \vec{t}} = \vec{f}(\vec{v})
\end{equation}
where the variables in the state vector $\vec{v}$ are
\begin{equation}
  \vec{v} = \begin{pmatrix} x(t) \ y(t) \ z(t) \end{pmatrix}
\end{equation}
and the function defining the ODE is
\begin{equation}
  \vec{f} = \begin{pmatrix} \sigma \left( y(t) - x(t) \right) \ x(t) \left( \rho - z(t) \right) - y(t) \ x(t) y(t) - \beta z(t) \end{pmatrix}.
\end{equation}
The parameters $\sigma, \rho, \beta$ are all real numbers.
Exercise 1
Write a function dvdt(v, t, params) that returns $\vec{f}$ given $\vec{v}, t$ and the parameters $\sigma, \rho, \beta$.
Solution
End of explanation
"""


import numpy
from scipy.integrate import odeint

v0 = [1.0, 1.0, 1.0]
sigma = 10.0
beta = 8.0/3.0
t_values = numpy.linspace(0.0, 100.0, 5000)
rho_values = [13.0, 14.0, 15.0, 28.0]
v_values = []
for rho in rho_values:
    params = (sigma, rho, beta)
    v = odeint(dvdt, v0, t_values, args=params)
    v_values.append(v)

%matplotlib inline
from matplotlib import pyplot
from mpl_toolkits.mplot3d.axes3d import Axes3D

fig = pyplot.figure(figsize=(12,6))
for i, v in enumerate(v_values):
    ax = fig.add_subplot(2,2,i+1,projection='3d')
    ax.plot(v[:,0], v[:,1], v[:,2])
    ax.set_xlabel(r'$x$')
    ax.set_ylabel(r'$y$')
    ax.set_zlabel(r'$z$')
    ax.set_title(r"$\rho={}$".format(rho_values[i]))
pyplot.show()

"""
Explanation: Exercise 2
Fix $\sigma=10, \beta=8/3$. Set initial data to be $\vec{v}(0) = \vec{1}$. Using scipy, specifically the odeint function of scipy.integrate, solve the Lorenz system up to $t=100$ for $\rho=13, 14, 15$ and $28$.
Plot your results in 3d, plotting $x, y, z$.
Solution
End of explanation
"""


t_values = numpy.linspace(0.0, 40.0, 4000)
rho = 28.0
params = (sigma, rho, beta)
v_values = []
v0_values = [[1.0,1.0,1.0],
             [1.0+1e-5,1.0+1e-5,1.0+1e-5]]
for v0 in v0_values:
    v = odeint(dvdt, v0, t_values, args=params)
    v_values.append(v)

fig = pyplot.figure(figsize=(12,6))
line_colours = 'by'
for tstart in range(4):
    ax = fig.add_subplot(2,2,tstart+1,projection='3d')
    for i, v in enumerate(v_values):
        ax.plot(v[tstart*1000:(tstart+1)*1000,0], 
                v[tstart*1000:(tstart+1)*1000,1], 
                v[tstart*1000:(tstart+1)*1000,2], 
                color=line_colours[i])
    ax.set_xlabel(r'$x$')
    ax.set_ylabel(r'$y$')
    ax.set_zlabel(r'$z$')
    ax.set_title(r"$t \in [{},{}]$".format(tstart*10, (tstart+1)*10))
pyplot.show()

"""
Explanation: Exercise 3
Fix $\rho = 28$. Solve the Lorenz system twice, up to $t=40$, using the two different initial conditions $\vec{v}(0) = \vec{1}$ and $\vec{v}(0) = \vec{1} + \vec{10^{-5}}$.
Show four plots. Each plot should show the two solutions on the same axes, plotting $x, y$ and $z$. Each plot should show $10$ units of time, ie the first shows $t \in [0, 10]$, the second shows $t \in [10, 20]$, and so on.
Solution
End of explanation
"""


import sympy
sympy.init_printing()

y, t = sympy.symbols('y, t')

sympy.dsolve(sympy.diff(y(t), t) + y(t)**2 - sympy.exp(-t), y(t))

"""
Explanation: This shows the sensitive dependence on initial conditions that is characteristic of chaotic behaviour.
Systematic ODE solving with sympy
We are interested in the solution of
\begin{equation}
  \frac{\text{d} y}{\text{d} t} = e^{-t} - y^n, \qquad y(0) = 1,
\end{equation}
where $n > 1$ is an integer. The "minor" change from the above examples mean that sympy can only give the solution as a power series.
Exercise 1
Compute the general solution as a power series for $n = 2$.
Solution
End of explanation
"""


for n in range(2, 11):
    ode_solution = sympy.dsolve(sympy.diff(y(t), t) + y(t)**n - sympy.exp(-t), y(t), 
                                ics = {y(0) : 1})
    print(ode_solution)

"""
Explanation: Exercise 2
Investigate the help for the dsolve function to straightforwardly impose the initial condition $y(0) = 1$ using the ics argument. Using this, compute the specific solutions that satisfy the ODE for $n = 2, \dots, 10$.
Solution
End of explanation
"""


%matplotlib inline

for n in range(2, 11):
    ode_solution = sympy.dsolve(sympy.diff(y(t), t) + y(t)**n - sympy.exp(-t), y(t), 
                                ics = {y(0) : 1})
    sympy.plot(ode_solution.rhs.removeO(), (t, 0, 1));

"""
Explanation: Exercise 3
Using the removeO command, plot each of these solutions for $t \in [0, 1]$.
End of explanation
"""


def all_primes(N):
    """
    Return all primes less than or equal to N.
    
    Parameters
    ----------
    
    N : int
        Maximum number
        
    Returns
    -------
    
    prime : generator
        Prime numbers
    """
    
    primes = []
    for n in range(2, N+1):
        is_n_prime = True
        for p in primes:
            if n%p == 0:
                is_n_prime = False
                break
        if is_n_prime:
            primes.append(n)
            yield n

"""
Explanation: Twin primes
A twin prime is a pair $(p_1, p_2)$ such that both $p_1$ and $p_2$ are prime and $p_2 = p_1 + 2$.
Exercise 1
Write a generator that returns twin primes. You can use the generators above, and may want to look at the itertools module together with its recipes, particularly the pairwise recipe.
Solution
Note: we need to first pull in the generators introduced in that notebook
End of explanation
"""


from itertools import tee

def pair_primes(N):
    "Generate consecutive prime pairs, using the itertools recipe"
    a, b = tee(all_primes(N))
    next(b, None)
    return zip(a, b)

"""
Explanation: Now we can generate pairs using the pairwise recipe:
End of explanation
"""


def check_twin(pair):
    """
    Take in a pair of integers, check if they differ by 2.
    """
    p1, p2 = pair
    return p2-p1 == 2

"""
Explanation: We could examine the results of the two primes directly. But an efficient solution is to use python's filter function. To do this, first define a function checking if the pair are twin primes:
End of explanation
"""


def twin_primes(N):
    """
    Return all twin primes
    """
    return filter(check_twin, pair_primes(N))

"""
Explanation: Then use the filter function to define another generator:
End of explanation
"""


for tp in twin_primes(20):
    print(tp)

"""
Explanation: Now check by finding the twin primes with $N<20$:
End of explanation
"""


def pi_N(N):
    """
    Use the quantify pattern from itertools to count the number of twin primes.
    """
    return sum(map(check_twin, pair_primes(N)))

pi_N(1000)

"""
Explanation: Exercise 2
Find how many twin primes there are with $p_2 < 1000$.
Solution
Again there are many solutions, but the itertools recipes has the quantify pattern. Looking ahead to exercise 3 we'll define:
End of explanation
"""


import numpy
from matplotlib import pyplot
%matplotlib inline

N = numpy.array([2**k for k in range(4, 17)])
twin_prime_fraction = numpy.array(list(map(pi_N, N))) / N

pyplot.semilogx(N, twin_prime_fraction)
pyplot.xlabel(r"$N$")
pyplot.ylabel(r"$\pi_N / N$")
pyplot.show()

"""
Explanation: Exercise 3
Let $\pi_N$ be the number of twin primes such that $p_2 < N$. Plot how $\pi_N / N$ varies with $N$ for $N=2^k$ and $k = 4, 5, \dots 16$. (You should use a logarithmic scale where appropriate!)
Solution
We've now done all the hard work and can use the solutions above.
End of explanation
"""


pyplot.semilogx(N, twin_prime_fraction * numpy.log(N)**2)
pyplot.xlabel(r"$N$")
pyplot.ylabel(r"$\pi_N \times \log(N)^2 / N$")
pyplot.show()

"""
Explanation: For those that have checked Wikipedia, you'll see Brun's theorem which suggests a specific scaling, that $\pi_N$ is bounded by $C N / \log(N)^2$. Checking this numerically on this data:
End of explanation
"""


class Polynomial(object):
    """Representing a polynomial."""
    explanation = "I am a polynomial"
    
    def __init__(self, roots, leading_term):
        self.roots = roots
        self.leading_term = leading_term
        self.order = len(roots)
        
    def __repr__(self):
        string = str(self.leading_term)
        for root in self.roots:
            if root == 0:
                string = string + "x"
            elif root > 0:
                string = string + "(x - {})".format(root)
            else:
                string = string + "(x + {})".format(-root)
        return string
    
    def __mul__(self, other):
        roots = self.roots + other.roots
        leading_term = self.leading_term * other.leading_term
        return Polynomial(roots, leading_term)
    
    def explain_to(self, caller):
        print("Hello, {}. {}.".format(caller,self.explanation))
        print("My roots are {}.".format(self.roots))
        return None

class Monomial(Polynomial):
    """Representing a monomial, which is a polynomial with leading term 1."""
    explanation = "I am a monomial"
    
    def __init__(self, roots):
        Polynomial.__init__(self, roots, 1)
        
    def __repr__(self):
        string = ""
        for root in self.roots:
            if root == 0:
                string = string + "x"
            elif root > 0:
                string = string + "(x - {})".format(root)
            else:
                string = string + "(x + {})".format(-root)
        return string

"""
Explanation: A basis for the polynomials
In the section on classes we defined a Monomial class to represent a polynomial with leading coefficient $1$. As the $N+1$ monomials $1, x, x^2, \dots, x^N$ form a basis for the vector space of polynomials of order $N$, $\mathbb{P}^N$, we can use the Monomial class to return this basis.
Exercise 1
Define a generator that will iterate through this basis of $\mathbb{P}^N$ and test it on $\mathbb{P}^3$.
Solution
Again we first take the definition of the crucial class from the notes.
End of explanation
"""


def basis_pN(N):
    """
    A generator for the simplest basis of P^N.
    """
    
    for n in range(N+1):
        yield Monomial(n*[0])

"""
Explanation: Now we can define the first basis:
End of explanation
"""


for poly in basis_pN(3):
    print(poly)

"""
Explanation: Then test it on $\mathbb{P}^N$:
End of explanation
"""


class Monomial(Polynomial):
    """Representing a monomial, which is a polynomial with leading term 1."""
    explanation = "I am a monomial"
    
    def __init__(self, roots):
        Polynomial.__init__(self, roots, 1)
        
    def __repr__(self):
        if len(self.roots):
            string = ""
            n_zero_roots = len(self.roots) - numpy.count_nonzero(self.roots)
            if n_zero_roots == 1:
                string = "x"
            elif n_zero_roots > 1:
                string = "x^{}".format(n_zero_roots)
        else: # Monomial degree 0.
            string = "1"
        for root in self.roots:
            if root > 0:
                string = string + "(x - {})".format(root)
            elif root < 0:
                string = string + "(x + {})".format(-root)
        return string

"""
Explanation: This looks horrible, but is correct. To really make this look good, we need to improve the output. If we use
End of explanation
"""


for poly in basis_pN(3):
    print(poly)

"""
Explanation: then we can deal with the uglier cases, and re-running the test we get
End of explanation
"""


def basis_pN_variant(N):
    """
    A generator for the 'sum' basis of P^N.
    """
    
    for n in range(N+1):
        yield Monomial(range(n+1))

for poly in basis_pN_variant(4):
    print(poly)

"""
Explanation: An even better solution would be to use the numpy.unique function as in this stackoverflow answer (the second one!) to get the frequency of all the roots.
Exercise 2
An alternative basis is given by the monomials
\begin{align}
  p_0(x) &= 1, \ p_1(x) &= 1-x, \ p_2(x) &= (1-x)(2-x), \ \dots & \quad \dots, \ p_N(x) &= \prod_{n=1}^N (n-x).
\end{align}
Define a generator that will iterate through this basis of $\mathbb{P}^N$ and test it on $\mathbb{P}^4$.
Solution
End of explanation
"""


from itertools import product

def basis_product():
    """
    Basis of the product space
    """
    yield from product(basis_pN(3), basis_pN_variant(4))

for p1, p2 in basis_product():
    print("Basis element is ({}) X ({}).".format(p1, p2))

"""
Explanation: I am too lazy to work back through the definitions and flip all the signs; it should be clear how to do this!
Exercise 3
Use these generators to write another generator that produces a basis of $\mathbb{P^3} \times \mathbb{P^4}$.
Solution
Hopefully by now you'll be aware of how useful itertools is!
End of explanation
"""


def basis_product_long_form():
    """
    Basis of the product space (without using yield_from)
    """
    prod = product(basis_pN(3), basis_pN_variant(4))
    yield next(prod)

for p1, p2 in basis_product():
    print("Basis element is ({}) X ({}).".format(p1, p2))

"""
Explanation: I've cheated here as I haven't introduced the yield from syntax (which returns an iterator from a generator). We could write this out instead as
End of explanation
"""


