02_NumberRepresentation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Integers\n",
    "\n",
    "Integrers number are represented by N bits words. Python 3 allows you to store integers with practically unlimited precision, the only limitation comes from the (contiguos) space available in memory.\n",
    "In python 2, N depends on the PC architercture, N=64 in modern computers.\n",
    "\n",
    "Typically 1 bit is dedicated to specifying the sign fo the number, thus 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 $i$ and the $b=\\alpha_{N-1}\\alpha_{N-2}..\\alpha_0$ is the binary representation of the number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9223372036854775807\n",
      "True\n",
      "9223372036854775808\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "print (sys.maxsize)\n",
    "\n",
    "# check \n",
    "print (2**63-1 == sys.maxsize)\n",
    "\n",
    "# python 3 doesn't have a limit for integers\n",
    "print (sys.maxsize+1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "incidentally note that python 3 recast int to floats when dividing them with `/`. To keep the old functionality of `/`, use instead `//`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6666666666666666\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "print (2/3)\n",
    "print (2//3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Binary and Hexadecimal representations\n",
    "\n",
    "numbers values (in python as all the other languages) are assumed to be expressed as decimal. Built-in functions allows explicitly to convert from one representation to another.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binary representation of 18 : 0b10010\n",
      "Hexadecimal representation of 18 : 0x12\n",
      "Decimal representation of 0b10010 : 18\n",
      "Decimal representation of 0x12 : 18\n"
     ]
    }
   ],
   "source": [
    "# an integer in decimal representation\n",
    "a=18\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",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bitwise operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logical operators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Logical AND  12\n",
      "Logical OR  61\n",
      "Logical XOR  49\n"
     ]
    }
   ],
   "source": [
    "a = 60            # 60 = 0011 1100 \n",
    "b = 13            # 13 = 0000 1101 \n",
    "\n",
    "c = a & b;        # 12 = 0000 1100\n",
    "print (\"Logical AND \", c)\n",
    "\n",
    "c = a | b;        # 61 = 0011 1101 \n",
    "print (\"Logical OR \", c)\n",
    "\n",
    "c = a ^ b;        # 49 = 0011 0001\n",
    "print (\"Logical XOR \", c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unary operators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Negation of a   -61\n"
     ]
    }
   ],
   "source": [
    "c = ~a;           # -61 = 1100 0011\n",
    "print (\"Negation of a  \", c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Masks and shifts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "c = a << 2;       # 240 = 1111 0000\n",
    "print (\"Left shift (towards most significant) of two positions \", c)\n",
    "\n",
    "c = a >> 2;       # 15 = 0000 1111\n",
    "print (\"Right shift (towards least significant) of two positions \", c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Floating point numbers\n",
    "\n",
    "non-integer number cannot be represent with infinite precision on a computer. Single (also known as *float*) and double precision numbers are assigned 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 in 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",
    "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\\times2^{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\\times2^{45}$.\n",
    "\n",
    "An example is given below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "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": 1,
     "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": [
    "### 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": 7,
   "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 could not store them properly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.00000000000001\n",
      "7.000000000000001\n",
      "7.0\n"
     ]
    }
   ],
   "source": [
    "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": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1\n",
      "True\n",
      "True\n",
      "False\n"
     ]
    }
   ],
   "source": [
    "# is 1/10 the same of 0.1?\n",
    "print (1/10)\n",
    "\n",
    "# but then whatch out!!\n",
    "print(round(0.1) + round(0.1) + round(0.1) == round(0.3))\n",
    "print(round(0.1+0.1+0.1) == round(0.3))\n",
    "print((0.1 + 0.1 + 0.1) == 0.3)\n"
   ]
  },
  {
   "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": 18,
   "metadata": {},
   "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": 26,
   "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!\n"
   ]
  },
  {
   "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 a corresponding number of decimal digits are concerned:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 23 bits are used for f in single precision floating points \n",
    "print (2**-23)\n",
    "\n",
    "# 53 bits are used for f in single precision floating points \n",
    "print (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": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 27,
     "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": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-67108863.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": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 3,
     "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": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "-1006.86942381\n"
     ]
    }
   ],
   "source": [
    "# (we'll see numpy soon, bear with me for the moment)\n",
    "\n",
    "# loss of precision can be a problem when calculating likelihoods\n",
    "import numpy as np\n",
    "probs = np.random.random(1000)\n",
    "print (np.prod(probs))\n",
    "\n",
    "# when multiplying lots of small numbers, work in log space\n",
    "print (np.sum(np.log(probs)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## From numbers to Functions: Stability and conditioning\n",
    "\n",
    "Suppose we have a computer algorithm $g(x)$ that represents the mathematical function $f(x)$. $g(x)$ is stable if for some small perturbation $\\epsilon$, $g(x+\\epsilon)\\simeq f(x)$\n",
    "\n",
    "A mathematical function $f(x)$ is well-conditioned if $f(x+\\epsilon)\\simeq f(x)$ for all small perturbations $\\epsilon$.\n",
    "\n",
    "That is, the functionf(x) is **well-conditioned** if the solution varies gradually as problem varies. For a well-conditinoed function, all small perutbations have small effects. However, a poorly-conditioned problem only needs some small perturbations to have large effects. For example, inverting a nearly singluar matrix is a poorly conditioned problem.\n",
    "\n",
    "A numerical algorithm $g(x)$ is numerically-stable if $g(x)\\simeq f(x′)$ for some $x′\\simeq x$. Note that stability is a property that relates the algorithm $g(x)$ to the problem $f(x)$.\n",
    "\n",
    "That is, the algorithm $g(x)$ is **numerically stable** if it gives nearly the right answer to nearly the right question. Numerically unstable algorithms tend to amplify approximation errors due to computer arithmetic over time. If we used an infitinte precision numerical system, stable and unstable alorithms would have the same accuracy. However, as we have seen (e.g. variance calculation), when using floating point numbers, algebrically equivaelent algorithms can give different results.\n",
    "\n",
    "In general, we need both a well-conditinoed problem and nuerical stabilty of the algorihtm to 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 properly redefying the function as in the example above with the likelihood. More examples follow  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t1 = 61249.008531503045\n",
      "t2 = 158057.9134162482\n",
      "% change in x = 0.0006366263894271296\n",
      "% change in tan(x) = 158.05791343536947\n"
     ]
    }
   ],
   "source": [
    "# The tangent function is poorly conditioned\n",
    "\n",
    "import math\n",
    "x1 = 1.57078\n",
    "x2 = 1.57079\n",
    "t1 = math.tan(x1)\n",
    "t2 = math.tan(x2)\n",
    "\n",
    "print ('t1 =', t1)\n",
    "print ('t2 =', t2)\n",
    "print ('% change in x =', 100.0*(x2-x1)/x1)\n",
    "print ('% change in tan(x) =', (100.0*(t2-t1)/t1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEJCAYAAABv6GdPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8HFeV6PHf6W7tiy1bkldZ8r7G\niWNlNXYSEoiz8gghLIFh2EImrMPAsAcmA+8xZAgwMzCQN8DMSyABQkISsjj7BnHi3bG8b7JkWZZs\nyVotqZf7/uiuVlvqlqpbVd2ydL6fDx8sqbp106o+fevcU+eKMQallFJjiyfTA1BKKeU8De5KKTUG\naXBXSqkxSIO7UkqNQRrclVJqDNLgrpRSY5AGd6WUGoM0uCul1BikwV0ppcYgX6Z+cWlpqamqqsrU\nr1dq7NmzJ/z/CxdmdhzKVZs2bTphjCkb7riMBfeqqio2btyYqV+v1Nhz+eXh/3/ppUyOQrlMRGrt\nHKdpGaWUGoM0uCul1BikwV0ppcYgDe5KKTUGaXBXSqkxSIO7UkqNQRrclVJqDNLgrlzz0p4mjpzs\nzvQwxoUDzZ20nfZnehhqFNHgrlwRChluv38Tv/7roUwPZVz4xcsH2NfUmelhqFFEg7tyxYnOXnr8\nIfoCoUwPZVzoDYQIBEMEQ7rhvQrT4K5cUdcaTscEghps0sF6nXv1w1RFaHBXrqhrOQ1AQGeSaREI\nhYN6TyCY4ZGo0UKDu3JFXUt45h4M6UwyHax0TK9fX28VpsFducJKy/h15p4W/mhaRmfuKkyDu3KF\nlZYJas49LXTmrgbS4K5cEV1Q1bRMWviD4ddZF1SVRYO7clwgGOJYW0/435qWSYvozD0QxBh9zZUG\nd+WCY2090WCjddfpEYh5vU91652qSoO7coFVKeOR/nSBclds+stKianxTYO7cpwVXGaU5OnMPU0C\nQUNOlhfoX8xW45sGd+W4+tbTeD1CRUm+5tzTJBgy5PrCb2eduSvQ4K5cUNfSzbQJueT4PNp+IE0C\nIYPP68Hn9UTTYmp80+CuHFfXepqKkny8Ho/O3NMkEAohQI7PQ12rpmWUBnflgrqWbiom5ZHlFQK6\noJoWgaBBJBzc63XmrtDgrhzW4w/S1NHLzJJ8vB7RBdU0CYQMIkJOlpf6U6cJ6es+7mlwV46qj6QE\nwjN3TcukSzBkommZvkCI5s7eTA9JZZgGd+Uoq1KjIjJz17RMeviDoUhaxiqH1NTMeKfBXTnKyvdW\nTMrH5xGduadJeOYu5GZpOaQK0+CuHFXXeppsn4eywhx8Xg3u6RLOuYfTMqA3MikN7sphdS3dzCzJ\nw+MRfB6PpmXSJBAMl0J6RCgvytG0jNLgrpxV19pNRUk+AD6tlkmLUMgQMiAiQDglpmkZpcFdOaqu\n5TQVk/IA8HpFd2JKAyv1FYntVJTkaVpGaXBXzmnv8dN22q8z9zQb+BpXTMqnsb1HU2LjnK3gLiJr\nRWSPiOwXka/G+fksEXlRRLaIyHYRudb5oarRrj4yW5wZDe4egiGjm0e4zGr364lM3WdGunFaG6ao\n8WnY4C4iXuCnwDXAEuADIrJkwGHfBH5vjFkBvB/4mdMDVaPfkWgZZDgt4/OEg41WzLhrYHO2iknh\nD9cjuqg6rtmZuV8I7DfGHDTG9AEPAu8acIwBiiP/ngA0ODdEdbY40NwJwJyyQiCccwfdjcltA3Pu\ncyOvv/X3UOOTneA+A6iL+bo+8r1Y3wE+JCL1wJPAZ+M9kYjcJiIbRWRjc3NzCsNVo9mB5k6mFudS\nmOMDIMsTPr10NyZ3BaPBPRzdy4tyKMzxcaBJg/t4Zie4S5zvDZyKfQD4b2PMTOBa4D4RGfTcxph7\njTHVxpjqsrKy5EerRrUDTZ3MLS+Ifu316Mw9HawPT+uNKiLMLS9kv87cxzU7wb0eqIj5eiaD0y4f\nB34PYIx5HcgFSp0YoDo7GGM40NzFvEhKACDLqzn3dAgOSMsAzC0r4EBTV4ZGpEYDO8F9AzBfRGaL\nSDbhBdPHBhxzBLgSQEQWEw7umncZR46399LZG2BueX9w90bSMrobk7usahmJucieV15IY3sPHT3+\nTA1LZdiwwd0YEwA+A6wDdhGuiqkRkbtE5MbIYf8AfFJEtgEPAH9rtP5tXLEW72Jn7v3VMppzd9PA\nBVXoX1Q92Kyz9/HKZ+cgY8yThBdKY793Z8y/dwKrnB2aOpvsjyzexc7cfVotkxbWlVHs4ti8yN9h\nf1Mn51ZMzMCoVKbpHarKEQeaOynK8VFelBP9nrWg6te0jKvizdxnRVouaznk+KXBXTlif1Mnc8oL\no+V4EL5DFXTm7raglXOPee2zvB4qJ+dHr6jU+KPBXTniQHPnGfl26E/LaJ27u+KlZSCcmtGZ+/il\nwV2NWEePn+PtvWfUuEP/gqrO3N0VLy0D4UXV2pPd+uE6TmlwVyN2IFKRMXfQzD1SCqnB3VXR4D5g\n7j6vvJBAyFB7UnvMjEca3NWIWbe5zysfENytUkidObrKen3jzdxBe8yMVxrc1Yjtb+7E5xFmRboR\nWrT9QHr0z9zPNDemHFKNPxrc1YgdaOqkqrSALO+Zp5PVfkB3Y3LXwMZhlsIcH1OLc3XmPk5pcFcj\ntr+5k7llBYO+742WQmpaxk3+BGkZiFTM6Mx9XNLgrkbEHwxx5GT3oHw7xObcdebupmCCtAxEGog1\nd+luWOOQBnc1IrUnuwiEzKBKGeivc9dqGXcFEqRlIDxz7+wNcLy9N93DUhmmwV2NyP6m+GWQoNvs\npUuim5hAK2bGMw3uakSsoDE3blpGc+7p0N9+YPDP5mnFzLilwV2NyIGmM7fWi6WNw9LDen3jpWXK\ninIoyvFpcB+HNLirEdl5rJ0FU4vi/kxb/qbHUK+viLBgahG7G9vTOCI1GmhwVynrDQTZ39TJ0unF\ncX/ui+7EpGkZN1lrGp54SXdg6fRidja0E9IP2XFFg7tK2d7GTgIhM0Rw1wXVdOj/8Iwf3ZdOL6ar\nL0hti/aYGU80uKuU1TS0AbB0+oS4P9e0THok6gppsf4+1t9LjQ8a3FXKahraKczxUTmgp4zFSsvo\ngqq7giGD1zOwJ2S/+VMK8XmEmgbNu48nGtxVymoa2lgyrRhPgmRvf+Mwzbm7yR8KRV/reHJ8XhZM\nKdLgPs5ocFcpCYYMu451sCRBvh36c+46c3dXMGjIGiK4g7Wo2qZtCMYRDe4qJYdOdHHaH0y4mArg\n8Qge0Zy72wKRtMxQlk4v5kRnH00d2oZgvNDgrlIy3GKqxef1aLWMywKh0KB2ywMtnaGLquONBneV\nkp0N7WR7PcyfMrjtQCyfR7TO3WVBGzP3xdOKEYGao5p3Hy80uKuU7GhoY8HUwmFnjF6P6MzdZf6g\nia5vJFKY46NqcgE7dOY+bmhwV0kzxlDT0M7SaUOnZACyvB7NubssGDLRzciHsmR6sVbMjCMa3FXS\nGtp6ONXtZ+mMxIuplvDMXdMybgqEhp+5Q3hRtb71NG3d/jSMSmWaBneVtJqj9hZTwcq568zdTYHg\n0HXulmXWnarHNDUzHmhwV0mraWhHBBZPi98NMpbPqzl3twVspmWsstWdmpoZFzS4q6TVNLQzp7SA\n/OzBPdwH8nm0FNJtgWDIVlpmcmEOU4tzNe8+TmhwV0nb2dBmKyUD4bSMth9wV3jmPnxwh/DsXWvd\nxwcN7iopJzt7aWjrGfLO1Fhej2j7AZcFbS6oQji472/qpLsv4PKoVKZpcB+Dmjp6aOroceW5N9W2\nArBiVomt431e0VJIlwWCw9/EZFkxq4SQga11p1wZy7G205zs1BYHo4Gt4C4ia0Vkj4jsF5GvJjjm\nFhHZKSI1IvJbZ4epkvGJ/9nIl/+w3ZXn3lTbSpZXWD7TblpGc+5us9N+wHJ+5EN50+FWV8by2d9u\n4bMPbHHluVVyhl0RExEv8FPgHUA9sEFEHjPG7Iw5Zj7wNWCVMaZVRMrdGrAa2pGT3Wyvb2PxNHtp\nk2RtrG1l2YwJ5GZ5bR2v7QfcZ6dxmGVCfhYLphSysdad4H68o4f61tM0d/RSVpTjyu9Q9tj5uL8Q\n2G+MOWiM6QMeBN414JhPAj81xrQCGGOanB2msuupHccA6Ox1/kaVHn+Qt+rbqK60l5IBbT+QDgEb\n7QdiraycxOYjra7sqdrZE8AYeGZno+PPrZJjJ7jPAOpivq6PfC/WAmCBiPxFRNaLyFqnBqiS89SO\n8Juqo8f5BbOahjb6giFWVk6y/Zgsr0dn7i4LL6jaXz6rriyhoyfAvqZOR8dhjImed0/v0OCeaXbO\niHhTgoEf+T5gPnA58AHgv0Rk4qAnErlNRDaKyMbm5uZkx6qGcaztNFvrTpGf7aWjJ+D4xgwbI3na\nlUnO3HVB1V3+UAivzVJIgOqq8N9vY22Lo+Po8YcIhAz52V5eP3CSU919jj6/So6d4F4PVMR8PRNo\niHPMo8YYvzHmELCHcLA/gzHmXmNMtTGmuqysLNUxqwSs2dL1y6cRDBl6/M7OmDfWtlI1OT+pXGqW\n3qHqumBo+J2YYs2alE9pYY7ji6odkVTg9cunEQgZnt153NHnV8mxE9w3APNFZLaIZAPvBx4bcMyf\ngCsARKSUcJrmoJMDHavae/ys/sELPPnWsRE/11M7GlkwpZDlM8MXTR09zuXdjTFsrm1NKiUDkZy7\n1rm7KlwKaT8tIyJUV5Y4vqhqpWQunVvKjIl5jqRm/rCxjrf/8CV6/MERP9d4M+wZYYwJAJ8B1gG7\ngN8bY2pE5C4RuTFy2DrgpIjsBF4EvmyMOenWoMeShzbWU9dymm31I6s7bu7oZcPhFtYum0ZRbrgI\nqqPXubz7oRNdnOzqi17S2xUuhdScu5sCIXvtB2JVV5VwpKXb0fshOiPBvSjXx9VLp/LqvhMjnmBs\nrTvFweYuHts2MFmghmPr494Y86QxZoExZq4x5nuR791pjHks8m9jjPmiMWaJMeYcY8yDbg56rAiF\nDPetrwWguX1kN348s7MRY+CaZVP7g7uDi6rWLC+ZShnQxmHpEEyi/YDFWjdxMjXTEQ3uWVxzzlT6\ngiFe2D2ywjlrz9f/+eth3dw7SXqH6ggEgiE++us3eXTr0ZQe/+r+Exw60YXXIzSP8K6+p3c0UjU5\nn0VTiyjMyQKcTctsOtzKhLws5pYNva3eQJqWcZ+dnZgGWjp9Ajk+j6OpGet8K8zxsXJWCWVFOSNO\nzTR39OL1CDUN7Ww+ktrV7YNvHuH2+zaNuw8HDe4j8OjWBl7c08zzu1Kbnfy/vx6mtDCHNfNLaRrB\nzP1Udx+vHzjJ2mXTEJHozL3T0Zl7CysrS/AkGUSyPLoTk9vs7sQUK9vn4dyZE50N7r39aRmPR7h6\n6RRe2tM8oj42zR29vGPxFIpyfPy/1w+n9BzP7jzO0zWNrKsZXwu8GtxTFAiG+PcX9gFQ29Kd9OOP\nnOzmhT1NfPDCCqZPzBvRzP3x7ccIhAzXL58G4HhaprWrjwPNXUmVQFq8Xt2JyW2p5NwBVlaVUHO0\njdN9zixWWudbcW74yvGG5dM57Q/yTIpB1RhDc0cvlaX53Fw9kyffOpbSGoH1/vzxc3tduXFrtNLg\nnqJHthzl8MluZk3K58jJrqQff/8btXhE+OBFlZQV5dDS1Yc/xZt9Htlcz8IpRdFOjUVWWsahBdVN\nKebbIdJ+YBy9oTIhmcZhsaorSwiEzIgX8y3WlWJBTrg1xQVVk5gxMY+Ht6SWtmw/HaAvGKKsMIcP\nX1yJP2h48M264R8YIxQyHGkJv093N3aMqztnNbinwB8M8e8v7GfZjGI+cOEsWrv9tCeR3z7dF+R3\nG+pYu3QqUyfkUl6UC8CJFGbvh050sfnIKW46fwYi4Td4YXTm7kzO/Y1DJ8n2eji3YtB9acPyeTwE\nNefuGmOM7Z2YBlpZWYIIrD/oTGFbR4+f/GxvdCwej3DT+TN4bV8zx9uTn3Fbs/Ty4lzmlBWyZkEZ\nv3mjNqlJ0PGOHvoCIT6xejZzSgv48XP7xs3sfVwH96OnUtss+JEtRznS0s0XrlxA1eR8IJxmsetP\nW4/SdtrP31xSCRC9Kai5I/ng/sjmejwC/2tFf0cIr0coiNyl6oRX953ggtkltpuFxfJ5Bb+mZVxj\nxalU0jIT87NZPmMCr+074chYOnoCFOac2Yvw3StmEDKkVHRgvR/KCsPvj49cUsnx9t5oiw07aiPv\ny9mlBXzuyvnsbuzg6ZrkZ+8tXX00trnTRtstYya4J7to19LVx3X/9irfebwmqcf5I7n2c2ZM4MrF\n5cyKBPdam8E9GDLc+8pBls0o5sLZ4RuCyiPBPdlF1VDI8PCWo6yaV8qU4twzflaY63NkQbWpvYfd\njR2snp/aHcXafsBd1iw2lbQMwOr5ZWypO5XUlWcinb2B6HqPZU5ZIStmTeThzckHd6sMsrw4/P64\nYmE5c0oL+MXLB2xXvliTrspJBdxw7nTmlhXwkxRm71/543Zu/I/Xkr4azuS5PyaC+4NvHuG8f3om\nqU0CfvD0bk51+5PON/4hctPRF66aj4hQObkAgNoWe3n3Z2oaOXSii9svmxtNo0Rn7kmmZTbWtlLf\nepqbzh/Yxy1ca9zhQGfIVyOzutXzS1N6fFZkJ6bxVoaWLlbwyEqyzt2yen4pwZDh9QMjT8209/gp\niiymxrppxQx2N3YkvTF3dOYeeX94PMJta+ZQ09DOa/vtXW3UtnTh8wjTJ+bi9Qifu3I+e4538Pj2\n5G6K2lZ3iqaOXn7y3D7bj2k4dZpzvrMuYzdgnfXBva6lm7v+vJOO3gBvHrLXCGnLkVZ+t7GOkvws\nDp3osl2q1dUb4EfP7aW6soS3Lwq3rC/M8TG5INtWWsYYw89fPkDl5HyuWTYt+v3SwtTSMg9vric/\n28vVS6cO+llhjs+RtMyr+5opLcxm8dTU+sNbt8Xr5N0d1j0EybQfiLViVgkF2V5e3TfyRn7xZu4A\n1y+fTpZXeHhzfVLP19zZS47PQ1FMqufd58+gvCiHn798wNZz1J7sZkZJXnQd4Ibl01k8rZh/fWYP\nvQF7VUInOntp6uilJD+LX//1MHsaO2w97i/7T9DdF+Rbf9pBUwprDiM16oL7fa8ftr0FWChk+PJD\n2/CIkOPzsMHG3XbBkOHOR2soL8rhW9cvwRhs/7H+69VDNHf08rVrF0dn3QCzJufbSsu8fvAk2+rb\nuG3NnDMuo7N9Hkrys5Iq8+rxB3li+zGuWTaN/OzBb6ii3JEH91DI8Nr+E7xtXmnS9e0W685JLYd0\nh/W6pjpzz/Z5uGTu5OgV2kh09MQP7iUF2bx9UTl/2tqQVPvnpvYeyotzzniv5fi8fPxts/nL/pNs\nt3HVbVXKWDwe4evXLqKu5TT3rz9iaxy7joWvOL77v86hKNfHnY/usHUluvFwK4U5Pnr8Qb7+yFu2\nr17fPNTC7zbYG9tQRlVwf3hzPd96tIZ7nt1r6/jfvFHL+oMtfPO6xZxbMdFWC9MH3jzCW0fb+MZ1\nS7igKpzz3nls+MvF5o5efvHKAa5ZNnVQvXflpHyO2Kh1/8+XDlBamMN7zp856GdlRTlJzdyf3Xmc\njt4A74mTkoFwrfFIq2V2NbZzorMv5Xw79C/06V2q7rDSMqnm3CGcd6892U1tCiW9sTp6/IMWVC03\nnT+TE529SX2INHf2RhdTY33wolkU5fpszd5rT3ZTOTn/jO+tnl/G6vml/PsL+2g7Pfx7xEonXTp3\nMv949SLeONRiK9WyobaFi2ZP4stXL+S5XU38yeai8vef2sVX/vgW61JY+I01aoL7geZOvvmnHXg9\nwvqDJ4e9saKupZv/89RuVs8v5X0XVHBBVQk1De1DplhOdvZy97o9XDJnMjcsn8bMkjyKcn3RT+ah\n/OT5vfQFQnz56oWDfjZrcgENbaeHvMzbcbSNV/ed4GNvq4pbdVJelBtdQLLjwQ1HmDExj4vnTI77\n88IcH50jrHMfab4d+oOO1rq7wx95XVOplrFYf9+Rzt47ewJxc+4QXgydVJDNg0nMSJvae6NlwrGK\ncrP48MWVPLUjvH6VSFu3n7bTfionFQz62VevWUTbaT//+dLwHxC7jrUzbUIuJQXZvO+CCpbPnMD3\nntg15CL0yc5eDjZ3UV01iY+ums3KyhK+89jOYdMzp7r72Fp3Cq9H+PIftlHfmvwNkpZREdx7/EE+\n/ZvN5Pg8fP+mc+gLhHj9YOITLRiTjvn+e5aHW5hWTSIYMmwdov/EXX/eSXdfgH9611JEBBFh8bTi\nYRd6DjR38sCbdXzwolnMidNbpXJSPsZAfevphM/x85cPUJjj49aLKuP+PJmZ+/6mTv6y/yQfvGhW\nwnSJE2mZV/c1s2hqEeXFg99gdlkbN2vFjDusewiS2YlpoNmlBcyYmDeivHswZOjqCyacuWf7PNxS\nXcGzO49zrC3x+yRWc2fifVg/umo2WV4P976SODhbRQ6zBszcIdxb590rZvCrvxzi6Kmhx7PzWDtL\nInsSez3CP79rGSc6e/mXp3YnfIx1498FVSV4PcLdNy+nxx/kaw8PnZ55Zd8JQgbuueVcQgY+98CW\nlG9uHBXB/btP7GR3Ywf33HIeN543nbwsLy/tSXyi3fvKQdYfbOHOG5YwY2IeEN7VXYSEefcXdh/n\n0a0N3HH5PBZMKYp+f8m0YnY3dgxZGvV/ntxNXpaXz105aP8RgOhlX6JF1T2NHTzx1jE+dHElE/Li\nz2ys4G4nL/ebN2rJ8gq3VFckPKYw10d3XzDloHq6L8iGQ60jmrVDzMxdt9pzhZVzT7YrZCwRYc2C\nUv66/2TKf6fOmL4yidx60SwM8MAbw8/eewNBTnX7Ewb3sqIc3lddEalei/++s9bBBqZlLP/wzvBV\n+L+u25NwHD3+IAeau87YcP7ciol8bNVsfvPGEd5IcAPYxtpWsn0ezpk5AQiXhH71mkU8v7uJ+4f4\n739pdxMl+Vlcv3w633/POWw+coofPmMvTT1QxoP749sauH/9EW5bM4crFpWT4/Ny6dzJvLinKW6g\ne6u+jR8+s4frzpnGe1f2564n5GWxcEpR3Lx7R4+fbzyygwVTCrnjirln/GzJtGK6+4IJ+8O8uKeJ\n53Yd544r5karWgbqr3WPf4n44+f2UpDt41Nr5sR/EQjXuvcGQrQPM9vu7gvw0KZ61i6bNuSOSNbl\ncaq17m8cOklfMDSifDv0L/RpWsYdgdDIZ+4QzkN39AZSbkVgre8UJ0jLAFRMyufyBWU8sKFu2Nno\nic7wFn3lQ5zjn75iHh6P8G/Pxy9PtNbBYhdUY82YmMdtq+fwyJajbDgcf71u3/FOgiHDkulnVot9\n8Z0LqJiUx1cffivuRiIbDrewfMYEcnz9Kdi/vbSKyxaU8d0/72Tf8cFFHKGQ4eW9zaxZUIbXI1y/\nfDofuHAWP3/5AC/sTr4/T0aD+57GDv7xoe1UV5bwpXf257IvX1ROXctpDg7Ip3X3Bfj8g1soK8rh\nf7/7nDNW0SG8AcGWI6cGzVZ/8PQeGtt7+P57lp/xYgPRP1q8vHtvIMhdj+9kTmkBn3hb4sBcVphD\nfrY37gdETUMbT+1o5GOrqigpyE78HDbvUn18WwMdPQE+fHH89I7FmkGlenPKq/tOkO3zRG+0SpVV\noqcLqu7oL4VMfeYO4cVCj8Are1PLu1spwMIhZu4AH76kkuaO3mGbiQ2scY9n6oRcbr1oFg9vORo3\n9157souyopy41WSWO66Yy/QJudz5aE3cq1wrLsTO3AHys318/6blHDrRxY8H1L73+IPsONpGddWZ\n7x0R4e73Lqcwx8fnHtw6aI1uR0MbJ7v6uGJhefR7375hCUunF/OFB7dyeIj1hXgyFtyDIcOn7ttI\nYa6Pn916Ptm+/qFcviA8WxyYmvnnP+/k0Mku7rnlPCbkD54hXFA1ic7eALsb+wP1m4dauG99LR+9\ndDbnzxrc+GpeeSFej8TNu//ytUMcOtHFd25cesb4BhIRZk3Kj3t5+KNn91Kc6+PjqxN/OED/STxU\nOaQx4c09Fkwp5IJhdkSyaoNTXVR9eW8zF1ZNSqnlQKxotYyWQroimpYZYXCfmJ/N8pkTeXlvanl3\nO2kZgMsWlDNjYh73RzapScRaeIy3oBrr7y6fS7bXw0+eG5y6qD3ZTWWCWbslP9vHN69fwq5j7fz2\njcFj2nmsnfxsb9znWTWvlPdVV/B/Xz3IW/Vt0e9vqzuFP2jivkfLi3L5wc3L2XWsnbufPjMd9OLu\nZkRgzYL+q+XcLC8//9BKPB7h9vs3JdU+OWPBva61m/rW0/znrecPWrCrmJTPvPJCXtrT3yf9z9sb\neODNOm6/bC6XzI1fIWKVKG6M5N27egN86Q/bmFmSx5euXhD3MblZXuaVFQ6auR9rO82/P7+fq5dO\nOePFTmTWpMG17lvrTvHcriZuWzMnYa7dUm5j5r6tvo0dR9v58MWVg65aBrLSMqksqh5o7mR/UydX\nLS4f/uBhWLlgXVB1RzQtM4Kcu+WqxeVsrTuVUg+V2I06huL1CLdePIvXD55kf1Pi+0usu7WH24y9\nvCiXv7m0kke3NbB3QKrjSEt33MXUga5ZNpVV8yZz97o9g+5y33msnUVTixIWLnz9usVMLsjmS3/Y\nFk3PWD3yE7XIvnLxFD58cSX/9dqhM9ItL+1t4tyZE5k04Aq/YlI+//b+Few53sFX/vjWsP89lowF\n946eAN+6fsmgSxfL5QvKeONgC919AQ40d/KVh7azsrKEL74jfpCGcA5t2oTcaP7su0/spK61m3tu\nOW/IS7PF04oG1bp/74ldhIzhm9ctsfXfUzk5XOseuzB7z7N7KcnP4m9XzR728WWF4Q+4oYL7/etr\nyc/2ntEkLBHr8rgzhRYE1u45Vy8bfOdrsqwZpV/TMq4IOpRzB1gbuWs6lba4sVvsDeeW6gqyvDLk\nTURN7b2IwOTCxKlMy6fWzCU/y8uPY2bvPf4gje09ccsgBxIRvnPDUrr7gvzrM/2zaWMMuxraB+Xb\nY03Iy+IHNy9nz/EOfhh57MbDLcwvL2RifuKxf+O6xSyZVszf/y5c7tjSFS6BvHxh/InkmgVlfOmd\nC3k8iVYGGQvuJfnZ0a6I8Vw5nXyXAAAZxklEQVS+sJy+YIjndzVxx/2bycny8h8fXBEtrYvHKonc\neLiV53Ye54E36/jUmrnD5o2XTC/mWFsPrV3hRZyX9zbz5+3HuOPyeVQMc1lnmTW5gN5AKFqrvv7g\nSV7Z28ztl80ddjYDUJznI9vnSRjcmzp6eGxbA+9eMcPWG2gkG3asq2nkvIqJTJuQl/RjB7Jy7jpz\nd8dIG4fFmldeyLzywpS2xuvfqGP4c720MIfrzpnGQ5vqE95E1NzZy6T87CHf75ZJBdl87G2zefKt\nxuhdq/Wt3RiTuFJmoPlTivjoqioe3FDHxsjksL71NB29gUH59oEuX1genYn/Zf8JNta2Jpy0WnKz\nvPzs1vMJhQyf/u0Wnt91HGPCz5XI3102l2uSmHBlLLjPLMkbMrVwwewS8rO9fOWP29nb1MGP33ee\nrWBTXVlCY3sPX/z9VhZPK+bv3xG/fDGW9cfbdaydrt4AX3/4LeaWFXD75UPnyWNZObnak12EQobv\nPrGT6RNy+cilVbYeLyKUFSaudf+fvx7GHwzxiWFy95b+BdXkgnt9azfb69tY68CsHfrTBdr21x0j\nbRw20NqlU3njUAstkYmOXXYXVC2fWD2Hzt4Av01QFtjckbjGPZ5PrpnD5IJsvvvELowx0RSpnbSM\n5QtXLWD6hHAFTG8gGL2aXzJMcAf42rWLmD25gL+7fxMdPQFbG9tUlRZw93uXs63uFN9+rIbJBeEW\nzIl4PMLPbj3f9n9PxkshEwmXRJbS3Rfks2+fbyvvDeGKGYCeQIgfv++8QdUx8VjBfeexdv71mT00\ntJ3mX+JU1gzFmiHUtnTzyJaj7Djazj+uXZTUgmR5cU7cu1Q7ewPc93ota5dOZXbp8JeZ0L8bU7Kl\nkNY+k2vjNCNLhZWW0Zm7OwIOtB+ItXbZVIIhw3M7kyu96+z14/UIeTbP92UzJrB6fim/+suhuHd2\nNyUZ3Itzs/j7dyzgzUMtrKs53l/jbvPKG6Agx8f/vukc9jd18tMX9rOzoR0RWDi1aNjH5mf7uOd9\n59EVubP+gmFm7pa1y6bxsVWz6e4LsmZB2bA9nIZba4s1aoM7wO2XzeFTl83h8wluHopn0dRiVsya\nyLdvWGLrjwLhy8Tyohwe3nyU//7rYT58ceWwl1UDTZ+Yh9cj7Gns4O51ezh35gRuPHd6Us+RaOb+\n4JtHaO8JcNsQdfID5WZ58Hkk6f4y63Y0smhqEVU2P0SG49NSSFdZr6ud9IUdS6cXM7MkL+kNLayN\nOpIJPp9aM5fmjl7+FGcbvhNJBneA919QwfzyQr7/1C72N3dSmOMbtDg5nMsWlHHTihn87KUDrKtp\nZHZpwZDrdbHOq5jI165ZxCVzJlMxyX5K86vXLOJjq2bzMRtrc8kY1cG9umoSX7tmcVKzEq9HeOSO\nVQlv809kyfRidh5rZ2pxbtz+McPJ8nqYMTGP+16vpbG9h29evyTpTorhmfuZlQr+YIhfvnaIi2ZP\nYkWcUs5ERCS8YUcSpZBNHT1sqG1xLCUD2hXSbcGQczl3CJ83a5dO5bV9J5KaGHQm6Ag5lFXzJrN0\nejG/eOXgGYUI1sbYw5VBDuTzevj6dYs5fLKbhzbWM2tSflIfNpZvXr+E4rwsdjd2DJtvH+gTq+fw\nwG0XJ/V7s30e7rxhSfRuVqeM6uCeTtbm0t979zJbC5bxVE7Opy8Y4tpzptq+LItVVphLa7efvkB/\nIHx8WwPH2nq4/fK5QzwyvmT7yzy7M7yo42Rw18Zh7vJHe8s4E9wh/PfvC4Z4YXfT8AdHtA/RNCwR\nEeH2y+ZysLmL53b1p4HaTvvDG2MnOXOHcJXd6vml9AVDthdTB5pUkM23bwhXyS0dolJmtNPgHvHR\nVbP55UeqefuiKSk/x9yyQrK8wlfWLkrp8dbJfLIrnJoxxvCLlw+ycEpR9MauZBTlJNf29+kdjVRN\nzmfhFHvpLDuyrGoZTcu4IloK6VBaBsJ9msqKcpJqOdvR4z9jUw27rlk2lYpJefw8Zus8O3enJiIi\nfPO6JXgE5pSlnlq88dzp/PxDK7n1wuQyAKOJBveI0sIcrlycemAH+NyV83nkjlXRrfeSNXAv1ed3\nNbHneAe3rZmT0uVlYRIz97ZuP68fOMnaZdNS+l2JePUOVVcFHGj5O5DHI1y9dAov7m4etvW2JdEu\nTMPxeT18cvUcNh85xfqD4RLE6N6pKQR3CC+A/unTq7htdfJXuxYRYe2yqXHvhD9baHB30KSCbJYN\nUco0nNj+MqGQ4Z5n91I5OZ8bz0tuYdZSnERwf+KtYwRChmvPcS4lA7E5d525uyHgYJ17rGuXTeO0\nP8izu+xVzXT0BGyXQQ50S3UFU4pz+NGze6P5dkht5m5ZPnPiWR2YnaDBfRSxdnlv6ujl6ZpGdh5r\n5wtXzU+5EqIwx2d7k+yHNtUxv7yQc0bw4RSP7sTkLifbD8S6eM5kpk/I5aFN9vY97ejxpzRzh/AN\nPZ+5Yh5vHm7h1X0nokUFqc7cVZgG91FkckH4ZG5s7+GeZ/cyr7yQG88dvtVAIkW5Wbbq3A80d7L5\nyCluXjnT0ZQMxJRC6szdFU62H4jl8QjvWTmT1/Y1D9trxhhDZ2+AwpzUZ8q3XFDBjIl5/PCZPTS1\n95Kb5bF1Z7dKTIP7KGJtlP27DUfY39TJ31+1YESX21a1zHAbgDy8uR6PwLtt9KxJVn/jMM25u8Fq\nP+D0zB3gPefPJGTgkTh16LF6AyH8QZPyzB3CNy1+/sr5bKtv4/HtDZQV5Tg+0RhvNLiPMuVFuRxv\n72XxtOKk+kjEU5jrIxAy9AYSB9ZgyPDw5qNctqBsRNvpJaKNw9wVdGFB1VJVWsAFVSU8tKluyAlC\nMn1lhnLT+TOompzP8QR7p6rkaHAfZaxFpC++Y0HSN0ENZNUdD7Vhx18PnOBYWw83r0y8Zd9IeLX9\ngKucbj8w0M0rZ3KguYutdYl3aIq2+x1hcPd5PXzhqnDX17IEu54p+2wFdxFZKyJ7RGS/iHx1iONu\nFhEjItXODXF8uXTeZK5aPMWRXurFNjpDPrSpngl5WVzpwO+Lx6q/TnWTXzW0aPsBh3PulmvPmUZu\nlmfIhdVou98R5NwtN5w7ndXzSxPu2aDsG/ajVkS8wE+BdwD1wAYRecwYs3PAcUXA54A33BjoeHHH\n5fMcey5rQSrRomp7j5+ndzRyS3XFiHdcSkQbh7krGAohwoiv8hIpys3immXTeGxbA9+6fknc88Rq\ncTHSmTuEr0Du+/hFI34eZW/mfiGw3xhz0BjTBzwIvCvOcf8M/ABIfhsX5YrhdmN6YvsxegMhbo7Z\naNxpWufuLn/IuDZrt9y8ciYdPQGeSdAp0krLjGRBVTnPzlkxA6iL+bo+8r0oEVkBVBhj/uzg2NQI\nRWfucWrdjTHcv76WhVOKWO5ww6JY2hXSXcGQcS3fbrlkzmRmluTxmwT7nvYvqI7vm4ZGGzvBPd6Z\nE32niogH+BHwD8M+kchtIrJRRDY2N6e2Ea+yb6gNO9481EJNQzsfubTK1ZIzK+5oKaQ7AkHjSqVM\nLI9H+JtLKnnjUAs1DW2Dfh7dqEPr0kcVO8G9HogtpZgJxG7kVwQsA14SkcPAxcBj8RZVjTH3GmOq\njTHVZWXJN8JSySkeIi3z678cZmJ+liu17bFEhCyvaFrGJYFQyJUa94HeVz2L/Gwvv/7L4UE/S3YX\nJpUedoL7BmC+iMwWkWzg/cBj1g+NMW3GmFJjTJUxpgpYD9xojNnoyoiVbQU54cWvgQuqdS3dPLOz\nkQ9cOIu8bHcWUmN5PRrc3RIImeg+tW6akJ/Fe86fyWNbGzjReeaGMp29fnKzPI5tGKKcMexfwxgT\nAD4DrAN2Ab83xtSIyF0icqPbA1Sp83k95Gd7B7X9vW99LSLChy9OTztTn8ejOXeXBIIh19Mylr9d\nVUVfMDRo39OOFHq5K/fZuo4yxjwJPDnge3cmOPbykQ9LOaUw58zdmLp6Azz45hHWLpvK9In2twIb\nCZ9XtOWvSwIhk5a0DIT3K7hsQRn3ra/l9svmku0Lzw07Umz3q9yl11Fj3MDdmB7eXE97T4CPrapK\n2xh8mpZxTTDk/oJqrI+uqqK5o5cn3upfduvoCaS0UYdylwb3Ma4oNyvafiAUMvz6r4dZPnMC5yex\nH+tI+Twe3YnJJYGgcXQXpuGsmV/GnLICfvXa4Wi/mXC7X03LjDYa3Me4ophNsp/a0cjB5i4+sTq1\nnZ1S5fUIfk3LuCIQSl/OHcJlkZ9cPYe3jrbxyr4TQHjBXssgRx8N7mOclZYJhQw/eT7cI/66c6al\ndQw+r2j7AZek4yamgd5z/kxmTMzjx8+Fd04KL6hqcB9tNLiPcYU5Pjp7Ajy1o5G9xzv53JXz0x4M\nNOfuHn+a0zIQ3nfg01fMY8uRU7yy70Rk/1RNy4w2GtzHOCvnnqlZO1ilkJqWcUO6F1QtN68Mz97v\neXZveBcmnbmPOhrcx7iiXB/dfUH2Hu/ks2+fl/ZZO2haxk3+NNa5x8r2ebjjirlsi/R5H+lGHcp5\nGtzHOGuha25ZAdcvn56RMfg8ojsxuSSYxjr3gd67soLpE8I7JumC6uijwX2Mm5ifDZCRXLvF69GZ\nu1vS1X4gnmyfh0+/Pbz/wIQ8zbmPNvpxO8a9c+kU7jbLMzZrh3AbBL1D1R2BUIisDH1oA9xSXYHP\nI1y2UBsBjjYa3Me44tws3lvtzv6odvk8Qt8Qm3Sr1AWC6S+FjJXl9fC+C2Zl7PerxDQto1ynXSHd\nk87eMursosFduS5L0zKuCZdC6ttYDaZnhXKd1yPa8tcl6W4/oM4eGtyV67K0zt014cZhGtzVYBrc\nleu8Ho/m3F2SyVJINbrpWaFcF+4tozl3N2Sq/YAa/TS4K9f5PKL93F3iD6Zng2x19tHgrlzn8wp+\nTcu4QmfuKhEN7sp12n7APeGbmPRtrAbTs0K5zufx4NeWv64IhEJkaVpGxaHBXbnOpzN3V4RChpAh\no+0H1OilwV25Ltw4TIO706zXNCvNOzGps4OeFcp1Po/oTkwusK6GdOau4tHgrlzn9QghE04jKOdY\n9w5otYyKR4O7cp214KepGWdZ/Xo0uKt4NLgr11mlerqo6izrw9KrOXcVh54VynX9M3fNuzvJej0z\nuROTGr00uCvXWQt+2vbXWdbrqQuqKh4N7sp1Vk5Yc+7OstJc2ltGxaPBXbnO59Wcuxv6q2X0bawG\n07NCuc5KG2gLAmdZV0JaLaPi0eCuXGcFH525O0tz7mooGtyV66y0jFbLOEvbD6ih2DorRGStiOwR\nkf0i8tU4P/+iiOwUke0i8ryIVDo/VHW20gVVdwQjH5Y6c1fxDBvcRcQL/BS4BlgCfEBElgw4bAtQ\nbYxZDjwE/MDpgaqzl09LIV3hD2q1jErMzsz9QmC/MeagMaYPeBB4V+wBxpgXjTHdkS/XAzOdHaY6\nm/m0/YAroqWQWi2j4rBzVswA6mK+ro98L5GPA0+NZFBqbOlvP6A5dycFtCukGoLPxjHxzpy4UzAR\n+RBQDVyW4Oe3AbcBzJo1y+YQ1dkuS9MyrrDaKOtOTCoeOzP3eqAi5uuZQMPAg0TkKuAbwI3GmN54\nT2SMudcYU22MqS4rK0tlvOos5NUFVVfozF0NxU5w3wDMF5HZIpINvB94LPYAEVkB/IJwYG9yfpjq\nbNZfCqnB3UlBLYVUQxj2rDDGBIDPAOuAXcDvjTE1InKXiNwYOexuoBD4g4hsFZHHEjydGof6q2U0\n5+4k645fnbmreOzk3DHGPAk8OeB7d8b8+yqHx6XGEE3LuCOo7QfUEPR6TrkuSxuHuSK6E5OmZVQc\nelYo12njMHdo4zA1FA3uynXaOMwd2n5ADUWDu3Jd9A5VrXN3lNV+IEvvUFVx6FmhXGfdHq8Lqs4K\nRjfI1pm7GkyDu3KdNXPX9gPO0py7GooGd+U6X3RBVWfuTrLuG9DgruLR4K5c59UFVVdo+wE1FA3u\nynVZ2n7AFYFQCJ9HENHgrgbT4K5c59X2A64IhIzO2lVCGtyV63SbPXcEg0bz7SohDe7KdSKC1yO6\nQbbDAiGjrQdUQnpmqLQIB3eduTvJyrkrFY8Gd5UWWR4hqKWQjgqGjG6OrRLS4K7SQmfuzvMHjW6O\nrRLSM0Olhc/r0Zy7w4JaLaOGoMFdpYXPI9o4zGH+YEjTMiohDe4qLXyalnFcMKSlkCoxDe4qLXxe\nj7YfcFggpDl3lZieGSotfB7RnZgcFtC0jBqCBneVFl6P6MzdYdp+QA1Fg7tKi3C1jAZ3JwVDRndh\nUgnpmaHSIlwto2kZJwWCOnNXiWlwV2mhNzE5LxDSnLtKTIO7Sossr9a5Oy2gpZBqCBrcVVrogqrz\nwmkZfQur+PTMUGmRpe0HHBcMGbI0LaMS0OCu0kJz7s7zh0K6oKoS0uCu0kJ7yzhP2w+ooWhwV2nh\n82j7AacFgroTk0pMzwyVFl6v4Necu6N0JyY1FA3uKi2ytFrGcboTkxqKBneVFl6PR3PuDtOdmNRQ\nbJ0ZIrJWRPaIyH4R+Wqcn+eIyO8iP39DRKqcHqg6u4X7uWtaxkm6E5MayrDBXUS8wE+Ba4AlwAdE\nZMmAwz4OtBpj5gE/Av7F6YGqs5vPq2kZp2n7ATUUOzP3C4H9xpiDxpg+4EHgXQOOeRfwP5F/PwRc\nKSJ61qmocD93De5OCgS1FFIl5rNxzAygLubreuCiRMcYYwIi0gZMBk44MUh19vN6PLT3+HnHPS9n\neihjRrifu+bcVXx2gnu8qcHAKZidYxCR24DbAGbNmmXjV6ux4oZzp3G8owdjdPbulIVTi1i7dGqm\nh6FGKTvBvR6oiPl6JtCQ4Jh6EfEBE4CWgU9kjLkXuBegurpa3+XjyIpZJfz0gyWZHoZS44ada7oN\nwHwRmS0i2cD7gccGHPMY8JHIv28GXjA6RVNKqYwZduYeyaF/BlgHeIFfGWNqROQuYKMx5jHgl8B9\nIrKf8Iz9/W4OWiml1NDspGUwxjwJPDnge3fG/LsHeK+zQ1NKKZUqXWpXSqkxSIO7UkqNQRrclVJq\nDNLgrpRSY5AGd6WUGoMkU+XoItIB7MnIL09OKWdHGwUdp3POhjGCjtNpZ8s4FxpjioY7yFYppEv2\nGGOqM/j7bRGRjTpO55wN4zwbxgg6TqedTeO0c5ymZZRSagzS4K6UUmNQJoP7vRn83cnQcTrrbBjn\n2TBG0HE6bUyNM2MLqkoppdyjaRmllBqDRkVwF5EviYgRkdJMjyUeEflnEdkuIltF5BkRmZ7pMcUj\nIneLyO7IWB8RkYmZHtNAIvJeEakRkZCIjLrKhOE2gx8NRORXItIkIjsyPZahiEiFiLwoIrsif/PP\nZ3pM8YhIroi8KSLbIuP8p0yPKRER8YrIFhH583DHZjy4i0gF8A7gSKbHMoS7jTHLjTHnAX8G7hzu\nARnyLLDMGLMc2At8LcPjiWcHcBPwSqYHMpDNzeBHg/8G1mZ6EDYEgH8wxiwGLgY+PUpfz17g7caY\nc4HzgLUicnGGx5TI54Fddg7MeHAHfgT8I3G25RstjDHtMV8WMErHaox5xhgTiHy5nvCuWaOKMWaX\nMWa03rxmZzP4jDPGvEKcnc5GG2PMMWPM5si/OwgHpRmZHdVgJqwz8mVW5H+j7j0uIjOB64D/snN8\nRoO7iNwIHDXGbMvkOOwQke+JSB1wK6N35h7rY8BTmR7EWSbeZvCjLhidjUSkClgBvJHZkcQXSXds\nBZqAZ40xo3GcPyY8EQ7ZOdj1O1RF5Dkg3i6+3wC+DrzT7THYMdQ4jTGPGmO+AXxDRL4GfAb4dloH\nGDHcOCPHfIPwJfFv0jk2i50xjlK2NnpXyRGRQuCPwBcGXAWPGsaYIHBeZJ3qERFZZowZNWsaInI9\n0GSM2SQil9t5jOvB3RhzVbzvi8g5wGxgm4hAOIWwWUQuNMY0uj2ugRKNM47fAk+QoeA+3DhF5CPA\n9cCVmdrHNonXcrSxsxm8SoKIZBEO7L8xxjyc6fEMxxhzSkReIrymMWqCO7AKuFFErgVygWIRud8Y\n86FED8hYWsYY85YxptwYU2WMqSL8xjo/E4F9OCIyP+bLG4HdmRrLUERkLfAV4EZjTHemx3MWsrMZ\nvLJJwrO2XwK7jDH3ZHo8iYhImVVZJiJ5wFWMsve4MeZrxpiZkVj5fuCFoQI7jI4F1bPB90Vkh4hs\nJ5xGGpUlXcB/AEXAs5GyzZ9nekADici7RaQeuAR4QkTWZXpMlshitLUZ/C7g98aYmsyOajAReQB4\nHVgoIvUi8vFMjymBVcCHgbdHzsetkZnnaDMNeDHy/t5AOOc+bKnhaKd3qCql1BikM3ellBqDNLgr\npdQYpMFdKaXGIA3uSik1BmlwV0opBzjd0E1EfhBpZLZLRP4tUlpqmwZ3pZRyxn/jUEM3EbmUcCnp\ncmAZcAFwWTLPocFdKaUcEK+hm4jMFZGnRWSTiLwqIovsPh3hO1GzgRzCzcyOJzMeDe5KKeWee4HP\nGmNWAl8CfmbnQcaY14EXgWOR/60zxthq9WtxvbeMUkqNR5GGaZcCf4hJl+dEfnYTcFechx01xlwt\nIvOAxfS37X5WRNZErg5s0eCulFLu8ACnIpv8nCHSRG2oRmrvBtZbfeZF5CnCG57YDu6allFKKRdE\n2hsfEpH3QriRmoica/PhR4DLRMQX6ax5GTZ3YLJocFdKKQckaOh2K/BxEdkG1GB/Z6+HgAPAW8A2\nYJsx5vGkxqONw5RSauzRmbtSSo1BGtyVUmoM0uCulFJjkAZ3pZQagzS4K6XUGKTBXSmlxiAN7kop\nNQZpcFdKqTHo/wOSWY385p5migAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7d7a3bd630>"
      ]
     },
     "metadata": {},
     "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",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "% matplotlib inline\n",
    "\n",
    "def f(x):\n",
    "    return (1 - np.cos(x))/(x*x)\n",
    "\n",
    "x = np.linspace(-4e-8, 4e-8, 100)\n",
    "plt.plot(x,f(x));\n",
    "plt.axvline(1.1e-8, color='red')\n",
    "plt.xlim([-4e-8, 4e-8]);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "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))) # exact answer is 6.05e-17\n",
    "print ('%2f' % ((1 - np.cos(1.1e-8))/(1.1e-8*1.1e-8)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEJCAYAAAB4yveGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAEKFJREFUeJzt3X+s3Xddx/Hnay3lxxgBXYmjLblV\nl8EyxnCXOiBCgw6Kkjb8kiEiVXCYWIcEopszI24xUTBijDNadDAV2WBKvBuDsiENRDfsLa7bulJW\nJtBL0RXq+BnYyt7+cU/n2d1pz/fenvbc7vN8JDe93+/5fO9532732W+/50dTVUiS2nDSuAeQJB0/\nRl+SGmL0JakhRl+SGmL0JakhRl+SGmL0JakhRl+SGmL0JakhS8c9wFynnnpqTUxMjHsM6dFj9+7Z\nX884Y7xz6Jjavn3716tq+bB1iy76ExMTTE9Pj3sM6dFj7drZX7duHecUOsaSfLnLOi/vSFJDjL4k\nNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcTo\nS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNcToS1JDjL4kNaRT9JOsS7I7yZ4kFw+4\nfWOS/Ulu6328ubf/nCS3JNmZ5PYkrx31NyBJ6m7psAVJlgBXAucDM8C2JFNVddecpddW1aY5+74H\n/EpV3Z3kacD2JFuq6r5RDC9Jmp8uZ/prgD1VdU9V3Q9cA2zo8sWr6gtVdXfv833AvcDyhQ4rSTo6\nXaK/Atjbtz3T2zfXq3qXcK5LsmrujUnWAMuALy5oUknSUesS/QzYV3O2rwcmqups4Gbg6od9geQ0\n4O+BX62qBx9xB8mFSaaTTO/fv7/b5JKkeesS/Rmg/8x9JbCvf0FVfaOqftDbfC9w7qHbkjwJ+Cjw\n+1V166A7qKrNVTVZVZPLl3v1R5KOlS7R3wacnmR1kmXABcBU/4Lemfwh64Fdvf3LgI8Af1dVHx7N\nyJKkhRr67J2qOphkE7AFWAJcVVU7k1wOTFfVFHBRkvXAQeAAsLF3+C8CLwR+NMmhfRur6rbRfhuS\npC5SNffy/HhNTk7W9PT0uMeQHj3Wrp39devWcU6hYyzJ9qqaHLbOV+RKUkOMviQ1xOhLUkOMviQ1\nxOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhL\nUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOM\nviQ1xOhLUkM6RT/JuiS7k+xJcvGA2zcm2Z/ktt7Hm/tu+3iS+5LcMMrBJUnzt3TYgiRLgCuB84EZ\nYFuSqaq6a87Sa6tq04Av8W7gCcBbjnZYSdLR6XKmvwbYU1X3VNX9wDXAhq53UFWfBL69wPkkSSPU\nJforgL192zO9fXO9KsntSa5Lsmok00mSRqpL9DNgX83Zvh6YqKqzgZuBq+czRJILk0wnmd6/f/98\nDpUkzUOX6M8A/WfuK4F9/Quq6htV9YPe5nuBc+czRFVtrqrJqppcvnz5fA6VJM1Dl+hvA05PsjrJ\nMuACYKp/QZLT+jbXA7tGN6IkaVSGPnunqg4m2QRsAZYAV1XVziSXA9NVNQVclGQ9cBA4AGw8dHyS\nzwDPAJ6YZAZ4U1VtGf23IkkaZmj0AarqRuDGOfsu6/v8EuCSwxz7M0czoCRpdHxFriQ1xOhLUkOM\nviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1\nxOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhLUkOMviQ1xOhL\nUkOMviQ1xOhLUkM6RT/JuiS7k+xJcvGA2zcm2Z/ktt7Hm/tue2OSu3sfbxzl8JKk+Vk6bEGSJcCV\nwPnADLAtyVRV3TVn6bVVtWnOsT8CvBOYBArY3jv2f0cyvSRpXrqc6a8B9lTVPVV1P3ANsKHj138p\ncFNVHeiF/iZg3cJGlSQdraFn+sAKYG/f9gzw0wPWvSrJC4EvAG+rqr2HOXbFke7snv3f5bV/fUuH\nsSR1cdm+bwFwuT9XotuZfgbsqznb1wMTVXU2cDNw9TyOJcmFSaaTTD/wwAMdRpIkLUSXM/0ZYFXf\n9kpgX/+CqvpG3+Z7gT/uO3btnGO3zr2DqtoMbAaYnJysa9/yvA5jSerkg08CwJ+rR7cP/Ua3dV3O\n9LcBpydZnWQZcAEw1b8gyWl9m+uBXb3PtwAvSfKUJE8BXtLbJ0kag6Fn+lV1MMkmZmO9BLiqqnYm\nuRyYrqop4KIk64GDwAFgY+/YA0muYPYPDoDLq+rAMfg+JEkdpOoRl9jHanJysqanp8c9hvTosXbt\n7K9bt45zCh1jSbZX1eSwdb4iV5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSFGX5Ia\nYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQl\nqSFGX5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSFGX5IaYvQlqSFGX5Ia0in6SdYl2Z1kT5KLj7Du\n1UkqyWRve1mS9yW5I8mOJGtHNLckaQGWDluQZAlwJXA+MANsSzJVVXfNWXcKcBHw2b7dvw5QVc9K\n8lTgY0meW1UPjuobkCR11+VMfw2wp6ruqar7gWuADQPWXQG8C/h+374zgU8CVNW9wH3A5FFNLEla\nsC7RXwHs7due6e17SJLnAKuq6oY5x+4ANiRZmmQ1cC6w6ijmlSQdhaGXd4AM2FcP3ZicBLwH2Dhg\n3VXAM4Fp4MvAvwMHH3EHyYXAhQBPf/rTO4wkSVqILmf6Mzz87HwlsK9v+xTgLGBrki8B5wFTSSar\n6mBVva2qzqmqDcCTgbvn3kFVba6qyaqaXL58+UK/F0nSEF2ivw04PcnqJMuAC4CpQzdW1Ter6tSq\nmqiqCeBWYH1VTSd5QpKTAZKcDxyc+wCwJOn4GXp5p6oOJtkEbAGWAFdV1c4klwPTVTV1hMOfCmxJ\n8iDwVeANoxhakrQwXa7pU1U3AjfO2XfZYdau7fv8S8AZCx9PkjRKviJXkhpi9CWpIUZfkhpi9CWp\nIUZfkhpi9CWpIUZfkhpi9CWpIUZfkhpi9CWpIUZfkhpi9CWpIUZfkhpi9CWpIUZfkhpi9CWpIUZf\nkhpi9CWpIUZfkhpi9CWpIUZfkhpi9CWpIUZfkhpi9CWpIUZfkhpi9CWpIUZfkhpi9CWpIUZfkhpi\n9CWpIZ2in2Rdkt1J9iS5+AjrXp2kkkz2th+T5OokdyTZleSSUQ0uSZq/odFPsgS4EngZcCbwuiRn\nDlh3CnAR8Nm+3a8BHltVzwLOBd6SZOLox5YkLUSXM/01wJ6quqeq7geuATYMWHcF8C7g+337Cjg5\nyVLg8cD9wLeObmRJ0kJ1if4KYG/f9kxv30OSPAdYVVU3zDn2OuC7wNeArwB/UlUHFj6uJOlodIl+\nBuyrh25MTgLeA7x9wLo1wA+BpwGrgbcn+fFH3EFyYZLpJNP79+/vNLgkaf66RH8GWNW3vRLY17d9\nCnAWsDXJl4DzgKneg7m/BHy8qh6oqnuBfwMm595BVW2uqsmqmly+fPnCvhNJ0lBdor8NOD3J6iTL\ngAuAqUM3VtU3q+rUqpqoqgngVmB9VU0ze0nnxZl1MrN/IHx+5N+FJKmTodGvqoPAJmALsAv4UFXt\nTHJ5kvVDDr8SeCJwJ7N/eLyvqm4/ypklSQu0tMuiqroRuHHOvssOs3Zt3+ffYfZpm5KkRcBX5EpS\nQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+\nJDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ4y+JDXE6EtSQ1JV457h\nYZJ8G9g97jk6OBX4+riH6MA5R8s5R+tEmPNEmBHgjKo6ZdiipcdjknnaXVWT4x5imCTTzjk6zjla\nzjk6J8KMMDtnl3Ve3pGkhhh9SWrIYoz+5nEP0JFzjpZzjpZzjs6JMCN0nHPRPZArSTp2FuOZviTp\nGFnU0U/yjiSV5NRxzzJIkiuS3J7ktiSfSPK0cc80SJJ3J/l8b9aPJHnyuGcaJMlrkuxM8mCSRfVs\niSTrkuxOsifJxeOe53CSXJXk3iR3jnuWw0myKsmnkuzq/fd+67hnGiTJ45L8R5IdvTn/YNwzHUmS\nJUn+M8kNR1q3aKOfZBVwPvCVcc9yBO+uqrOr6hzgBuCycQ90GDcBZ1XV2cAXgEvGPM/h3Am8Evj0\nuAfpl2QJcCXwMuBM4HVJzhzvVIf1fmDduIcY4iDw9qp6JnAe8JuL9PfzB8CLq+rZwDnAuiTnjXmm\nI3krsGvYokUbfeA9wO8Ai/ZBh6r6Vt/mySzSWavqE1V1sLd5K7BynPMcTlXtqqrF+MK8NcCeqrqn\nqu4HrgE2jHmmgarq08CBcc9xJFX1tar6XO/zbzMbqhXjneqRatZ3epuP6X0syp/xJCuBXwD+Ztja\nRRn9JOuBr1bVjnHPMkySP0yyF3g9i/dMv9+vAR8b9xAnmBXA3r7tGRZhpE5ESSaA5wCfHe8kg/Uu\nmdwG3AvcVFWLck7gz5g9SX5w2MKxvSI3yc3Ajw246VLg94CXHN+JBjvSnFX1L1V1KXBpkkuATcA7\nj+uAPcPm7K25lNm/Wn/geM7Wr8uci1AG7FuUZ3wnkiRPBP4J+O05f2teNKrqh8A5vcfBPpLkrKpa\nVI+XJHk5cG9VbU+ydtj6sUW/qn5u0P4kzwJWAzuSwOyliM8lWVNV/30cRwQOP+cA/wh8lDFFf9ic\nSd4IvBz42Rrj83Tn8fu5mMwAq/q2VwL7xjTLo0KSxzAb/A9U1T+Pe55hquq+JFuZfbxkUUUfeAGw\nPsnPA48DnpTkH6rqlwctXnSXd6rqjqp6alVNVNUEsz9wPzWO4A+T5PS+zfXA58c1y5EkWQf8LrC+\nqr437nlOQNuA05OsTrIMuACYGvNMJ6zMns39LbCrqv503PMcTpLlh57pluTxwM+xCH/Gq+qSqlrZ\n6+UFwL8eLviwCKN/gvmjJHcmuZ3Zy1GL8qlnwF8ApwA39Z5e+lfjHmiQJK9IMgM8D/hoki3jngmg\n9yD4JmALsw86fqiqdo53qsGSfBC4BTgjyUySN417pgFeALwBeHHv/8fbemepi81pwKd6P9/bmL2m\nf8SnQ54IfEWuJDXEM31JaojRl6SGGH1JaojRl6SGGH1JOoZG/SZ4Sd7VewO4XUn+vPcU2M6MviQd\nW+9nRG+Cl+T5zD7l9WzgLOC5wIvm8zWMviQdQ4PeBC/JTyT5eJLtST6T5Bldvxyzr7pdBjyW2TeB\n+5/5zGP0Jen42wz8VlWdC7wD+MsuB1XVLcCngK/1PrZU1dC3U+43tvfekaQW9d5o7vnAh/suxz+2\nd9srgcsHHPbVqnppkp8Ensn/vz36TUle2PvbRCdGX5KOr5OA+3r/+NLD9N587khvQPcK4NZD7/Of\n5GPM/kM0naPv5R1JOo56byP9X0leA7NvQJfk2R0P/wrwoiRLe+9U+iI6/GtZ/Yy+JB1Dh3kTvNcD\nb0qyA9hJ93+J7Trgi8AdwA5gR1VdP695fMM1SWqHZ/qS1BCjL0kNMfqS1BCjL0kNMfqS1BCjL0kN\nMfqS1BCjL0kN+T8CK3Vez+/6+gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7d780a3ac8>"
      ]
     },
     "metadata": {},
     "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-8, 4e-8, 100)\n",
    "plt.plot(x,f1(x));\n",
    "plt.axvline(1.1e-8, color='red')\n",
    "plt.xlim([-4e-8, 4e-8]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stable and unstable version of the variance:\n",
    "\n",
    "$s^2 = \\frac{1}{n-1} \\sum (x-\\bar{x})^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# sum of squares method (vectorized version)\n",
    "# watch out! big number minus big number!\n",
    "def sum_of_squers_var(x):\n",
    "    n = len(x)\n",
    "    return (1.0/(n*(n-1))*(n*np.sum(x**2) - (np.sum(x))**2))\n",
    "\n",
    "# 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",
    "\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 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0804168234022\n",
      "0.0\n",
      "0.080496349716\n",
      "0.0798072579566\n"
     ]
    }
   ],
   "source": [
    "# check the performances with an array \n",
    "# of randomly distributed data around 1e12\n",
    "\n",
    "x_ = np.random.uniform(0,1,int(1e3))\n",
    "x = 1e12 + x_\n",
    "\n",
    "# correct answer\n",
    "print (np.var(x_))\n",
    "\n",
    "print (sum_of_squers_var(x))\n",
    "print (direct_var(x))\n",
    "print (welford_var(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}