--- jupytext: formats: ipynb,md:myst text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.13.6 kernelspec: display_name: Python 3 (phys-581) language: python name: phys-581 --- ```{code-cell} :tags: [hide-cell] # This cell contains some initializations for the notebook so that # the code can be executed for the documentation. import sys; print(sys.executable) # This is useful for debugging the installation #import mmf_setup;mmf_setup.nbinit(quiet=True) # This loads some useful tools for me. # Sometimes matplotlib logging gets tedius: we limit this to critical errors. import logging;logging.getLogger('matplotlib').setLevel(logging.CRITICAL) # Now we initialize plotting with the inline (image) backend %matplotlib inline # To actually do plotting, we import numpy and matplotlib. These shorthands -- # `np` and `plt` -- are quite standard. import numpy as np, matplotlib.pyplot as plt ``` :::{margin} For more details, see {ref}`sec:CoCalc`. ::: (sec:Assignment0)= # 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 {py:mod}`math`. I.e., do not use {py:mod}`numpy`, {py:mod}`scipy` etc. in this section. ::: ### Functions 1. 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. 2. Use your function to compute $5^2$, $\pi^2$, and $(1+\I)^2$. ```{code-cell} :tags: [hide-input] import math def sq(x): """Return the square of x.""" return x**2 assert sq(5) == 25 assert sq(math.pi) == math.pi**2 assert sq(1.0+1.0j) == 2.0j # Anonymous function that returns the square of x. f = lambda x: x**2 assert sq(5) == 25 assert sq(math.pi) == math.pi**2 assert sq(1.0+1.0j) == 2.0j ``` :::{margin} Depending on your language, writing such a function might be quite challenging. Please let you instructor know if your language makes this difficult. ::: 3. 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*} :::{margin} Note that we use the reasonable convention $\mathrm{apply}(f, x, 0) = x$ here. ::: 4. 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. ```{code-cell} :tags: [hide-cell] from IPython.display import display, Latex def apply(f, n, x): """Return f applied to x n times.""" res = x for i in range(n): res = f(res) return res def compose(f, n): """Return a function that applies f, n times.""" def fn(x): res = x for i in range(n): res = f(res) return res return fn for x in [5, math.pi, 1+1j]: res = [] for n in [0, 2, 5]: assert compose(sq, n)(x) == apply(sq, n, x) res.append(f"f^{{{n}}}({{x}}) = {apply(sq, n, x)}") display(Latex("${}$".format(r",\quad ".join(res)))) ``` :::{admonition} 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. 1. Use recursion to compute $F_{100} = 354224848179261915075$. 2. Use iteration to compute a list $[F_n | n \in \{0, 1, 2, \dots, 100\}] = [0, 1, 1, 2, 3, 5, \dots, 354224848179261915075]$. ```{code-cell} :tags: [hide-cell] def fib(n): if n <= 1: return n return fib(n-2) + fib(n-1) print(fib(30)) print([fib(n) for n in range(31)]) ``` :::{admonition} 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 {py:mod}`timeit` library and IPython `%timeit` magic will help, but now you will need to be clever in your approach to finding the maximum `n`. ::: ```{code-cell} # 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) ``` Here is an iterative solution ```{code-cell} 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) ``` Here is another approach (don't peak until you give up). It gives the skeleton of a very fast algorithm. ```{code-cell} :tags: [hide-input] def fib2(log2_n): """Return fib(2**n).""" f00, f01, f10, f11 = [0,1,1,1] for i in range(log2_n): f00, f01, f10, f11 = (f00*f00 + f01*f10, f00*f01 + f01*f11, f10*f00 + f11*f10, f10*f01 + f11*f11) return f01 assert fib1(2**14) == fib2(14) print("fib1(2**14)") %timeit fib1(2**14) print("fib2") %timeit fib2(14) ``` ### Helpers Here is some code for more accurately timing. ```{code-cell} 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)) ``` ## Assignment Preparation Complete the following in one of two ways: 1. A simple brute force solution that you completely understand but which may not be fully accurate. 2. 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](https://numpy.org/doc/stable/) or [SciPy](https://docs.scipy.org/doc/scipy/reference/) libraries. (You may use [mpmath](https://mpmath.org) 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} $$ ```{code-cell} ``` $$ \sum_{n=1}^{\infty} \frac{1}{n^{2}} = \zeta(2) = \frac{\pi^2}{6} $$ where $\zeta(s)$ is the [Riemann zeta function](https://en.wikipedia.org/wiki/Riemann_zeta_function). ```{code-cell} ``` ### 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} $$ ```{code-cell} ``` What about infinite limits? $$ \int_{-\infty}^{\infty} e^{-x^2} \d{x} = \sqrt{\pi} $$ ```{code-cell} ``` 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}) $$ ```{code-cell} ``` This should be easy if you could do the previous one: $$ \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 $$ where $G_{mn}^{pq}\bigl(\begin{smallmatrix}\vect{a}_p\\ \vect{b}_q\end{smallmatrix}\big|z\bigr)$ is the [Meijer G-function](https://en.wikipedia.org/wiki/Meijer_G-function). ```{code-cell} ``` ### 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}\}$. ```{code-cell} ``` $$ x^4 + 3x^2 - 6x + 10 = 0, \qquad x = \{1\pm\I, -1\pm 2\I\} $$ ```{code-cell} ``` #### 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](https://en.wikipedia.org/wiki/Lambert_W_function). +++ ## More Details +++ Write a function {func}`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 ```{code-cell} from phys_581.assignment_0 import quadratic_equation np.allclose(quadratic_equation(a=1, b=-3, c=2), [1, 2]) ``` 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. 1. When will this become a problem? 2. How can you overcome this issue? 3. 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. ```{code-cell} 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() ``` ### Floating Point Numbers +++ **Readings:** * {cite:p}`Gezerlis:2020` Chapter 2. Try some of the "experiments" Alex suggests. * {cite:p}`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. ```{code-cell} import numpy as np print(np.finfo(float)) ``` ```{code-cell} ``` [MATLAB]: [CoCalc]: