Assignment 0: Introduction

Assignment 0: Introduction#

The purpose of this assignment to assess your current program skills and to get you started programming as quickly as possible. Complete this assignment using your CoCalc project for the course so you can immediately start programming without having to install anything.

If you have a strong preference for a language, the please feel free to use this, but note that we can only provide full support for Python at the moment. (Ask your instructor if a kernel does not exist on CoCalc.)

Basic Programming Constructs#

Note

For the tasks in this section, please only use standard features of the language and libraries like math. I.e., do not use numpy, scipy etc. in this section.

Functions#

Write a function \(\mathrm{sq}(x)\) that returns \(x^2\). Be sure to document your function, with comments if necessary. Are there multiple ways to define functions in your language? Does your language support keyword arguments? Default values? Can you document your function? If your language is strongly typed, be sure to implement a version that works on integers, floating point numbers (double precision), and complex numbers.
Use your function to compute \(5^2\), \(\pi^2\), and \((1+\I)^2\).

Write another function apply(f, n, x) that applies the function f to x n times:

\[\begin{gather*} \mathrm{apply}(f, n, x) = \overbrace{f(f(\cdots f}^{n \text{ times}}(x)\cdots)). \end{gather*}\]

Can you solve this by writing a higher-order function compose(f, n) that returns another function which applies f n times?

\[\begin{gather*} \mathrm{compose}(f, n)(x) = \overbrace{f(f(\cdots f}^{n \text{ times}}(x)\cdots)). \end{gather*}\]

Use your solutions to compute \(5\), \(5^{4}\), \(5^{32}\), etc., applying your \(\textrm{sq}\) function for \(n\in {0, 2, 5}\) to the previous numbers.

Discussion.

Different languages have different support for functions. For example:

C and C++ require you to define separate functions for different types – you must define a different function for integers and floating point numbers – although there are often mechanisms like templates that allow you to not have to write duplicate code;
Most languages allow you to pass functions around as variables, often with function pointers, but some, like MATLAB, requires that you define functions in a separate file. This might make it hard to write compose.
Play a little bit with your language: does it support arbitrary precision integers? (Can you compute \(5^{100}\)?) Does it support different precision floating point numbers (single, double, quadruple, …)?

Iteration and Recursion#

The Fibonacci numbers are defined by the following recurrence relationship

\[\begin{gather*} F_0 = 0, \qquad F_1 = 1, \qquad F_n = F_{n-1} + F_{n-2}. \end{gather*}\]

I.e., each successive number is the sum of the previous two numbers.

Use recursion to compute \(F_{100} = 354224848179261915075\).
Use iteration to compute a list \([F_n | n \in \{0, 1, 2, \dots, 100\}] = [0, 1, 1, 2, 3, 5, \dots, 354224848179261915075]\).

Challenge

What is the largest Fibonacci number you can compute in 1ms on CoCalc? Note: this is a bit of a rabbit hole. The naïve implementation is slow, but more careful implementations will easily start challenging your ability accurately measure the performance of your method. the timeit library and IPython %timeit magic will help, but now you will need to be clever in your approach to finding the maximum n.

# Here is a simple code to approximate n
import time
for n in range(100):
    tic = time.time()
    fn = fib(n)
    t = time.time() - tic
    if t > 1e-3:
        n -= 1
        break
print(n, fn)

# Check that it was not a fluke:
%timeit fib(n)

20 10946
740 μs ± 3.4 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Here is an iterative solution

def fib1(n):
    fn = [1, 1]
    while len(fn) < n:
        fn.append(fn[-2] + fn[-1])
    return fn[-1]

assert fib1(20) == fib(20)

for n in range(100):
    tic = time.time()
    fn = fib1(n)
    t = time.time() - tic
    if t > 1e-3:
        n -= 1
        break
print(n, fn)

# Check that it was not a fluke:
%timeit fib1(n)

99 218922995834555169026
6.34 μs ± 20.1 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

Here is another approach (don’t peak until you give up). It gives the skeleton of a very fast algorithm.

fib1(2**14)
6.9 ms ± 132 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
fib2
179 μs ± 597 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

Helpers#

Here is some code for more accurately timing.

import timeit

def mytimeit(n, fib=fib, number=10, samples=5):
    """Return an estimate of the time it takes to compute `fib(n)`"""
    ts = [
        timeit.timeit('fib(n)', globals=dict(fib=fib, n=n), number=number)
        for _n in range(samples)
    ]
    return min(ts) / number

def get_n(fib, target=1e-3):  # 1ms
    """Return the largest `n` where `fib(n)` takes < target."""
    n0 = 1
    t0 = mytimeit(n0, fib=fib)
    n1 = 2
    t1 = mytimeit(n1, fib=fib)
    while t1 < target:
        n0, t0 = n1, t1
        n1 *= 2
        t1 = mytimeit(n1, fib=fib)
    
    # Now n0 and n1 should bracket our desired n.  We can use bisection.
    while n1 - n0 > 1:
        n = (n1 + n0)//2
        t = mytimeit(n, fib=fib)
        if t < target:
            n0 = n
        else:
            n1 = n
    return n0

print(get_n(fib))
print(get_n(fib1))

20
7366

Assignment Preparation#

Complete the following in one of two ways:

A simple brute force solution that you completely understand but which may not be fully accurate.
A more sophisticated solution that is accurate to close to machine precision.

For the second method, please feel free to use any tools in the NumPy or SciPy libraries. (You may use mpmath to check, but don’t use it as part of your solution.)

Floating Point#

Series#

Numerically check the following formula:

\[ \sum_{n=1}^{M} n = \frac{M(M+1)}{2} \]

\[ \sum_{n=1}^{\infty} \frac{1}{n^{2}} = \zeta(2) = \frac{\pi^2}{6} \]

