02_NumberRepresentation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Integers\n",
    "\n",
    "Integer numbers are represented by N bit words. Python3 allows you to store integers with practically **unlimited precision**, the only limitation comes from the (contiguous) space available in memory.\n",
    "In Python2 (deprecated), N depends on the PC architercture, N=64 in modern computers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9223372036854775807\n",
      "Is your system a 64 bit one? True\n",
      "9223372036854775808\n"
     ]
    }
   ],
   "source": [
    "# Check the largest integer\n",
    "import sys\n",
    "print(sys.maxsize)\n",
    "\n",
    "# Check also that corresponds to a 64-bit integer\n",
    "print(\"Is your system a 64 bit one?\", 2**63 - 1 == sys.maxsize)\n",
    "\n",
    "# Python 3 doesn't have a limit for integers\n",
    "maxint = sys.maxsize+1\n",
    "print(maxint)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Binary and Hexadecimal representations\n",
    "\n",
    "The common assumption is that numbers (in Python as in all the other languages) are expressed as decimal numbers. Built-in functions allows explicitly to convert from one base to another.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the binary representation, typically 1 bit ($j$) is dedicated to specifying the sign of the number, and the conversion between binary and decimal representation is:\n",
    "\n",
    "$$d = (-1)^j\\sum_{i=0}^{N-1} \\alpha_i ~ 2^i$$\n",
    "\n",
    "where $\\alpha_i$ are either 0 or 1. \n",
    "$b=\\alpha_{N-1}\\alpha_{N-2}..\\alpha_0$ is the binary representation of the number.\n",
    "\n",
    "Example: an 8-bit integer in binary representation with one bit for the sign:\n",
    "\n",
    "|  j | 6 | 5 | 4 | 3 | 2 | 1 | 0  |\n",
    "|---|---|---|---|---|---|---|---|\n",
    "|  0 | 0 | 0 | 1 | 0 | 1 | 1 | 1  |\n",
    "\n",
    "corresponds to:\n",
    "\n",
    "$$d = (-1)^j\\sum_{i=0}^{N-1} \\alpha_i ~ 2^i = (-1)^{0} [ (1) \\cdot 2^0 + (1) \\cdot 2^1 + (1) \\cdot 2^2 + (0) \\cdot 2^3 + (1) \\cdot 2^4 + (0) \\cdot 2^5 + (0) \\cdot 2^6] = 0 + 1 + 2 + 4 + 16 = 23$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binary representation of 12 : 0b1100\n",
      "Hexadecimal representation of 12 : 0xc\n",
      "Decimal representation of 0b1100 : 12\n",
      "Decimal representation of 0xc : 12\n",
      "Decimal representation of 88 : 88\n"
     ]
    }
   ],
   "source": [
    "# an integer in decimal representation\n",
    "a = 12\n",
    "\n",
    "# its binary representation\n",
    "a_bin = bin(a)\n",
    "print('Binary representation of', a, ':', a_bin)\n",
    "\n",
    "# its hexadecimal representation\n",
    "a_hex = hex(a)\n",
    "print('Hexadecimal representation of', a, ':', a_hex)\n",
    "\n",
    "# converting back to integer\n",
    "print('Decimal representation of', a_bin, ':', int(a_bin, 2))\n",
    "print('Decimal representation of', a_hex, ':', int(a_hex, 16))\n",
    "print('Decimal representation of',  \"88\", ':', int( '88', 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bitwise operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logical operators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bitwise AND  12\n",
      "Bitwise OR   61\n",
      "Bitwise XOR  49\n"
     ]
    }
   ],
   "source": [
    "a = 60           # 60 = 0011 1100 \n",
    "b = 13           # 13 = 0000 1101 \n",
    "                 #    = 0000 1100\n",
    "c = a & b        # 12 = 0000 1100\n",
    "print(\"Bitwise AND \", c)\n",
    "\n",
    "c = a | b        # 61 = 0011 1101 \n",
    "print(\"Bitwise OR  \", c)\n",
    "\n",
    "c = a ^ b        # 49 = 0011 0001\n",
    "print(\"Bitwise XOR \", c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unary operators\n",
    "\n",
    "#### bitwise NOT\n",
    "given a integer a:\n",
    "\n",
    "`\n",
    "~a = ~bin(a)\n",
    "   = -(bin(a)+1)\n",
    "`\n",
    "\n",
    "i.e. it returns the complement to (-) 1 of that number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bitwise NOT  -13 False\n"
     ]
    }
   ],
   "source": [
    "c = ~a           # -61 = 1100 0011\n",
    "d = not a\n",
    "print(\"Bitwise NOT \", c, d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Shifts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Left shift (towards most significant) of two positions  240\n",
      "Right shift (towards least significant) of two positions  15\n",
      "\n"
     ]
    }
   ],
   "source": [
    "c = a << 2       # 240 = 1111 0000\n",
    "print(\"Left shift (towards most significant) of two positions \", c)\n",
    "                 #    = 0011 1100\n",
    "c = a >> 2       # 15 = 0000 1111\n",
    "print(\"Right shift (towards least significant) of two positions \", c)\n",
    "\n",
    "print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more details check the python [documentation](https://realpython.com/python-bitwise-operators/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Floating point numbers\n",
    "\n",
    "Non-integer number **cannot be represented with infinite precision** on a computer. Single precision (also known as *float*) and double precision numbers use 32 and 64 bits respectively. \n",
    "Note that all floating point numbers in python are double precision (64 bits).\n",
    "A standard has been developed by IEEE such that the relative precision (see later) is the same across the whole validity range.\n",
    "\n",
    "The 32 or 64 bits are divided among 3 quantities uniquely characterizing the number:\n",
    "\n",
    "$x_{float} = (-1)^s \\times 1.f \\times 2^{e-bias}$\n",
    "\n",
    "where *s* is the sign, *f* the fractional part of the mantissa and *e* the exponent. In order to get numbers smaller than 1, a constant *bias* term is added to the exponent, such *bias* is typically equal to half of the max value of *e*.\n",
    "The mantissa is defined as:\n",
    "\n",
    "${\\rm mantissa}=1.f=1+m_{n-1}2^{-1}+m_{n-2}2^{-2}+..+m_{0}2^{-n}$\n",
    "\n",
    "where $n$ is the number of bits dedicated to *f* (see below) and $m_i$ are the binary coefficients. \n",
    "\n",
    "Numbers exceeding the maximum allowed value are *overflows* and the calculations involving them provide incorrect answers. Numbers smaller in absolute value than the minimum allowed value are *underflows* and simply set to zero, also in this case incorrect results are yielded.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Single\n",
    "\n",
    "For single precision floating point numbers, $0\\le e \\le 255$ and $bias=127$. Bits are arranged as follows:\n",
    "\n",
    "|   | *s* | *e* | *f* |\n",
    "|---|---|---|---|\n",
    "| Bit position | 31 | 30-23 | 22-0 |\n",
    "\n",
    "An example is given below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"http://www.dspguide.com/graphics/F_4_2.gif\"/>"
      ],
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "Image(url='http://www.dspguide.com/graphics/F_4_2.gif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Special values are also possibiles. N.B.: those are not numbers that can be used in the mathematical sense!\n",
    "\n",
    "|   |  conditions | value |\n",
    "|---|---|---|\n",
    "|  $+\\infty$ | s=0, e=255, f=0 | +INF  |\n",
    "|  $-\\infty$ | s=1, e=255, f=0 | -INF  |\n",
    "|  not a number | e=255, f>0  | NaN  |\n",
    "\n",
    "The largest value is obtained for $f\\sim 2$ and $e=254$, i.e. $2\\times2^{127}\\sim 3.4\\times10^{38}$.\n",
    "\n",
    "The value closest to zero is obtained instead for $f=2^{-23}$ and $e=0$, i.e. $2^{-149}\\sim 1.4\\times10^{-45}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Double\n",
    "\n",
    "For double precision floating point numbers, $0\\le e \\le 2047$ and $bias=1023$. Bits are arranged as follows:\n",
    "\n",
    "|   | *s* | *e* | *f* |\n",
    "|---|---|---|---|\n",
    "| Bit position | 63 | 62-52 | 51-0 |\n",
    "\n",
    "Special values are also possibiles. N.B.: those are not numbers that can be used in the mathematical sense!\n",
    "\n",
    "|   |  conditions | value |\n",
    "|---|---|---|\n",
    "|  $+\\infty$ | s=0, e=2047, f=0 | +INF  |\n",
    "|  $-\\infty$ | s=1, e=2047, f=0 | -INF  |\n",
    "|  not a number | e=2047, f>0  | NaN  |\n",
    "\n",
    "The validity range for double numbers is $2.2^{-308} - 1.8^{308}$\n",
    "\n",
    "Serious scientific calculations almost always requires at least double precision floating point numbers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Floating point numbers on your system\n",
    "\n",
    "Information about the floating point reresentation on your system can be obtained from sys.float_info. Definitions of the stored values are given on the python doc [page](https://docs.python.org/2/library/sys.html#sys.float_info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sys.float_info(max=1.7976931348623157e+308, max_exp=1024, max_10_exp=308, min=2.2250738585072014e-308, min_exp=-1021, min_10_exp=-307, dig=15, mant_dig=53, epsilon=2.220446049250313e-16, radix=2, rounds=1)\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "print(sys.float_info)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Accuracy and the perils of calculations with floats\n",
    "\n",
    "\n",
    "Floats can only have a limited number of meaningful decimal places, on the basis of how many bits are allocated for the fractional part of the mantissa: 6-7 decimal places for singles, 15-16 for doubles. In particular this means that calculations involving numbers with more than those decimal places involved do not yield the correct result, simply because the binary representation of those numbers does not allow to store them with sufficient accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.00000000000001\n",
      "7.000000000000001\n",
      "7.0\n"
     ]
    }
   ],
   "source": [
    "# Addding an increasingly small number to 7\n",
    "for e in [14, 15, 16]: print (7+1.0*10**-e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It should never been forgotten that computers store numbers in binary format. In the same way it is not possible to express the fraction 1/3 with a finite decimal places, analogously fraction well represented in the decimal base cannot be represented in binary, e.g. 1/10 is the infinitely repeating number:\n",
    "\n",
    "$0.0001100110011001100110011001100110011001100110011...$\n",
    "\n",
    "corresponding to $3602879701896397/2^{55}$ which is close to but not exactly equal to the true value of 1/10 (even though it is even printed to be like that!).\n",
    "Similarly 0.1 is not 1/10, and making calculations assuming that exactly typically yield to wrong results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1 0.1\n",
      "True\n",
      "True\n",
      "True\n",
      "True\n",
      "False\n",
      "True\n",
      "-2.7755575615628914e-17\n",
      "2.7755575615628914e-17\n",
      "False\n"
     ]
    }
   ],
   "source": [
    "print (1 / 10, 0.1)\n",
    "\n",
    "# is 1/10 the same of 0.1?\n",
    "print (1 / 10 == 0.1)\n",
    "\n",
    "# but then watch out! Does it work for 0.3, too?\n",
    "print(0.1 + 0.1 == 0.2)\n",
    "print(3/10 == 0.3)\n",
    "print(0.1 + 1/10 == 0.2)\n",
    "print(0.1 + 0.1 + 0.1 == 0.3)\n",
    "print(0.1 + 0.1 + 0.1 + 0.1 == 0.4)\n",
    "print(0.3 - 0.1 - 0.1 - 0.1)\n",
    "print(0.4 - 0.1 - 0.1 - 0.1 - 0.1 )\n",
    "print(0.3 == 0.1 + 0.1 + 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A lesson of paramount importance is that you must **never** compare floating point numbers with the \"==\" operator as *what is printed is not what is stored*!\n",
    "\n",
    "The function ```float.hex()``` yield the exact value stored for a floating point number:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.141592653589793\n",
      "0x1.921fb54442d18p+1\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "x = math.pi\n",
    "print(x)\n",
    "print(x.hex())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are ways to print floats (e.g. filling data into an output file) controlling the number of decimals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.1415926535898\n",
      "1.00000000000000022\n"
     ]
    }
   ],
   "source": [
    "print(format(math.pi, '.13f'))  # give 13 significant digits\n",
    "\n",
    "print('%.17f' % (0.1 * 0.1 * 100)) # give >15 significant digits\n",
    "\n",
    "# now repeat trying with >15 digits!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is no bug here, this is simply due to the fact that the mantissa is represented by a limited amount of bits, therefore calculations can only make sense if an appropriate number of decimal digits are concerned:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Single precision: 1.1920928955078125e-07\n",
      "Double precision: 1.1102230246251565e-16\n"
     ]
    }
   ],
   "source": [
    "# 23 bits are used for f in single precision floating point\n",
    "print(\"Single precision:\", 2**-23)\n",
    "\n",
    "# 53 bits are used for f in double precision floating point\n",
    "print(\"Double precision:\", 2**-53)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A typical case is subtraction of numbers very close by in value (e.g. when dealing with spectral frequencies). The same happens with functions evaluated near critical points (see later)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1 + 6.022e23 - 6.022e23"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Associative law does not necessarily hold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "print (6.022e23 - 6.022e23 + 1)\n",
    "print (1 + 6.022e23 - 6.022e23)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Distributive law does not hold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import math\n",
    "a = math.exp(1)\n",
    "b = math.pi\n",
    "c = math.sin(1)\n",
    "a*(b + c) == a*b + a*c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also identities after casting large numbers may not yield the expected result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "x = 287475839859383374\n",
    "print(x == int(float(x)))\n",
    "\n",
    "x = 2287475839859383\n",
    "print(x == int(float(x)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## From numbers to functions: conditioning and stability\n",
    "\n",
    "#### Function conditioning\n",
    "\n",
    "A mathematical function $f(x)$ is well-conditioned if $f(x+\\epsilon)\\simeq f(x)$ for all small perturbations $\\epsilon$.\n",
    "\n",
    "In other words, the function $f(x)$ is **well-conditioned** if the solution varies gradually as the input varies. For a well-conditioned function, small pertubations in the input result in small effects in the output. However, a poorly-conditioned problem only needs some small perturbation to have large effects. For example, inverting a nearly singluar matrix (a matrix whose determinant is close to zero) is a poorly conditioned problem.\n",
    "\n",
    "#### Algorithm stability\n",
    "\n",
    "Suppose we have a computer algorithm $g(x)$ that represents the mathematical function $f(x)$. $g(x)$ is **numerically stable** if $g(x+\\epsilon) \\simeq f(x)$ and it is called **unstable** if large changes in the output are produced. Analyzing an algorithm for stability is more complicated than determining the condition of an expression, even if the algorithm simply evaluates the expression. This is because an algorithm consists of many basic calculations and each one must be analyzed and, due to roundoff error, we must consider the possibility of small errors being introduced in every computed value.\n",
    "\n",
    "Numerically unstable algorithms tend to amplify approximation errors due to computer arithmetic over time. If we used an infinite precision numerical system, stable and unstable alorithms would have the same accuracy. However, as we see below (e.g. variance calculation), when using floating point numbers, algebrically equivalent algorithms can give different results.\n",
    "\n",
    "In general, we need both a well-conditioned problem and an algorihtm with sufficient numerical stabilty to obtain reliably accurate answers. In this case, we can be sure that $g(x) \\simeq f(x)$.\n",
    "\n",
    "In most of the cases, the solution to stability issues is solved by properly redefining the function as in the example above and below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Example of a poorly conditioned function: the tangent of an angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tan(x1) = 61249.008531503045\n",
      "tan(x2) = 158057.9134162482\n",
      "% change in x = 0.0006366263894271296 %\n",
      "% change in tan(x) = 158.05791343536947 %\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "# Define two numbers x and x + epsilon very close to pi/2\n",
    "x1 = 1.57078\n",
    "x2 = 1.57079\n",
    "# Calculate the tangent of the x1 and x2 angles\n",
    "t1 = math.tan(x1)\n",
    "t2 = math.tan(x2)\n",
    "\n",
    "print ('tan(x1) =', t1)\n",
    "print ('tan(x2) =', t2)\n",
    "print ('% change in x =', 100.0*(x2-x1)/x1, '%')\n",
    "print ('% change in tan(x) =', (100.0*(t2-t1)/t1), '%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Example of a numerically unstable algorithm: the limit      $\\lim_{x \\to 0} \\frac{1-\\cos(x)}{x^2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7f3d7143b160>"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEFCAYAAADzHRw3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1OUlEQVR4nO3deXyb1Zno8d8jyUvixLFjO5vtxNn3hSSErZAQUgjLAKXQAbrS9nKZW5gu0ym0nXa4bbm93RimLTSXS2kv0wKlBVqWQIACISyBLJCExCRxFi+JE2+xEy+xLencPyQ5sq3llfxKsuXn+/nwAUuvXh1eSY+OnvOcc8QYg1JKqaHPkeoGKKWUsocGdKWUShMa0JVSKk1oQFdKqTShAV0ppdKEK1VPXFhYaMrKylL19EoNbXv3+v49e3Zq26GSbtu2bQ3GmKJQ96UsoJeVlbF169ZUPb1SQ9uqVb5/v/56KluhUkBEKsPdpykXpZRKExrQlVIqTUQN6CLysIjUiciHYe4XEfmliFSIyE4RWWp/M5VSSkVjpYf+e2BthPsvB2b6/7kV+M3Am6WUUipWUQO6MeYNoCnCIdcAjxifzUCeiEy0q4FKKaWssSOHXgxUB/1d47+tHxG5VUS2isjW+vp6G55aKaVUgB0BXULcFnIJR2PMg8aY5caY5UVFIcsolVJKxcmOgF4DlAb9XQIcteG8Kg3VtnTwyp7jqW7GkLa96gRtne5UN0MNQnYE9GeAz/mrXc4FWowxtTacV6Whtfdt4suP6ISygbjn+XIqm9pT3Qw1CEWdKSoijwGrgEIRqQH+HcgAMMasA9YDVwAVQDtwS6Iaq4a+lo7uVDdhyHN7DZ1ub6qboQahqAHdGHNTlPsN8BXbWqSUiswYutxeDKEHsNTwpTNFlRqCjDG4Pbp9pOpNA7pKGt2/1l5dHk+qm6AGGQ3oKmlOtGv+3A6Br8UuzaOrPjSgq6SpbelIdRPSigZ01ZcGdJU0tc2nU92EtBDIXHV5NKCr3jSgq6TRHrq9tIeu+tKArpLmaIv20O2kAV31pQFdJU1ts/bQ7WD8w6I6uUj1pQFdJY320O3V5fFqKajqRQO6ShrNodsjEMO9XsMpXaRLBdGArpLC6zUc0x667fSaqmAa0FVSNLR20u0xjM3JTHVT0oLL6fvo1mpAV0E0oKukCOTPJ+VlA7oMwEAYA5n+gH5M01gqiAZ0lRSBCpdJY0akuCXpIcPl++ge1claKogGdJUUZ3rovoCuHfT4GXwf3AynQweaVS8a0FVS1DZ3kJ3hIG9kRqqbkjayMhzUnNCArs7QgK6SorblNJPGjEB0SwZ7CGS5nBrQVS8a0FVSHGnuYKJ/QBTOLAGrYhcYUM5yOTja3IHHq1dT+WhAV0lR29Lh66FrB902WS4Hbq+h7pQOjCofDegq4brcXupOdfYMiIKWLdohK8MJoGkX1UMDukq4I80dGAOlY0emuilpJctfulhzoj3FLVGDhQZ0lXDVTb6AU5o/QodEbdQT0Ju0h658NKCrhKv29yCDe+iacIlfIFvlEKFwVBZHdFli5acBXSVcVVM7mU4H43OzdVDUZiX5IzSHrnpoQFcJV9PUQXH+CJwOjeZ2CVxJX0DXHLry0YCuEq76RDsl+b3XcNEil/iZoIRVSf5Ijjafxqu16AoN6CoJqpvae/LnojkXW5Xkj6DL46W+tTPVTVGDgAZ0lVCnTndzor2b0vzeJYtGh0XjFvzrJvDLR9MuCjSgqwSr9pfUlY7VZXMTocT/RVmtpYsKDegqwXpKFvN1UlEilOT7llOobNQeutKArhIsMKlocp9ZojooGr/gS5ed4WRibjaVjW0pa48aPCwFdBFZKyJ7RaRCRO4Kcf8YEXlWRHaIyG4RucX+pqqhqOZEB6OyXD3roOuYqP2mFORwWAO6wkJAFxEncD9wOTAPuElE5vU57CvAHmPMYmAV8AsR0d2AFdVN7f60gEZyu/gWNjtzPcsKR2rKRQHWeugrgApjzEFjTBfwOHBNn2MMMFp8n9pRQBPgtrWlakiqbGrvlW7RDS7sN6Ugh8a2Lk6e7k51U1SKWQnoxUB10N81/tuC/RqYCxwFdgFfNcZ4bWmhGrI8XkNVYztTi3JS3ZT0E/S9WFbg+8Ks0l76sGcloIfqUvUd0roM+ACYBCwBfi0iuf1OJHKriGwVka319fUxNlUNNUebO+jyeJla0D+g66Bo/Ppeuin+66t5dGUloNcApUF/l+DriQe7BXjK+FQAh4A5fU9kjHnQGLPcGLO8qKgo3jarIeJggy/ATC08E9A1lW6/Kf4euubRlZWAvgWYKSJT/QOdNwLP9DmmCrgEQETGA7OBg3Y2VA09h0ME9ACdKToApvfP5pGZLsaNzuq53mr4ckU7wBjjFpHbgQ2AE3jYGLNbRG7z378O+CHwexHZhe+9dqcxpiGB7VZDwKGGNnIynRSNzuq5TTvoiVFWkKM9dBU9oAMYY9YD6/vcti7ov48Cl9rbNDXUHWpoo6wwR0sWk2BKwUg27tNxqeFOZ4qqhDnc2BYy3QI6KDoQoS5dWWEOdac6ae/SauHhTAO6Sogut5fqpvZ+AV0764kRGBg93KBpl+FMA7pKiOoT7XhN6AFR0D1FB8KE+HkzrXAUAAfqW5PdHDWIaEBXCXGoPnSFi84UTYxpRTmIaEAf7jSgq4QITHIJ10NX9srOcFKSP4KKOg3ow5kGdJUQBxvayBuZQd7I0Gu0hUobKGvCXbkZRaM4UK+16MOZBnSVEBXHW5k5blS/23VQNHGmF43iYH0rHt0wetjSgK5sZ4xhX90pZowbHf6YJLYn3RgTeoLWjHGj6HR7OXJCt6MbrjSgK9s1tnXR3N4dsoeuEmfGOK10Ge40oCvb7T/uCygzx2tAT6bpRb7rrQOjw5cGdGW7/XWnAJgZKeWiOZe4hVvYLD8nk4KcTO2hD2Ma0JXt9h9vZXSWi/G5Wf3u03VdEmt60SjtoQ9jGtCV7fbXnWLm+FGRg7f20OMW6dfN9HGjqKhv1bLQYUoDurLd/uOtYdMt2j9PrBnjRtHc3k19a2eqm6JSQAO6slVjayeNbV06IJpoYb4Z507wfZHuPXYqiY1Rg4UGdGWrQP52RpSSRd2xKH6RsimzNaAPaxrQla321QVKFsOkXDTnklAFo7IoGp1Fea0G9OFIA7qy1Ue1Jxmd7WLSmOyIx+mYXeLMmTCavcdPproZKgU0oCtbldeeZO7EXC1PTKE5E0az73grbo831U1RSaYBXdnG6zV8dOwU8ybmhj1Gw7w9Il3HORNy6XJ7OaybRg87GtCVbSqb2mnv8jB3YvgZogGacYlftBrzwMDoR8c07TLcaEAfxh555zBPbK227Xzltb4AMm/imLDHaCom8WaMG4XTIbZWumw+2MgPnt1j2/lUYmhAH6Y8XsMvXtrHX98/Yts5y2tP4nSIpRp0nckYv2hXLjvDybTCHFsrXV79qI6H3zpEZaNuoDGYaUAfpnYdaaGlo9vWapM9R08yrTCH7AynfSdVcZkzMbfnF5MdAl/Ab+yrt+2cyn4a0IepTf4Ppp0TfMprTzJvUvgBUdA69GRZMCmXI80dnGjrsuV8gS/+jfsabDmfSgwN6MPUG/v9Ad2meN7c3sXRltPMjVDhEkwTLvGz8potLPaNY+w60mLPc/r//c6BBrq1HHLQ0oA+DJ063c32qmbAvsC6x//zPlpA1w56csy3O6D73yhtXR62V56w5ZzKfhrQh6F3DjTi8RqyMxy2RfRdNb7AEegZRqNjovHzpckifzWOGZFBWcHIntfFjufMdDpwOoRN+zXtMlhpQB+GNu1vYGSmkwWTxuC1KbLurGmhJH8EY3MybTmfGrgFxWNs7aFnZThYUprXk65Tg48G9GHGGMMb++s5d1oBWRkO21IuO2qaWVySF/1AHRW1hZWruKhkDEeaO2iyYWDUGINDhItmFrHrSIst51T204A+zByob6WysZ2LZxchiC314I2tndSc6GBRibV0C+jyuQNh9SVbYGMe3eD7Lr5oViHGwCbtpQ9KlgK6iKwVkb0iUiEid4U5ZpWIfCAiu0Vko73NVHZ5pbwOgEvmjkfEnhR6IGAsstBD1/558vQE9JrmAZ/LGN9rt6gkj7E5mbz6Ud2Az6ns54p2gIg4gfuBjwM1wBYRecYYsyfomDzgAWCtMaZKRMYlqL1qgP5efpx5E3OZlDcCsGdwcmdNCyKwoNhayaLviQf+vMOVAUvfjLnZGUwtzLGph24QEZwOYfWccby0+xjdHi8ZTv2RP5hYeTVWABXGmIPGmC7gceCaPsfcDDxljKkCMMbo1/cgdKKti22VJ1gz1/d9K2JPymVnTTPTCnMYnZ0x4HMpey0oHmNLpYvXnPkOWTN3PCdPu9lyuGnA51X2shLQi4HgFZxq/LcFmwXki8jrIrJNRD5nVwPVGXc9uZM/D2Axrdf21uE1vnQL+D6gAw3nxhh21LRYGxBFx0STbXHJGI62nKbu5OkBnceYMwurXTSrkEyXg1f2xN9v21Z5gs/+9l3au9wDapfqzUpAD/UR7BsHXMAy4ErgMuB7IjKr34lEbhWRrSKytb5eB1VisffYKR7fUs07BxvjPsffy+soGp3VUyvukIGnXI6dPE39qc6YBkRBMy4DEctrtnRKPgDbqwY6Gcj0fBmPzHTxsRmFvFx+LO5feDuqm9m0v4H1u44NsF0qmJWAXgOUBv1dAhwNccyLxpg2Y0wD8AawuO+JjDEPGmOWG2OWFxUVxdvmYelPW3w983gDcJfby8Z99VwyZxwOh++TKSIDrjbZetgXKJZNGWvpeNFh0aSaPymXTJeDbQOc3WlM757dmrnjqW7qYL9/D9lYBeY/PP5e1YDapXqzEtC3ADNFZKqIZAI3As/0OeZvwIUi4hKRkcA5QLm9TR2+utxenn6/BiDuiUDvHWqitdPdk24Bf8plgF3lrYebGJnptLSphUq+LJeTxSVj2GpHQA+K6Jf4x2Fe3nM87vMBbK08QUWdbmhtl6gB3RjjBm4HNuAL0k8YY3aLyG0icpv/mHLgRWAn8B7wkDHmw8Q1e3h5pfw4J9q7Ad/gVDzWf1jLyEwnH5tR2HOb2JBy2Vp5giWlebhirHbQqf8DEdvFWzolnw+PtHC62zOAZzS9fl2Nz81mcckYXoozoAc6JiJnfn2qgbP0KTTGrDfGzDLGTDfG3OO/bZ0xZl3QMT8zxswzxiwwxtyXoPYOaR1d8X2gHt9SzcQx2UwtzImrh+72eNnw4TFWzxnHiMzgtcplQFP/WzvdlNeeZLk/T2uFDoom37LJ+XR7DB8OoHzRa/q/dpctmMCO6mZqTsS+d2mgY7JyVhFPbj9Clzv2FRxPd3t0o5Q+tIg0SV7YVcvCuzdwPMZqg8MNbbyxr54bz56MyxFfmeF7h5pobOviyoUTe93uGGBw/aCqGa+BZWXW8ufBdKZo/GJ9CwQGRgeSRzcGHH0ieuD99EIcA5uBjsRNKybT1NbFK+Wx9/RX//x1fv7S3pgfl840oCeB12u475X9uL0m5jUw/vhuJS6HcNOKUhwieONYivr5XbWMyHCyanbv+V4DTblsrWxCBM6anGf5MdpBT77CUVmUFYwcWEAP8QU8pSCHhcVjeG5Xbezn87/xVs4qYuKYbB6LY3D0+KlOHn7zMI2tnTE/Nl1pQE+CVz+qY+9x38CPJ4Yk+OluD09sreGyBRMYl5uNSOyDom6Plxc/PMbquX3TLb6Kk4H0lLdVnmD2+NHk6oSipIv1i3HplHy2V52IP0URIuUCcOWiieyobqa6Kba0S+BjkOF0cOPZk9m0v4ED9bFVzHi8ho5uDw+/dSimx6UzDegxevr9Gu5/rcLy8cYYfh10fCwB+dkdR2np6Oaz504BfD95Yx0UDZdugYH10N0eL9srT3B2HOkW0EHRgYjn0i2bkk9DaxeHG2PPdweeM2RA97+v1sfYSw98DhwCN58zmQyn8F/vVFp/vPfMoOr/e7uSFn/RgBXP76zl3jRN1Qz7gB7LyH/dydN89+kPYxqV37ivng+qm1njLxeMpYf+X5srmTluFOdM9QVNh4OYe1iBdMvFs/svrzOQxbk+OnaKti4Py8usD4gGnlMl37nTCgDf5ibxMMaEnENQOnYki0vG8NzO2AJ64G0sIhSNzuKqRZP4y7YaTp22Fpg9/hNcvmACrZ1u/t87hy0/9/pdtfzy1QrePmB9o46BVAgl07AO6JsPNjLney+y0eJO5ve+vI/2Lo/loGyM4T9e2U9x3ghuPNs3N8tqD31HdTM7a1r47HlTeqZcO0RiCsBdbi8vfHiMS0KkW8D3YYq3yiUQGM6ZWhDX47WDHr940ibTCnMYn5sV90xjrwk/iH7VoknsOtJCZWOb5fP51lc/8/fnzy+jtdPNk9tqLD0+8BmcP2kMa+aO4+G3DtHaaW0ZgcBjf/RcuaXPct2p0yy6+6UBLbuRLGkX0B97r4o/vhv9p5vb4+Wf/rANgAMWZrt9dOwkT2ytRsR6L/u1vXXsqG7mjtUzyMrwXWqrHfSH3zpETqaTT5x1ZtmcWAPwxn31NLV1cd3Svkvv+M8HcUfWtw80MK0ohwljsmN6nM4UTQ0R4bxpBbxzoDGuLwTDmbVc+rpy0URE4K/v951AHp63T9XMktI8lpTm8cg7lT3plIjt8R/idAhfuXgGze3d/GGztZSN2+v7MtlTe5Int0f/Amls7aLL4+WHz+2h7lT0KrU39zfwkxc/SklJZVoF9G2VJ/j2U7v47tPR5zT97q3DPZN1QvVegxljuOf5ckZlubh8wYSen3vRHnPfK/spHTuCTy4rwel/81r5Mjja3MFzO2u5acXkXisYOiS2iUVPba+hcFQmF84MvcyCxNjjD+j2eHnvUBPnT4+vd65S47zpBTS0dlIRx3R9X8oltEl5IzhvWgFPvV9jOYh5/TsgBfvC+WUcbGhjU0X0VIgnKAd/1uR8Vs4qYt3GA5y0kLLxGsO8SbmcNTmPn23YS1uUnr3b43uuk6fd/PC56BPgn99Vy29eP8Af3k3+sgZpE9DbOt1844kPACjJHxHx2KrGdn7x8l6WlOYBvm/sSDbsPsam/Q18bc0sxuZkWgrKfy+vY2dNC3esnkmG09GzfoqV3sfv/KP2t3xsaq/bHTEsd9vS3s3fy+u4enFx2DWrfVP/Yw/pO2uaaevycMH0wugHh6ETQuIX75U73/96vR1HHj3aGuyfXFpCZWO75SUGQk1UumLhRIpGZ/HQpoNRHx/4DAa+FP71stk0t3fz0BvRH+v2GlwOB9+7ah71pzpZt/FAlON9tcLLpuTz7I6jvBZlc49uj+/4e57fw8EYK3cGKm0C+g+f20NVUztFo7PIzgjf4zbG8J2nd+FyOPjxdQsBX/olnI4uDz98rpw5E0bzufOm4HI4Ih4PvqB978v7mFIwkuv8KROnP6BH692fPN3NY+9Vc9WiiRTn9f5icsRQtvjcrqN0ebxh0y0Q/6DoWxWNiJwZaIuJZlxSpnTsSIrzRsQ3MGoiv3RrF0xgZKbTcg7chOihZ7ocfPGCqWza3xB1VmugYxT4XC0oHsOVCyfy0JuHaIhSl+7xenE6hKWT87l68SQefONgxNmugS+Pr1w8nRnjRvFvf/0w4rK/bo+X/JEZZGc4+fqfPugJ8MmQFgH95T3HeXxLNbetnM550woiBtwntx/hzYoG7lw7u6cnH6nH/cDrFRxp7uB/Xj0fl9OB0xG9dPCZHUfZU3uSr6+Z1bPGicNiyuXx96po7XTz3y6c1u8+wfrEoqe2H2HW+FHMnxR+FyFHnIOib1U0MG9iLvk5mTE/NkA76PEbyLU7f3oBmw81WvqlGCxUiiRYTpaLtQsm8PzOWksVId4+g6IBnz53MqOzXFF7zYGOkTPoJN+4dBadbm/UsmKP1/Q87s7L5+AQ4e5n9oQ9vtufchmR4eLH1y3kSHMHv3hpX/jjvYb8nEzuuXYhO2paYipzHqghH9DrT3Vy15M7mT8p1x9ApecF6Kvu1Gl+9Pwelk3J59Pn+HrbQNjjDze08X82HuSaJZM4x98bdTqk5ydYKKe7Pfxsw14WFOdy9eJJPbcH3neRPozdHi+/e+sw500r6NkPMpjViUWHG9rYVnmC65aWhB3IgvhWW+zo8vB+VbPmz4eo86YX0NzezZ7akzE9ru9qi6Fcv7SEU51uSwt29R0UDcjNzuDmcyezfldtxKqZM4t7nTnH9KJRXL+0hD9urora4w6MaRXnjeBra2bySvnxsCtHBjphLqdwdtlYPnPuZB5+6xDvhqkYcnu8ZDgcXLloItedVcyvXq3gg+rmsO2x05AO6F6v4Zt/3kFrp5v7/nEJmS4HGQ5HyJ84xhju/MtOOro8/OSTC3E45EwaJEyA/sFze8hwCt+5Ym7PbU5H5F7yI+8c5khzB9+5fG5P3jzwON9zhY+gT79/hNqW09x6Uf/eOQRy6OGfO+CxLVU4HdKrQiakOCYWba1sosvj5fwZ8eXPNeNij3jr+QOrbVot1Q3ou9piKOdOK6A4b4Sl8j5v3wXWg3zpgqm4HA4ejJAPD3wGnX0uxFfXzATxlRiH4/YaXM4zj/vix6Yye/xo7n5md8hUSqADF/gMf/vyuZTmj+Sbf9kRckDV7Tlz/ruvmc+E3GzueGw7LR3WJz/Fa0gH9N9sPMDGffX821XzmDnetx53hktCDnI+9l41r+2t5861c5gxznesy/8ChTr+hV21vPpRHV9dM5PxuWdK85wSvofe3N7Fr1+t4OLZRf0CXk/KJUwEdXt8PxUXFo9h1ezQVSkOR/Qeepfby1+21rBm7rhe7Q4lnhLC1/fWk+l0sCLOGaJq4AYyoDwuN5v5k3LZuDfGgG6hh+5wCDcsL2HT/gaqosxIDbXYV3AbP7msmD9vqwlbJngm5dL79kl5I7jlgjKe2n6EnTXNIR/r9fZOH2U4HfzoEws40tzBL//ePz0SqHLJ8P+iz8ly8YtPLabmRAf3rO9f9dIVtHl2bnYGv7zpLGqbT3PXkzsTXgwwZAP6uwcb+cVLe7lq0UQ+c87knttdIXrolY1t/Oj5PVwwo4AvnF/Wc7vDITjkzAsW0NLezfef2c28ibncckHvSpNADj3UC/OrVyto7XRz1+Vz+93njFLl8rcPjlLZ2M4/XzIzbJrESs77pT3HaGzr4qYVkyMeB4Gp/7G9wV7bW8c508aSk+WK6XFnnlP76Km2anYR26pOxNRjtPouufHsyTgdwqNRFtsKl0MPuPWi6bg9Xv5vmF66t0+VS7DbL55B4ahMfvDsnpDvb1+VS+/HnV02lhuWlfDQpoPsO36q3/HQO19/dtlYbr1wGo++W8Xre3tXvbg9hoygXwDLpuRz59o5vPDhMR6JYXmDeAzJgN7Q2skdj73PlIIcfnzdwl5BIsMpvQK0x2v4xhM7cDqEn12/uFcaBMDldPTrod+zfg9NbV389PpF/Ur+XGFSJwfrW3nkncPcsKyU2RP6794TqYfu8frWe5k7MZc1c/tP0Q8QC2u5PPpuFSX5I7goTO157zbFVuVS2djGwfq2kMsIxEoHReM30Et38exxeLyGN/dbn/oeqUcdbMKYbC6ZM44/b62OuMZ5tEHWqYU5XHtWMY+8Uxlyg2tviEHRgNHZGXzz0tlsrTwRckmC4EHRYN++Yi6js13c9eTOXp/v4Bx6sK9/fBazxo/izid30tx+ZhVVt9fbMz4X8OULp3LJnHHc83x52F8OdhhyAd3jNXz9Tx/Q0tHN/Tcv7TXxBgIB+swbad3GA2yrPMEPrpnPpLz+9ekuh/SqinmrooEnttbw3y6cFnJg0hEiTWOM4e5n95DtcvLNy2aHbHfgeyFUQH52x1EONbTx1UtmROzBOqL0qA81tPH2gUZuWjG53xdXKBLjBheB+tvVcwYe0FXqLCnNIzfb1a9nGYkxxnLe/tPnTqGxrYsNu8Ovk+6rQ498wn9ePRO31/DA6/0rXvrWofd1w/JS5k3M5X+/8FG/qptwAX1sTiZ3Xz2f7VXNPPzmmRUcA/Gkb68+O8PJvZ9aQmNrF3c9uavns9ntMf2Cv4jw8xsWUzgqk688mrh8+pAL6H/eWs2m/Q3cffV85oUoyctw+KpcjDFsq2zi3pf3ceXCiVy7JPQAoa9qxfdCdHR5+PZTu5hamMPX1swMeXzgRQ0OhC/tOc4b++r5+sdnUTQ6K+TjAm+8vikXt8fLL1/dz+zxo7l03oSI/+/RVlt87L0qXA7hhmUlEc8TEOtqi6/urWdaYQ5lhTnWH9T3OeN+pLKLy+ngollFvL6v3nL5YrjVFkO5cEYhpWNHRFyCo+9aLqGUFeZw/dISHn2vitqWjl739azWGOYkTofw/X+Yx5Hmjn6Dq+ECOsDViydx6bzx/PylvT3L+QZ+8fftdYOv/v1ba2fz4u5jPUsPuL3ekJP58nMy+dXNS6ltPs1PX/wo0v963IZcQL9uaQn3/eOSnsWu+gpcyIbWLu549H0m5WXz408uDNsbyHA6er7tf7rhI6qa2vnxdQvDTk5y9umhn+728INn9zB7vG/iUTjhqlye3F7Dwfo2vrZmZtRedaSJRW2dbh5/r4pL549nXJTB0IBYJha1d7nZfLCx3yYZ8dIdiwbAhku3avY46k91Wi5fDLfaYigOh3DziilsPtjULx8d4PVaS+HcvnoGXq/hgdd699I9Yapcgp07rYArF07k/tcqepVAekz4gC4i/OgTC8jOcPKtv/hSL+FSLgFf/tg0Vs0u4ofPl7P7aAvd7t459GDLpuTzm88s41/D/JIfqCEX0DNdDq49qzhsgA5M5PnGEx9Q39rJ/TcvjbgBQ6CufNP+en731mG+cH5ZxBmQfQc3H3j9gG/i0TXzI26UHCqHfrrbw3+8vJ8lpXmsXRC5dw6Rc+hPba/h5Gk3X+qzXECUM1ruob9d0UiX2zvgdIuOiQ4OK2cVIYLlrd9i6aED/OPZpWRnOHqlLoJFGxQNKB07kk+dXcrjW3rXlnt6BiojP/57V80jw+nge3/b3ZMScXvCB3SAcaOzufvqeWyrPMHv3jpEd5iUS4DDIfzihsXkjcjgjkffp6WjO2Is+Pi88eSNjH9SXiRDLqBHE/hm3LS/ge9cMZdFJXkRj3c5hIbWLr755x1ML8rhrsvnRD0efD30A/WtrNt4gH9YPCnqNPhQa7n8/u3DHDt5mrsun2Op+iNcDt3rNfzurcMsLs1j6eT8qOcJ8D2ltYj+SvlxRmW5OHuq9fNHooOi8fNPqRnQOYpGZ7F8Sj4bdlsM6FGm/vc1NieTTy4t4an3j1B/qv9UfCs59IDbL/aNLd0bNDsz1MSiUCaMyeZfLp3FG/vqewZIPSGqXPq6dkkxa+aO46cb9lLu/xUT6UugYFQW9924hEONbRw7eZqMgW7YG6e0C+iBF+qy+eN7lSiGPd4pvLznOI2tXfznjWdFXAcGggZFPYZvP7mLbJeD713Vv0yxr8BPw0A8b2nv5oHXKlg1u8jymijhyhZf31fHwYY2vnhBWUxlgVZXb3R7vGzYfYzVc8aR5Yp8fdTQcdn8CZTXnrS0jrnXmJhLTr/4sal0ub0hl7U1xhAiJR3SpLwRfOljU3nq/SPsqmnpaQ9ETrkEfO68MhYWj+EHz+2hpaM7YsolQET48XWLGJ3l4g+bfSWYkXrd4Fv87J9X+8beDjVYXxveTmkX0FdM9eXNfvrJxZbegIGBjq9/fFbIqpb+x/vO+V+bD/Pe4Sb+7cp5jBsdPWcdePMGUi4PbKzgVKebb10W+RdBr3OESbk8/OZhJuRmc0WIbeYiEayt3vjeoSZOtHdzuYW0UNTn1JTLoHHZfN/rGakaJVisr930olFcMmccf9hc2a/SJFrZYl//Y9V0CnIy+dHzvtryMymX6OdwOoT/9YmFNLZ28vMNeyMOigYrGp3FT69f1PN3tF49wD9fMpPPnDuZf1o1PeqxiZB2AX3epFzu//RSxoy0tnFx4ahMVpSN5baV1l6AwJtw3caDXDCjgBuWW6socQZVuVQ3tfP7tw5z7ZLikJU64YRay2X30RberGjgc+dPCbtMbqTzWcl8vPDhMUZkOG0bEAXdsWgg7JptWDp2JAuKc3nxw+gBPdaUS8CXL5xGY1sXf33/SK/bw63lEs7o7Ay+tmYm7x5q4pXyuqhVLn0tLBnD584r4w/vVtLS0W2pZw9wydzxfP68KYzKcpHpiv75cjqEH127kLULYutc2SXtAnqsfn/LCv7ryyssfWPDmZHuTKeDH39ikeWfocGrLd7zfDkOEb61NraR7lBrudz/WgWjs1x8+pzwFTbhCNFz2V6v4cXdx1g1uyjqRiDWnlO76IPJ2vkT2F7VzPEQk3eCGWJPuQCcO20s8yfl8n83Hew1fuSNoa494KYVk5lelMOP15fT2e0bqIwlVf2vl82mOG+Ev4duPfTdffV83rzz4pg7TKkw+FuYYDlZrpjywoFj/+XSWUwuGGn5cYGexKb99by4+xi3r57BxDGRN+Lod44+PfSKulO88OExPn9+GWNGWPtFEkwsbJixreoE9ac6LVXhxEI3uIifwb56/sDrGi3tEm8PXUS4beV0DtS38ULQLwGrM0+DuZwOvnPFXA42tPVsCm21pw2+z/pPPulLoVjpbQeISMKqUuw27AN6rFbOKuI/b1zSb42XaAK/AF7bW8/ksSNjLC/06Tso+sDrB8h2OfliHOcCaymX9btqyXQ5dHZompoxbjSzxo/imQ8i7wcaTwAOuGLhRKYX5fCrV/f39NKtli32tXrOOC6aVdRTnWM15RJwwYxCfvv55RHnjAxlGtBjlJPl4polxZZTNAHBPYnvXzUvajVNKCJnlu6tamznbx8c5eZzJjM2zo0mJEodutvj5bmdtayaVdRviYV46aDo4HPNkmK2Vp6guin8ComRlruNxukQ7lg9k4+OnepZKz3WQdEAEeHuf5jXU54c6+cQfHnxUMuApAMN6Eni9K/suHJWEZdEWIArkuA69HVvHMApEnbtdCuirbb41oFG6k91RtzGLl6acImf3dmqa5b4NmJ5Zkf4XvpA0zxXLZpIWcFIfvn3/RhjYqpD72ta0Si+7N/RK3MI5LWTSa9GkmS6HDz0+eXc+ylr5ZShBMoWqxrb+fPWaj51dknUNc8jESIH1qe31zBmRAYXa7olrZXkj+Tssnyefv9I+C94C+uhR+JyOvjKxTPYU3uSv5fXWVrLJZKvrZnJr246i4UWSo2HEw3oSbR6zngKRoVevMuKwAYX//HKvp6fsQMRaXGutk43G3Yf58pFExMymUjHRONnou3YHIdrlhRTUdfK7qOh13axsmNRNNeeVczksSO59+V9uL3xpVwCslxO/mHxpJhz6OlOA/oQIiKcaO/irx8c4QvnTx1Q7xz8ZZBh+ugbdh+jo9sTfRs7lRauXDiRDKf0qxcPMAbLMzvDyXA6+MbHZ7Gn9iSbDzYOqIeuQrP0EonIWhHZKyIVInJXhOPOFhGPiFxvXxNVgEN8ay2PynLxTxYnQkUUYer/U9uPUJI/guVT7Fm7pecpdVR0UMrPyWT1nHE8/f6RkBtTeGNYbTGSqxdPYu7EXE53e/W9kABRA7qIOIH7gcuBecBNIjIvzHE/ATbY3UjlE/iJetvK6ZZnwkYihK5brGps582KBq5fVpLAD53mXOKVqHTVjSsm09jWFXIFxlhXWwzH4TgzoU576Paz0kNfAVQYYw4aY7qAx4FrQhx3B/AkYH0bFBWTSXkjKB3r2wTXDr469P7R4fEtVTjEtwSq3fQzPHhdNLOI4rwRPBZiP1A7v0RWzSpi1ewiJqZp6WAqWdnptxioDvq7Bjgn+AARKQY+AawGzg53IhG5FbgVYPLk6JsYq95uWzmdL31sqm1TkENN/e/2eHliaw2r54yPeSZrLHRQNH6JunROh3DD8hL+8+/7qW5qp3TsmZnQvh66PV/HIsJDn1s+oEFRFZqVyBDqqvd9T90H3GmM8YQ49syDjHnQGLPcGLO8qCj6JsaqPzvXk/ANivb2yp7jNLR2cvM59vfO1eB3w3Lf6/6nLdW97xhgmWFfLqdDK1QSwEp0qAGCP90lQN8ZCMuBx0XkMHA98ICIXGtHA1XihJpY9Oh7VUwak83KWYmpPddOmT0SdRmL80awalYRj2+pptN9pn9m5/oxKnGsBPQtwEwRmSoimcCNwDPBBxhjphpjyowxZcBfgP9hjPmr3Y1V9hJ6V7lU1J1i0/4GbloxOa4p1bHQjMsAJPji3XLBVBpaO3luR23PbfFscKGSL2pAN8a4gdvxVa+UA08YY3aLyG0icluiG6gSqM8H9LdvHibL5eDmcxI3vqHL5w5+F84sZMa4Ufz2zUM9v+AGsJSLSiIrg6IYY9YD6/vcti7MsV8YeLNUMgQ+oMYYTrR389T2Gq5bWjyg2axW6aBo/MJNBrOLiPDFC6bynad38e6hJs6dVuAL6BrRBz2dKTqMBaoMjIE/bq6k0+3lizEuC6zS03VLi8kfmcHDbx4C7K1yUYmjAX0YC3w+T7s9PLK5kotmFTFz/OikPKca3LIznNx8zmReLj/O4YY2jDGachkCNKAPY4EP6J+2VFN/qpP/PoCleGOV6LRBOktWuurz55WR4XTwm9cPaMpliNCAPowFPqDrNh5g2ZR8zp9ekPjnTPgzKLuMy83mprNLeXJ7DfWtnTqgPQRoQB/GAjnR4yc7uWP1DM2Rqn7++8rpiEBTW5f20IcADeiKxSVjWDkruTN3tcolfsm8dJPyRvTMHtWAPvhpQB/GApOH7lg9M2m9cw0K9kjmZfynldNxOUR/wQ0BlurQVXq6fMEEHELce5wOhPbQ4xdpH9hEKB07kp/fsDhtN1ZOJxrQh7EpBTncepENG2WotHet7lw1JGjKRSWZ/my3hV5GFYIGdJUSWoceP71yKhwN6CqpdFxNqcTRgK5SQgdF46fXToWjAV0ppdKEBnSVVJpxsYteSdWfBnSllEoTGtBVUulsQ6USRwO6Sgkd2ItPsmeJqqFFA7pSSqUJDegqqTThYg+9jioUDegqJXSmaHw046Ii0YCukkrHRJVKHA3oKiW0pxkfvWwqEg3oSimVJjSgq6TSlItSiaMBXaWEpg7io3XoKhIN6CqpRAvulEoYDegqJbSnGR+9aioSDehKKZUmNKCr5NKMiy10cFmFogFdpYSmDuKjmSoViaWALiJrRWSviFSIyF0h7v+0iOz0//O2iCy2v6kqHWjHUqnEiRrQRcQJ3A9cDswDbhKReX0OOwSsNMYsAn4IPGh3Q1V60Z5mfHQNHBWJlR76CqDCGHPQGNMFPA5cE3yAMeZtY8wJ/5+bgRJ7m6mUUioaKwG9GKgO+rvGf1s4XwJeCHWHiNwqIltFZGt9fb31Vqq0oTsWKZU4VgJ6qE9gyN99InIxvoB+Z6j7jTEPGmOWG2OWFxUVWW+lSkOaOoiHpqpUJC4Lx9QApUF/lwBH+x4kIouAh4DLjTGN9jRPpRvtnyuVOFZ66FuAmSIyVUQygRuBZ4IPEJHJwFPAZ40x++xvpko32tNUyn5Re+jGGLeI3A5sAJzAw8aY3SJym//+dcD3gQLgAX+O1G2MWZ64ZiullOrLSsoFY8x6YH2f29YF/feXgS/b2zSVjnRMVKnE0ZmiKiU04xIfTVWpSDSgq6TS5XOVShwN6ColtKcZn8BMUf1aVKFoQFdKqTShAV0llQ6KKpU4GtBVSuiORfHRy6Yi0YCulFJpQgO6SqpAxkU7mvHpuW6aulIhaEBXSqk0oQFdJZf2LG2iF1L1pwFdpYQO7sVHB5NVJBrQlVIqTWhAV0kVmPqve2PGJ3DVNOGiQtGArpRSaUIDukoqnSmqVOJoQFepoRmXuOiYqIpEA7pSSqUJDegqqTTjolTiaEBXKaGZgzjphVMRaEBXSSU6KqpUwmhAVymhg3vx0fp9FYkGdKWUShMa0FVSacZFqcTRgK5SQlMH8dFUlYpEA7pKKu2g20OvowpFA7pKCe1pxkcvm4pEA7pSSqUJDegqqXRQ1CZ6HVUIGtBVSmjqID66Y5GKRAO6SjLtWiqVKBrQVUpoTzM+etVUJJYCuoisFZG9IlIhIneFuF9E5Jf++3eKyFL7m6qUUiqSqAFdRJzA/cDlwDzgJhGZ1+ewy4GZ/n9uBX5jcztVmtBBUaUSx2XhmBVAhTHmIICIPA5cA+wJOuYa4BHj+x29WUTyRGSiMabW9hartPCdp3aRk2Xl7aeCebyadFHhWflEFQPVQX/XAOdYOKYY6BXQReRWfD14Jk+eHGtbVRqYM2E0/7i8lFOd3aluypC1sGQMeZsyUt0MNQhZCeihfiT37SZYOQZjzIPAgwDLly/XrsYwNDLTxU+uX5TqZgx965ypboEahKwMitYApUF/lwBH4zhGKaVUAlkJ6FuAmSIyVUQygRuBZ/oc8wzwOX+1y7lAi+bPlVIquaKmXIwxbhG5HdgAOIGHjTG7ReQ2//3rgPXAFUAF0A7ckrgmK6WUCsVSmYExZj2+oB1827qg/zbAV+xtmlJKqVjoTFGllEoTGtCVUipNaEBXSqk0oQFdKaXShKRq1TsRqQcq43x4IdBgY3PsMljbBYO3bdqu2Gi7YpOO7ZpijCkKdUfKAvpAiMhWY8zyVLejr8HaLhi8bdN2xUbbFZvh1i5NuSilVJrQgK6UUmliqAb0B1PdgDAGa7tg8LZN2xUbbVdshlW7hmQOXSmlVH9DtYeulFKqDw3oSimVJoZ8QBeRb4qIEZHCVLcFQER+6N8o+wMReUlEJqW6TQAi8jMR+cjftqdFJC/VbQIQkRtEZLeIeEUk5eVl0TZETxUReVhE6kTkw1S3JUBESkXkNREp97+GX011mwBEJFtE3hORHf52/c9UtymYiDhF5H0Rec7ucw/pgC4ipcDHgapUtyXIz4wxi4wxS4DngO+nuD0BLwMLjDGLgH3At1PcnoAPgeuAN1LdEIsboqfK74G1qW5EH27gX4wxc4Fzga8MkuvVCaw2xiwGlgBr/fs0DBZfBcoTceIhHdCB/wC+RYjt7lLFGHMy6M8cBknbjDEvGWPc/j8349tVKuWMMeXGmL2pbodfz4boxpguILAhesoZY94AmlLdjmDGmFpjzHb/f5/CF6SKU9sq33LexphW/58Z/n8GxedQREqAK4GHEnH+IRvQReRq4IgxZkeq29KXiNwjItXApxk8PfRgXwReSHUjBqFwm52rKESkDDgLeDfFTQF60hofAHXAy8aYQdEu4D58nVBvIk5uaYOLVBGRV4AJIe76LvAd4NLktsgnUruMMX8zxnwX+K6IfBu4Hfj3wdAu/zHfxfdT+Y/JaJPVdg0SljY7V72JyCjgSeBrfX6hpowxxgMs8Y8VPS0iC4wxKR1/EJGrgDpjzDYRWZWI5xjUAd0YsybU7SKyEJgK7BAR8KUPtovICmPMsVS1K4RHgedJUkCP1i4R+TxwFXCJSeIEhBiuV6rpZucxEpEMfMH8j8aYp1Ldnr6MMc0i8jq+8YdUDyhfAFwtIlcA2UCuiPzBGPMZu55gSKZcjDG7jDHjjDFlxpgyfB/EpckI5tGIyMygP68GPkpVW4KJyFrgTuBqY0x7qtszSFnZEF35ia839Vug3Bhzb6rbEyAiRYEqLhEZAaxhEHwOjTHfNsaU+GPWjcCrdgZzGKIBfZD73yLyoYjsxJcSGhSlXMCvgdHAy/6SynXRHpAMIvIJEakBzgOeF5ENqWqLf9A4sCF6OfCEMWZ3qtoTTEQeA94BZotIjYh8KdVtwtfj/Cyw2v+e+sDf+0y1icBr/s/gFnw5dNtLBAcjnfqvlFJpQnvoSimVJjSgK6VUmtCArpRSaUIDulJKpQkN6EopZQO7F1ATkZ/6FxcrF5Ff+stEI9KArpRS9vg9Ni2gJiLn4ysLXQQsAM4GVkZ7nAZ0pZSyQagF1ERkuoi8KCLbRGSTiMyxejp8s0kzgSx8C4wdj/YgDehKKZU4DwJ3GGOWAd8EHrDyIGPMO8BrQK3/nw3GmKhL7g7qtVyUUmqo8i9adj7w56D0d5b/vuuAH4R42BFjzGUiMgOYy5llrl8WkYv8vwLC0oCulFKJ4QCa/Zvd9OJfyCzSYmafADYH1nUXkRfwbSISMaBrykUppRLAv5TwIRG5AXyLmYnIYosPrwJWiojLv6LlSizscqQBXSmlbBBmAbVPA18SkR3AbqzvgPUX4ACwC9gB7DDGPBu1Dbo4l1JKpQftoSulVJrQgK6UUmlCA7pSSqUJDehKKZUmNKArpVSa0ICulFJpQgO6Ukqlif8PmBoShhThxHwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Catastrophic cancellation occurs when subtracitng\n",
    "# two numbers that are very close to one another\n",
    "# Here is another example\n",
    "\n",
    "# We'll see numpy and matplotlib in the next lectures: forget about the technical details, for now\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def f(x):\n",
    "    return (1 - np.cos(x))/(x*x)\n",
    "\n",
    "#x = np.linspace(-4e-1, 4e-1, 1000)\n",
    "x = np.linspace(-4e-8, 4e-8, 1000)\n",
    "plt.plot(x, f(x))\n",
    "plt.axvline(1.1e-8, color='red')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.999999999999999888977697537484\n",
      "0.000000000000000111022302462516\n",
      "0.917540\n"
     ]
    }
   ],
   "source": [
    "# We know from L'Hopital's rule that the answer is 0.5 at 0\n",
    "# and should be very close to 0.5 throughout this tiny interval\n",
    "# but errors arisee due to catastrophic cancellation\n",
    "\n",
    "print('%.30f' % np.cos(1.1e-8))\n",
    "print('%.30f' % (1 - np.cos(1.1e-8))) # failure point: the exact answer is 6.05e-17\n",
    "print('%2f' % ((1 - np.cos(1.1e-8))/(1.1e-8*1.1e-8)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solution: rewrite the function using $\\sin$ instead of $\\cos$: $1-\\cos(x)$ = $2 \\sin^2 (\\frac{x}{2})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7f3d37e88f70>"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEFCAYAAADt1CyEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQFUlEQVR4nO3df6zdd13H8eeLtuM30dBLJv2xNtpMymQDL5WBwgJMuzk3R0bS8SMSiXWE4jAQGSzZAsZExCDBjtRGlv0BriKDOUZnmXHITAa2Ky20lGodQu86bQHZaCDOO97+cU/H2eG093vbe3tuPzwfyc3O93u+33Pet9t97tvvOed7U1VIktr1pFEPIEmaW4Zekhpn6CWpcYZekhpn6CWpcQtHPcAwixcvrhUrVox6DOnMtH//1D/PPXe0c+i0uv/++79dVWPD7puXoV+xYgU7duwY9RjSmemii6b++fnPj3IKnWZJvnm8+zx1I0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mN6xT6JGuT7E9yIMl1Q+6/KMnDSXb1vm7orV+W5J4k+5LsTXLtbH8DkqQTWzjdBkkWADcBFwMTwPYkd1TV1wY2vbeqLhtYNwm8o6p2JnkmcH+Su4fsK0maI12O6NcAB6rqgap6FNgCXNHlwavqoara2bv9fWAfsORkh5UkzVyX0C8BDvYtTzA81hcm2Z3kriTPH7wzyQrghcCXhj1JkvVJdiTZceTIkQ5jSZK66BL6DFlXA8s7gXOq6nzgL4Hbn/AAyTOA24C3V9Ujw56kqjZX1XhVjY+NjXUYS5LURZfQTwDL+paXAof6N6iqR6rqaO/2VmBRksUASRYxFfmPV9WnZmVqSVJnXUK/HViVZGWSs4B1wB39GyQ5O0l6t9f0Hvc7vXUfBfZV1Qdnd3RJUhfTvuumqiaTbAC2AQuAm6tqb5JrevdvAq4C3pJkEvghsK6qKsmvAm8EvppkV+8h39M76pcknQbThh4ePx2zdWDdpr7bG4GNQ/b7F4af45cknSZ+MlaSGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxnUKfZG2S/UkOJLluyP0XJXk4ya7e1w19992c5HCSPbM5uCSpm2lDn2QBcBNwCbAauDrJ6iGb3ltVF/S+3te3/hZg7WwMK0mauS5H9GuAA1X1QFU9CmwBruj6BFX1BeC7JzmfJOkUdQn9EuBg3/JEb92gC5PsTnJXkufPynSSpFO2sMM2GbKuBpZ3AudU1dEklwK3A6tmMkiS9cB6gOXLl89kV0nSCXQ5op8AlvUtLwUO9W9QVY9U1dHe7a3AoiSLZzJIVW2uqvGqGh8bG5vJrpKkE+gS+u3AqiQrk5wFrAPu6N8gydlJ0ru9pve435ntYSVJMzdt6KtqEtgAbAP2AZ+oqr1JrklyTW+zq4A9SXYDHwbWVVUBJLkVuA84N8lEkjfPxTciSRquyzn6Y6djtg6s29R3eyOw8Tj7Xn0qA0qSTo2fjJWkxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxhl6SWqcoZekxnUKfZK1SfYnOZDkuiH3X5Tk4SS7el83dN1XkjS3Fk63QZIFwE3AxcAEsD3JHVX1tYFN762qy05yX0nSHJk29MAa4EBVPQCQZAtwBdAl1qey74y99zN7+dqhR+bioaUzxg29n4H3/dV9I55EM7X6uc/ixt96/qw/bpdTN0uAg33LE711gy5MsjvJXUmOTdp1X5KsT7IjyY4jR450GEuS1EWXI/oMWVcDyzuBc6rqaJJLgduBVR33nVpZtRnYDDA+Pj50m+nMxf8JpTPOrc8C4G9//8IRD6L5ossR/QSwrG95KXCof4OqeqSqjvZubwUWJVncZV9J0tzqEvrtwKokK5OcBawD7ujfIMnZSdK7vab3uN/psq8kaW5Ne+qmqiaTbAC2AQuAm6tqb5JrevdvAq4C3pJkEvghsK6qChi67xx9L5KkIbqcoz92OmbrwLpNfbc3Ahu77itJOn38ZKwkNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjDL0kNc7QS1LjOoU+ydok+5McSHLdCbZ7cZLHklzVt+7aJHuS7E3y9lmYWZI0A9OGPskC4CbgEmA1cHWS1cfZ7v3Atr515wG/B6wBzgcuS7JqdkaXJHXR5Yh+DXCgqh6oqkeBLcAVQ7Z7G3AbcLhv3fOAL1bVD6pqEvhn4MpTnFmSNANdQr8EONi3PNFb97gkS5gK+KaBffcAL0/y7CRPAy4Flp38uJKkmVrYYZsMWVcDyx8C3lVVjyU/3ryq9iV5P3A3cBTYDUwOfZJkPbAeYPny5R3GkiR10eWIfoInHoUvBQ4NbDMObEnyn8BVwEeS/DZAVX20ql5UVS8Hvgv8+7AnqarNVTVeVeNjY2Mz+y4kScfV5Yh+O7AqyUrgQWAd8Lr+Dapq5bHbSW4B7qyq23vLz6mqw0mWA68BLpyd0SVJXUwb+qqaTLKBqXfTLABurqq9Sa7p3T94Xn7QbUmeDfwf8Naq+p9THVqS1F2XI3qqaiuwdWDd0MBX1ZsGln/tZIeTJJ06PxkrSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY0z9JLUOEMvSY3rFPoka5PsT3IgyXUn2O7FSR5LclXfuj9MsjfJniS3JnnKbAwuSepm2tAnWQDcBFwCrAauTrL6ONu9H9jWt24J8AfAeFWdBywA1s3O6JKkLroc0a8BDlTVA1X1KLAFuGLIdm8DbgMOD6xfCDw1yULgacChU5hXkjRDXUK/BDjYtzzRW/e43pH7lcCm/vVV9SDw58C3gIeAh6vqc8OeJMn6JDuS7Dhy5Ej370CSdEJdQp8h62pg+UPAu6rqsSfsmPwsU0f/K4HnAk9P8oZhT1JVm6tqvKrGx8bGOowlSepiYYdtJoBlfctL+cnTL+PAliQAi4FLk0wCi4BvVNURgCSfAl4KfOwU55YkddQl9NuBVUlWAg8y9WLq6/o3qKqVx24nuQW4s6puT/IrwEuSPA34IfAqYMcszS5J6mDa0FfVZJINTL2bZgFwc1XtTXJN7/5NJ9j3S0k+CewEJoEvA5tnZXJJUiddjuipqq3A1oF1QwNfVW8aWL4RuPEk55MknSI/GStJjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktS4VNWoZ/gJSY4A3zzJ3RcD357FcWaLc82Mc82Mc81Mi3OdU1Vjw+6Yl6E/FUl2VNX4qOcY5Fwz41wz41wz89M2l6duJKlxhl6SGtdi6DePeoDjcK6Zca6Zca6Z+amaq7lz9JKkJ2rxiF6S1MfQS1Ljmg59kncmqSSLRz0LQJI/TvKVJLuSfC7Jc0c9E0CSDyT5em+2Tyf5mVHPBJDktUn2JvlRkpG+FS7J2iT7kxxIct0oZ+mX5OYkh5PsGfUs/ZIsS3JPkn29f4fXjnomgCRPSfKvSXb35nrvqGc6JsmCJF9OcudsP3azoU+yDLgY+NaoZ+nzgap6QVVdANwJ3DDieY65Gzivql4A/Bvw7hHPc8we4DXAF0Y5RJIFwE3AJcBq4Ookq0c5U59bgLWjHmKISeAdVfU84CXAW+fJn9n/Aq+sqvOBC4C1SV4y2pEedy2wby4euNnQA38B/BEwb15trqpH+hafzjyZrao+V1WTvcUvAktHOc8xVbWvqvaPeg5gDXCgqh6oqkeBLcAVI54JgKr6AvDdUc8xqKoeqqqdvdvfZypgS0Y7FdSUo73FRb2vkf8cJlkK/Cbw13Px+E2GPsnlwINVtXvUswxK8idJDgKvZ/4c0ff7XeCuUQ8xzywBDvYtTzAPonWmSLICeCHwpRGPAjx+imQXcBi4u6rmw1wfYurA9Edz8eAL5+JBT4ck/wicPeSu64H3AL9+eieacqK5qurvq+p64Pok7wY2ADfOh7l621zP1F+5P346Zuo61zyQIetGfhR4JkjyDOA24O0Df6Mdmap6DLig91rUp5OcV1Uje40jyWXA4aq6P8lFc/EcZ2zoq+rVw9Yn+SVgJbA7CUydhtiZZE1V/deo5hrib4DPcppCP91cSX4HuAx4VZ3GD1fM4M9rlCaAZX3LS4FDI5rljJFkEVOR/3hVfWrU8wyqqu8l+TxTr3GM8sXslwGXJ7kUeArwrCQfq6o3zNYTNHfqpqq+WlXPqaoVVbWCqR/SF52OyE8nyaq+xcuBr49qln5J1gLvAi6vqh+Mep55aDuwKsnKJGcB64A7RjzTvJapo6yPAvuq6oOjnueYJGPH3lWW5KnAqxnxz2FVvbuqlvZ6tQ74p9mMPDQY+nnuT5PsSfIVpk4tzYu3nAEbgWcCd/fe+rlp1AMBJLkyyQRwIfDZJNtGMUfvheoNwDamXlT8RFXtHcUsg5LcCtwHnJtkIsmbRz1Tz8uANwKv7P03tat3xDpqPwfc0/sZ3M7UOfpZfzvjfOMlECSpcR7RS1LjDL0kNc7QS1LjDL0kNc7QS9Icmu0LzyX5s94F2fYl+XDvrawnZOglaW7dwixdeC7JS5l66+oLgPOAFwOvmG4/Qy9Jc2jYheeS/HySf0hyf5J7k/xi14dj6tOzZwFPZuqibP893U6GXpJOv83A26rql4F3Ah/pslNV3QfcAzzU+9pWVdNe2viMvdaNJJ2Jehd6eynwd32n15/cu+81wPuG7PZgVf1Gkl8AnsePLyV+d5KX9/7WcFyGXpJOrycB3+v9AqIn6F387UQXgLsS+OKxa+onuYupX+xywtB76kaSTqPe5Zq/keS1MHUBuCTnd9z9W8ArkizsXR30FXT4rVSGXpLm0HEuPPd64M1JdgN76f4byz4J/AfwVWA3sLuqPjPtDF7UTJLa5hG9JDXO0EtS4wy9JDXO0EtS4wy9JDXO0EtS4wy9JDXu/wFdQG4fabRuGQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Numerically stable version of funtion using simple trignometry\n",
    "\n",
    "def f1(x):\n",
    "    return 2*np.sin(x/2)**2/(x*x)\n",
    "\n",
    "#x = np.linspace(-4e-1, 4e-1, 1000)\n",
    "x = np.linspace(-4e-8, 4e-8, 1000)\n",
    "plt.plot(x, f1(x))\n",
    "plt.axvline(1.1e-8, color='red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Another common example of a numerically unstable algorithm. The stable and unstable version of the (unbiased sample) variance:\n",
    "\n",
    "$s^2 = \\frac{1}{n-1} \\sum (x-\\bar{x})^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numpy: 0.08404450403604484\n",
      "Direct: 0.08412869649009781\n",
      "Sum of squares: 0.0\n",
      "Welford's: 0.08357807855169336\n"
     ]
    }
   ],
   "source": [
    "# direct method\n",
    "# squaring occuring after subtraction\n",
    "def direct_var(x):\n",
    "    n = len(x)\n",
    "    xbar = np.mean(x)\n",
    "    return 1.0/(n-1)*np.sum((x - xbar)**2)\n",
    "\n",
    "# sum of squares method (vectorized version)\n",
    "# pay attention to the subtraction of two large numbers\n",
    "def sum_of_squares_var(x):\n",
    "    n = len(x)\n",
    "    return (1.0/(n*(n-1))*(n*np.sum(x**2) - (np.sum(x))**2))\n",
    "\n",
    "# Welford's method\n",
    "# an optimized method\n",
    "def welford_var(x):\n",
    "    s = 0\n",
    "    m = x[0]\n",
    "    for i in range(1, len(x)):\n",
    "        m += (x[i]-m)/i\n",
    "        s += (x[i]-m)**2\n",
    "    return s/(len(x) - 1)\n",
    "\n",
    "\n",
    "# check the performances with an array \n",
    "# of randomly distributed data around 1e12\n",
    "x_ = np.random.uniform(0, 1, int(1e3))\n",
    "x = 1e12 + x_\n",
    "\n",
    "# correct answer from a purpose-built function in numpy\n",
    "print(\"Numpy:\", np.var(x_))\n",
    "print(\"Direct:\", direct_var(x))\n",
    "print(\"Sum of squares:\", sum_of_squares_var(x))\n",
    "print(\"Welford's:\", welford_var(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. The example of the Likelihood: $\\mathcal{L} = \\prod_{i=0}^{N} prob(x, \\mu)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L = 0.0\n",
      "log L = -955.5990036886642\n"
     ]
    }
   ],
   "source": [
    "# loss of precision can be a problem when calculating Likelihoods\n",
    "probs = np.random.random(1000) # Generating 1000 random numbers between 0 and 1, as if they were probabilities\n",
    "#print(probs)\n",
    "print(\"L =\", np.prod(probs))\n",
    "\n",
    "# when multiplying lots of small numbers, work in log space\n",
    "print(\"log L =\", np.sum(np.log(probs)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Summary:\n",
    "\n",
    "- Well-/ill-Conditioned refers to the problem; Stable/Unstable refers to an algorithm or numerical process.\n",
    "- If the problem is well-conditioned then there is a stable way to solve it.\n",
    "- If the problem is ill-conditioned then there is no reliable way to solve it in a stable way.\n",
    "- Mixing roundoff-error with an unstable process is a recipe for disaster.\n",
    "- With exact arithmetic (no roundoff-error), stability is not a concern."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}