where \(\zeta(s)\) is the Riemann zeta function.

Integrals#

Numerically check the following integrals for all viable values of \(p\) (both positive and negative):

\[ \int_0^1 x^p \d{x} = \frac{1}{p+1} \]

What about infinite limits?

\[ \int_{-\infty}^{\infty} e^{-x^2} \d{x} = \sqrt{\pi} \]

This one is probably challenging to do accurately:

\[ \int_{-\infty}^{\infty} e^{-x^2}\sin\frac{1}{x^2}\d{x} = \sqrt{\pi}e^{-\sqrt{2}}\sin(\sqrt{2}) \]

This should be easy if you could do the previous one:

\[\begin{split} \int_0^{\infty} xe^{-x^2}\sin^2\frac{1}{x}\d{x} = \frac{\sqrt{\pi}}{4} G_{14}^{13}\left( \left. \begin{matrix} \tfrac{1}{2}\\ \tfrac{1}{2} & \tfrac{1}{2} & 0 & -\tfrac{1}{2} \end{matrix} \right| z=1 \right) = 0.32006330909018418888\cdots, %37810082287661243051390948375855688739690521501945667936210716620524245639515708 \end{split}\]

where \(G_{mn}^{pq}\bigl(\begin{smallmatrix}\vect{a}_p\\ \vect{b}_q\end{smallmatrix}\big|z\bigr)\) is the Meijer G-function.

Roots#

Find all solutions \(x\) to the following equations:

Polynomials#

\[ x^2 - (1+\epsilon)x + \epsilon = 0, \qquad x = \{1, \epsilon\} \]

for \(\epsilon \in \{1, 10^{-10}, 10^{-20}\}\).

\[ x^4 + 3x^2 - 6x + 10 = 0, \qquad x = \{1\pm\I, -1\pm 2\I\} \]

Lambert W function#

\[ xe^x = w, \qquad x = L_0(1) \]

for \(w = 1\).

The function \(w = W_k(x)\) is the Lambert W function.

More Details#

Write a function phys_581.assignment_0.quadratic_equation() which returns the roots of the equation

\[ ax^2 + bx + c = 0. \]

For example:

\[ x^2 - 3x + 2 = (x-2)(x-1) \]

so we expect

from phys_581.assignment_0 import quadratic_equation
np.allclose(quadratic_equation(a=1, b=-3, c=2), [1, 2])

True

Note: if you attempt to blindly use the quadratic formula as is done in the default version, you will encounter errors when \(b \approx \pm \sqrt{b^2 - 4ac}\) because the two terms can cancel. This is the main source of error associated with floating point comptations.

When will this become a problem?
How can you overcome this issue?
How will you test your code to make sure it works well?

The goal should be for every function to return an answer that has a relative error comparable to the machine precision of the computer: sometimes called \(\epsilon = \)eps. This is not always possible if the problem is ill-conditioned.

for epsilon in [1, 1e-10, 1e-20]:
    roots = quadratic_equation(a=1, b=-(1+epsilon), c=epsilon)
    exact_roots = [1, epsilon]
    print(f"epsilon = {epsilon}")
    print(f"roots =       {roots}")
    print(f"exact roots = {exact_roots}")
    print(f"Roots allclose? {np.allclose(roots, exact_roots)}")
    print()

epsilon = 1
roots =       (array(1.+0.j), array(1.+0.j))
exact roots = [1, 1]
Roots allclose? True

epsilon = 1e-10
roots =       (array(1.00000008e-10+0.j), array(1.+0.j))
exact roots = [1, 1e-10]
Roots allclose? False

epsilon = 1e-20
roots =       (array(0.+0.j), array(1.+0.j))
exact roots = [1, 1e-20]
Roots allclose? False

Floating Point Numbers#

Readings:

[Gezerlis, 2020] Chapter 2. Try some of the “experiments” Alex suggests.
[Goldberg, 1991] “What Every Computer Scientist Should Know About Floating-Point Numbers.” This is a rather technical, but complete account about floating point numbers, and contains a detailed analysis of this assignment.

We can use NumPy to see the properties of the floating point numbers. The defult float in python is the IEEE double-precision floating point number which has \(64 = 52 + 11 + 1\) bits. 52 of these are used for the mantissa, 11 for the exponent and 1 for the sign. Roughly, the relative error in a floating point number is eps=\(\epsilon = 1/2^{52}\approx 2.22\times 10^{-16}\) (16 digits of precision in decimal), with an exponent that can range from \(\pm 2^{10} = \pm 1024\) with a smallest value of tiny\(\approx 2^{-1024} \approx 5.56\times 10^{-309}\) and largest value of max\(=2^{1024}\approx 1.798\times 10^{208}\).

Try to compute these numbers! In particular, 2**1024.0 causes an overflow… how can you compute \(2^{1024}\approx 1.798\times 10^{208}\) using floating point numbers? Hint:

\[ 2^n = a\times 10^m, \qquad n\log_{10}(2) = \log_{10}(a) + m. \]

The actual values are a little different as we see below (e.g. tiny\(=2^{-1022}=2.225\times 10^{-308}\)) because of some subtleties in the actual implementation of the IEEE standard.

import numpy as np
print(np.finfo(float))

Machine parameters for float64
---------------------------------------------------------------
precision = 15   resolution = 1e-15
machep = -52   eps =        2.220446049250313e-16
negep =  -53   epsneg =     1.1102230246251565e-16
minexp = -1022   tiny =       2.2250738585072014e-308
maxexp = 1024   max =        1.7976931348623157e+308
nexp =   11   min =        -max
smallest_normal = 2.2250738585072014e-308   smallest_subnormal = 5e-324
---------------------------------------------------------------