{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [ "# Set up environment and import packages\n", "import numpy as np\n", "import math\n", "import struct" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## MA934\n", "\n", "## Floating point arithmetic\n", "\n", "### How computers approximate arithmetic" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split", "slideshow": { "slide_type": "slide" } }, "source": [ "### Representation of unsigned integers\n", "\n", "Binary representation of (3 bit) integers:\n", "$$ b_2b_1b_0 = b_2\\,2^2 + b_1\\,2^1 + b_0\\,2^0 $$\n", "\n", "Finite maximum and minimum integers.\n", "Modern computers use 64 bits but smaller (and larger) integer types can often be accommodated for depending on the programming language of choice.\n", "\n", "Here are the unsigned integer types:\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split" }, "source": [ "|Type |Bits |Minimum |Maximum |\n", "|:- |--- | --- |---|\n", "|UInt8 | 8 |0 | $2^8 - 1$|\n", "|UInt16 | 16 |0 | $2^{16} - 1$|\n", "|UInt32 | 32 |0 | $2^{32} - 1$|\n", "|UInt64 | 64 |0 | $2^{64} - 1$|" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split", "slideshow": { "slide_type": "slide" } }, "source": [ "### Representation of signed integers : two's complement\n", "The negative of $x>0$ is encoded using *two's complement*,\n", "\n", "> $\\overline{x}$ = flip all bits of $x$ and add 1.\n", "\n", "Example: \n", "6 = 00000110, \n", "$\\overline{6}$ = 11111001 + 1 = 11111010 " ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split" }, "source": [ "Signed integer types in Python:\n", "\n", "|Type |Bits |Minimum |Maximum |\n", "|:- |--- | --- |---|\n", "|int8 | 8 |$-2^7$ | $2^7 - 1$|\n", "|int16 | 16 |$-2^{15}$ | $2^{15} - 1$|\n", "|int32 | 32 |$-2^{31}$ | $2^{31} - 1$|\n", "|int64 | 64 |$-2^{63}$ | $2^{63} - 1$|\n", "\n", "You should (in time) become comfortable with checking documentation and pages such as the [docs.python.org](https://docs.python.org/3/library/stdtypes.html) pages and [other resources](https://data-apis.org/array-api/latest/API_specification/data_types.html)." ] }, { "cell_type": "markdown", "metadata": { "cell_style": "center", "slideshow": { "slide_type": "slide" } }, "source": [ "### Advantages of two's complement \n", "\n", "Why not just use a \"sign\" bit?\n", "* Avoids two representations of zero.\n", "* Subtraction can be performed using the same hardware as addition:\n", "> Subtraction of y from x : add two's complement of y to x and drop leading (\"overflow\") bit.\n", "\n", "* Example: $7 - 6 = 7 + \\overline{6}$ = (check) " ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split", "slideshow": { "slide_type": "slide" } }, "source": [ "### Representation of real numbers\n", "\n", "Use a binary version of *normalised* scientific notation:\n", "$ x = -1^S \\times (1.0+0.M) \\times 2^E $\n", "\n", "e.g. IEEE 754 32-bit (Float32): \n", "\n", "|Field |Size|Bits |\n", "| :-- | :-: | :-: |\n", "| Sign (S) | 1 | 31 |\n", "| Exponent (E + 127) | 8 | 23 - 30 |\n", "| Mantissa (M) | 23 | 0 - 22 |" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split", "slideshow": { "slide_type": "-" } }, "source": [ "\"Bias\": if exponent is $E$, we store E+127 (for Float32). This makes comparisons easier..\n", "\n", "![Float32](files/images/590px-Float_example.png)\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "cell_style": "split", "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "\n", "1100100\n", "0000000001100100\n" ] } ], "source": [ "x1 = 2\n", "print(type(x1))\n", "\n", "x2 = 20.1\n", "print(type(x2))\n", "\n", "# Hint: have a look at Python datatypes and in particular at the 'float32' 'and float64' datatypes\n", "\n", "# We can also see our data as binary strings, either in compact form ...\n", "x3 = 100\n", "binary3 = \"{0:b}\".format(x3)\n", "print(binary3)\n", "\n", "# or with a pre-specified precision\n", "binary4 = '{0:016b}'.format(x3)\n", "print(binary4)" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split", "slideshow": { "slide_type": "slide" } }, "source": [ "### Round-off error\n", "\n", "* Finite mantissae introduce errors when truncating real numbers whose binary expansion is longer than 23 (Float32) or 52 (Float64).\n", "* Round-off error is a feature of the hardware and cannot be avoided." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "cell_style": "split" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "a = 1.0000000000000001e-01 b = 2.0000000000000001e-01 c = 2.9999999999999999e-01\n" ] } ], "source": [ "a = np.float64(0.1);\n", "b = np.float64(0.2);\n", "c = np.float64(0.3);\n", "\n", "print('a = %.16e' % a, 'b = %.16e' % b, 'c = %.16e' % c)" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "split" }, "source": [ "This can lead to counter-intuitive behaviour:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "cell_style": "split" }, "outputs": [ { "data": { "text/plain": [ "False" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a + b == c" ] }, { "cell_type": "markdown", "metadata": { "cell_style": "center", "slideshow": { "slide_type": "slide" } }, "source": [ "### Floating point arithmetic\n", "\n", "Rules for adding two floating point numbers:\n", "\n", "1. Rewrite the smaller number so its exponent matches that of the larger number.\n", "2. Add the mantissae.\n", "3. Normalise the sum.\n", "4. Round the sum.\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Machine precision " ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "cell_style": "center" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "False\n", "False\n", "True\n" ] }, { "data": { "text/plain": [ "1.000000013351432e+22" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a = np.float32(1.0); \n", "b = np.float32(pow(10.0,-10));\n", "\n", "print(a == 0)\n", "print(b == 0)\n", "print((a+b) == a)\n", "\n", "b*pow(10.0, 32)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The *machine precision* (or *machine epsilon*) is the smallest floating point number which when added to 1 gives an answer larger than one." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Machine precision\n", "\n", "In Python ``.eps`` gives the machine precision:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.1920929e-07\n", "2.220446049250313e-16\n" ] } ], "source": [ "print(np.finfo(np.float32).eps)\n", "print(np.finfo(np.float64).eps)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These values are $2.0^{-23}$ and $2.0^{-52}$, respectively." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Loss of significance\n", "Subtraction of \"nearby\" numbers leads to loss of precision.\n", "\n", "Calculate sum and difference of 2 nearby numbers in 32-bit precision:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "32 bit: a = 2.000001907348633 b = 2.0\n", "32 bit: sum = 4.000001907348633 diff = 1.9073486328125e-06\n" ] } ], "source": [ "# Interestingly - single precision is not actually supported in Python in earnest, see:\n", "# https://python-reference.readthedocs.io/en/latest/docs/float/\n", "\n", "a32o = np.float32(2.0) + pow(np.float32(3.0),-12.0); \n", "b32o = np.float32(2.0);\n", "\n", "a32 = struct.unpack('f', struct.pack('f', a32o))[0];\n", "b32 = struct.unpack('f', struct.pack('f', b32o))[0];\n", "\n", "sum32 = a32 + b32;\n", "diff32 = a32 - b32;\n", "\n", "print('32 bit: a = ', a32, 'b = ', b32)\n", "print('32 bit: sum = ', sum32, 'diff = ', diff32)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Loss of significance\n", "\n", "Now calculate the same sum and difference in 64 bit precision:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "64 bit: a = 2.000001881676423 b = 2.0\n", "64 bit: sum = 4.000001881676424 diff = 1.8816764231210925e-06\n" ] } ], "source": [ "a64 = np.float64(2.0) + pow(np.float64(3.0),-12.0); \n", "b64 = np.float64(2.0);\n", "\n", "sum64 = a64 + b64;\n", "diff64 = a64 - b64;\n", "\n", "print('64 bit: a = ', a64, 'b = ', b64)\n", "print('64 bit: sum = ', sum64, 'diff = ', diff64)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Loss of significance\n", "\n", "We can treat the 64 bit answer as \"exact\" compared to the 32 bit answer and calculate the relative error in the 32 bit result:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Sum relative error = 6.41804929265656e-09\n", "Difference relative error = 0.013643264790885564\n" ] } ], "source": [ "relErrSum = abs(sum64 - sum32)/sum64;\n", "relErrDiff = abs(diff64 - diff32)/diff64;\n", "print('Sum relative error = ', relErrSum)\n", "print('Difference relative error = ', relErrDiff)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Numerical instability\n", "\n", "Iterative calculations can lose precision by accumulation of round-off error:\n", "\n", "* Assuming inputs $\\sim 1$, each FP addition introduces an error of $\\sim \\epsilon_m$.\n", "* Might expect, after $n$ iterative steps, total error $\\sim \\sqrt{n}\\, \\epsilon_m$?\n", "\n", "However, some iterations can produce total error $\\sim e^n$ due to *instability*." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Numerical instability: simple example\n", "\n", "Consider the following iterative procedures:\n", "\n", "P1 : $a_{n+1} = \\phi\\,a_n$ with $a_0 = 1$, \n", "\n", "P2 : $a_{n+1} = a_{n-1} - a_{n}$ with $a_0 = 1$ and $a_1 = \\phi$,\n", "\n", "where $\\phi = (\\sqrt{5}-1)/2$.\n", "\n", "Both have the exact solution (check): \n", "$$a_n = \\phi^n.$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Numerical instability: simple example\n", "\n", "However their numerical behaviour is very different for large $n$:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "phi = np.float64(0.5*(math.sqrt(5.)-1.)); \n", "n = 70;\n", "\n", "# Allocate some arrays\n", "P1 = [np.float64(0.0)] * n;\n", "P2 = [np.float64(0.0)] * n;\n", "\n", "# Set initial conditions\n", "P1[0] = P2[0] = np.float64(1.); \n", "P2[1] = phi;\n", "\n", "# Note that indexing for lists start from 0 in Python!" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "for i in range(1, n):\n", " P1[i] = phi*P1[i-1]" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "for i in range(2, n):\n", " P2[i] = P2[i-2] - P2[i-1]" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Numerical instability: simple example" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwYAAAHzCAYAAACaBgpXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABHyElEQVR4nO3df2xUZ57v+c9jl8GEiqkQ4woXE+4U444Q08pg0g1eMWrcbW7/k9awc0lYtXbZXnFjB93bpkEQT6S7YPinlwaZbnc0xE6j1dK6yhKYGaLMH5FwL0RCgtD8GG24jCKW6glxejA47kqlECYYnv2jTpmyfewq23Xq1I/3SyqFOqdc9fVjB+pb53mej7HWCgAAAEB5q/C7AAAAAAD+ozEAAAAAQGMAAAAAgMYAAAAAgGgMAAAAAIjGAAAAAICkgN8FlLJQKGT//M//3O8yyta9e/c0f/58v8soS4y9vxh/fzH+/mHs/cX4++vy5cuD1tpFs3kOGgMPhcNhXbp0ye8yytbZs2e1fv16v8soS4y9vxh/fzH+/mHs/cX4+8sY89lsn4OpRAAAAABoDAAAAADQGAAAAAAQjQEAAAAA0RgAAAAAEI0BAAAAANEYAAAAABCNAQAAAADRGAAAAAAQjQEAAAAA0RgAAAAAEI0BAAAAANEYAAAAABCNAQAAAADRGAAAAAAQjQEAAAAA0RgAAAAAEI0BAAAAANEYAAAAABCNAQAAAADRGAAAAAAQjQEAAAAA0RgAAAAAEI0BAAAAANEYAAAAABCNAQAAAABJAb8LKBbGmDckRSUtlBS11vb5XBIAAACQMzQGWTDG9Eg6kWoGjDGnjTGXrLUxfysDAAAAcoOpRNl5ddwVgtOSXvWrGAAAACDXaAwyMMY0ShoadzgmaUP+qwEAAAC8UfJTiYwxmyRttta+Msn5kKQ3Jd10Di231nakPSSiZCOQbsg5DgAAAJSEkm0MnHUBUvIN/MIpHnpCUpu1Nup8XcQYc9pam+mKQGj2VQIAAACFoWSnEllr26y1bUq+8XflXE2IppoC5+uiaedSQuO+dKEmXkXIynP1z8sYk/H2XP3zM3l6AAAAYEZK9opBljYruZB4vNOS2iSdVLIBGH/FIaTk1qXTNvDF59LO30nf3FfF5eOa+8kHGr73larnL9CDb/9Ij1dvlubM00DXD2by9AAAAMCMlOwVgyy1yP0NflTSS5I0SV7Bs3JvKLLzzX3Ne++nWt9/Su+uCejaxrDeXRPQ9z4/pXnv/VT65v6MnxoAAACYibJtDJxFxyFN3HFISl4lCKXd73N2J0pplPTeTF+74vJxNc0d1N99Z75WhKoUqDBaEarSke/O19o5g6q4fHymTw0AAADMSNk2Bpp6QbKk0eZBkl6TtNkY0+IkIB+YTbjZ3E8+UPu35soYM/711P7CXM299sFMnxoAAACYkXJfY5AVpwlIbWHqNrVolDGmVVKrJC1atEhnz56d8Jjhe1+poSbs+vUNNQENJ76UJNevRfYSiQRj6BPG3l+Mv78Yf/8w9v5i/IsfjUGOWWt7JfVK0gsvvGDXr18/4THV8xfoRnxEK0JVE87diI+oOrhA97+Oye1rkb2zZ88yhj5h7P3F+PuL8fcPY+8vxr/4lfNUIre1BWPMZrrQVB58+0f61acPZK0d/3rq/vSBHvzFj7x4WQAAAGBSZdsYOG/6Y3JPMHZLO86Zx6s368I3tdp28Z6uxx7q4WOr67GH2nbxni58U5vcshQAAADIo3KfSnRJ7ouQlyvDWoJZmTNP91/9tT66fFwXLn6g4cSXqg4u0IO/2JhsCo5tlaQJi5PHCy9Zqtv9tzwrEwAAAOWj3BuDE5I2yFkTkKZF0s+9eMHwkqWj4WWPJaUSC+5/HZPO/zZ5kwhBAwAAQF6Vw1SikMZmEoxyFgpHjDGj04mcvIIha+1JL4q53X9L1topb5IIQQMAAEBelewVA2PMASUbglclhYwxJ5RccNxjrb2S9tAfSHrTGHPTub/cWrshr8W6SA9BS00pSoagBbTt4qA+unxcj32uEQAAAKWjZBsDa20qd6Atw+NiepJRUDDmfvKB2tdMHoJ24eIH4poBAAAAcqUcphIVpWQImnvflgxB+yrPFQEAAKCU0RgUqFQImptUCBoAAACQKzQGBYoQNAAAAOQTjUGByioErbJKxpiMt+fqn/f72wEAAECBK9nFx0UvUwjanHnSo4fkHQAAACAnaAwKUNYhaNJo3kHT3EG1r5mrhpqwbsRH9KtPT+nCzXO6/+qv810+AAAAihBTiQpQNiFoqbUH6XkHK0JVClQYJ+9gvtbOGVTF5eM+fzcAAAAoBjQGRW7uJx+o/VuT5x3MvfaBT5UBAACgmNAYFDnyDgAAAJALNAZFjrwDAAAA5AKNQZEj7wAAAAC5QGNQ5DLmHVz7UJLIOgAAAMCU2K602GXKO3jrZbIOAAAAkBGNQRHLOu+ArAMAAABkwFSiIpZN3oFE1gEAAAAyozEoA2QdAAAAIBMagzJA1gEAAAAyoTEoA2QdAAAAIBMagzJA1gEAAAAyoTEoAxmzDlZv9rtEAAAA+IztSstBpqyDOfOkyqoJi5PdhJcs1e3+W3koGgAAAPlEY1Diss46kAhCAwAAKGM0BiUu20/3jTEEoQEAAJQx1hhgFEFoAAAA5YvGAKMIQgMAAChfNAYYRRAaAABA+aIxwCiC0AAAAMoXjQFGEYQGAABQvmgMMCqrIDQn7yDT7bn65/3+dgAAADANbFeKJ7IJQnv0kLwDAACAEkRjAEnTDEIj7wAAAKDkMJUIkpJBaNbajDeJvAMAAIBSRGOAaSPvAAAAoPTQGGDayDsAAAAoPTQGmDbyDgAAAEoPjQGmjbwDAACA0kNjgGnLmHdw7UNJIusAAACgiLBdKaYvU97BWy+TdQAAAFBkaAwwLVnnHZB1AAAAUFSYSoRpySbvQCLrAAAAoNjQGMATZB0AAAAUFxoDeIKsAwAAgOJCYwBPkHUAAABQXGgM4AmyDgAAAIoLjQE8kTHrYPVmqbIqY9YBeQcAAAD5wXal8EamrIM586RHD8k7AAAAKBA0Bsi5rLMOJPIOAAAACgRTiZBz2WQdkHcAAABQWGgM4CvyDgAAAAoDjQF8Rd4BAABAYaAxgK/IOwAAACgMNAbwFXkHAAAAhYHGAL7KKu8AAAAAnmO7UvgrU97Bsa2SNGFx8njhJUt1u/9WPioGAAAoSTQG8E3WeQeEoAEAAHiOxgC+yeYTfmMMIWgAAAB5wBoDFDxC0AAAALxHY4CCRwgaAACA92gMUPAIQQMAAPAejUEWjDEhY0yrMeYNY8wJY0yL3zWVE0LQAAAAvEdjkJ03rbW91tpfSHpN0mljTMTvosoFIWgAAADeozHIwGkAQqn71tqYpJOSOnwqqexkFYJWWSVjzJhbc3PzhGPP1T/v97cDAABQkNiuNDuvGmM6nKYgZaFfxZSdTCFoc+ZJjx6SdwAAADALRdcYGGM2SdpsrX1lkvMhSW9KuukcWm6tnfGn+9baqKRnxh1ulNQz0+dE9rIOQZPIOwAAAJiFomkMjDGpN+IRTf1p/QlJbc4behljIsaY09baDTmqo1GSnPUG8Fg2IWhScuvS9LyD1NamybyDgLZdHNRHl4/rsZfFAgAAFLGiWWNgrW2z1rYp+cbflXM1IZpqCpyvi6ady4V3JOWkyUBukXcAAAAwc0VzxSBLmyWddjl+WlKbkouGZYxplbQ6w3OdttaeTD9gjDkg6bX0xgOFI5l3EHY9l8w7+DLPFQEAABSPUmsMWuQ+9z8q6aXUHWtt73Sf2Gkmjltrrzj3G1N/RmFI5R2sCFVNOJfKO7j/dSz/hQEAABSBoplKlImz6DgkacjldExpW47O4LlbnOeNOmFnIaU1GigM5B0AAADMXMk0Bspi+1DnDf20ODkGp5Vc2/CntBsKTMa8g2sfStKEbAOyDgAAAEpvKlHOOesJTMYHwn+Z8g7eepmsAwAAgEnQGOSYsxahVZIWLVqks2fP+ltQmXhmUVh/yibvIMusA35us5NIJBhDHzH+/mL8/cPY+4vxL36l1Bi4rS0YY1xysSechc29kvTCCy/Y9evXe/2SkDR05/aEY2fPnlX6+E8n64Cf2+yMH3vkF+PvL8bfP4y9vxj/4lcyawycN/0xJQPQxos451DmyDoAAABwVzKNgeOS3BchL5fUl+daUICSWQfuF8qSWQdf5bkiAACAmUkkEtrfuVfLFtdJmTO6Miq1xuCE3FOJWyQdz3MtKECprAM3qawDAACAQpdIJNS8rkkfH+tW90qrlaHZrxAoxsYgpEkyCZz5/RFni1FJySAySUPjU4xRnsg6AAAApaDr0EHVxfvVtaraCXed/SaaRdMYGGMOGGN6JL2p5Jv/E8aYHueNf7ofSGozxrQ6OwRttta6XUVAGcqYdbB6s1RZlTHrgLwDAADgp6M9R9QaCUxYNzkbRbMrkbW2w/ljW4bHxSR1TPUYlLFMWQdz5kmPHpJ3AAAAClr/wKAamsI5fc6iaQyA2QovWTr6Zn7KrAMp67wDAAAAP9SHa3UjPuJMI8qNoplKBMzW7f5bstZmvEkak3ewIlSlQIVx8g7ma+2cQVVcZi07AADwz9a2beqJjkxYNzkbNAaAC/IOAABAIdu5a7fu1tRrx9VhXY89lNXsGwQaA8AFeQcAAKCQBYNBnTl3Xmu3tGv79Qpdj7lvxz4dNAaAC/IOAABAoQsGg9rTuU+f/XFAki7P9vloDAAX5B0AAIByQ2MAuMgq7wAAAKCEsF0p4CabvAMnCC2T8JKlut1/Kw9FAwAAzByNATDOtPIOCEIDAAAlgsYAGCfbT/eNMQShAQCAksEaA2AWCEIDAAClgsYAmAWC0AAAQKmgMQBmgSA0AAAwHYlEQvs792rZ4jpVVlRo2eI67e/cq0Qi4XdpNAbAbBCEBgAAspVIJNS8rkkfH+tW90qrTzaG1b3S6sKxbjWva/K9OaAxAGaBIDQAAJCtrkMHVRfvV9eq6jFrEw+vqtaieL+6Dh30tT4aA2AWMgahXftQUnLNwVS35+qf9/k7AQAAXjvac0StkYDr2sTWSEBHe9/2qbIktisFZiNTENpbL5N1AAAAJEn9A4NqaAq7nmuoCeiLgTt5rmgsGgNghrIOQiPrAAAASKoP1+pGfEQrQlUTzt2Ij2hJuNaHqp5gKhEwQ7f7b8laO+VNIusAAAAkbW3bpp7oiOvaxN7oiLa2vu5TZUk0BoDHyDoAAACStHPXbt2tqdeOq8Nj1ibuuDqsuzX12rlrt6/10RgAHiPrAAAASFIwGNSZc+e1dku7tl+v0Ivv39H26xVau6VdZ86dVzAY9LU+1hgAHktlHUw2n7A6uCC5LgEAAJS8YDCoPZ37tKdzn9+lTMAVA8BjZB0AAIBiQGMAeCxj1sHqzVJlVcasA/IOAACAl5hKBHgtU9bBnHnSo4fkHQAAAF/RGAAeyjrrQCLvAAAA+IqpRICHssk6IO8AAAAUAhoDoECQdwAAAPxEYwAUCPIOAAAoTIlEQvs792rZ4jpVVlRo2eI67e/cq0Qi4XdpOUVjABSIVN6Bm1TeAQAAyK9EIqHmdU36+Fi3uldafbIxrO6VVheOdat5XVNJNQc0BkCBIO8AAIDC03XooOri/epaVT1mDeDhVdVaFO9X16GDfpeYMzQGQIHImHdw7UNJIusAAIA8OtpzRK2RgOsawNZIQEd73/apstxju1KgUGTKO3jrZbIOAADIs/6BQTU0hV3PNdQE9MXAnTxX5B0aA6AAZJ13QNYBAAB5VR+u1Y34iFaEqiacuxEf0ZJwrQ9VeYOpREAByCbvQCLrAACAfNvatk090RHXNYC90RFtbX3dp8pyj8YAKCJkHQAAkF87d+3W3Zp67bg6PGYN4I6rw7pbU6+du3b7XWLO0BgARYSsAwAA8isYDOrMufNau6Vd269X6MX372j79Qqt3dKuM+fOKxgM+l1izrDGACgiqayDyeY5VgcXJNclAACAnAkGg9rTuU97Ovf5XYqnuGIAFBGyDgAAgFdoDIAikjHrYPVmv0sEAABFiqlEQDHJlHUwZ55UWTVhcbKb8JKlut1/Kw9FAwCAYkBjABSJrLMOJILQAADAtNEYAEUi20/3jTEEoQEAgGljjQFQgghCAwAA00VjAJQggtAAAMB00RgAJYggNABAOUskEtrfuVfLFtepsqJCyxbXaX/nXiUSCb9LK2g0BkAJSgWhuUkFoQEAUIoSiYSa1zXp42Pd6l5p9cnGsLpXWl041q3mdU00B1OgMQBKEEFoAIBy1XXooOri/epaVT1mnd3hVdVaFO9X16GDfpdYsGgMgBKUMQjt2oeSkmsOpro9V/+8z98JAADTc7TniFojAdd1dq2RgI72vu1TZYWP7UqBUpQpCO2tl8k6AACUpP6BQTU0hV3PNdQE9MXAnTxXVDxoDIASk3UQGlkHAIASVB+u1Y34iFaEqiacuxEf0ZJwrQ9VFQemEgEl5nb/LVlrp7xJZB0AAErT1rZt6omOuK6z642OaGvr6z5VVvhoDIAyRdYBAKAU7dy1W3dr6rXj6vCYdXY7rg7rbk29du7a7XeJBYvGAChTZB0AAEpRMBjUmXPntXZLu7Zfr9CL79/R9usVWrulXWfOnVcwGPS7xILFGgOgTKWyDiabg1kdXJBclwAAQJEJBoPa07lPezr3+V1KUeGKAVCmyDoAAADpaAyAMpUx62D1ZqmyKmPWAXkHAACUBqYSAeUqU9bBnHnSo4fkHQAAUCZoDIAylHXWgUTeAQAAZYKpRNNkjGk0xrzhdx3AbGSTdUDeAQAA5YXGYPoO+F0AkE/kHQAAUB5oDKbBGNMiKep3HUA+kXcAAMiHRCKh/Z17tWxxnSorKrRscZ32d+5VIpHwu7SyQWMwPSFJN/0uAsinVN6Bm1TeAQAAs5FIJNS8rkkfH+tW90qrTzaG1b3S6sKxbjWva6I5yJOiawyMMZuMMSemOB8yxhwwxrQ6t5xM/THGbLLWnszFcwHFhLwDAIDXug4dVF28X12rqsesZzu8qlqL4v3qOnTQ7xLLQtE0BsaYHmNMj6Q2SZEpHnpCUo+1ttda2yupxxhzepavHRFTiFCmMuYdXPtQksbkGjQ3N5N1AADI2tGeI2qNBFzXs7VGAjra+7ZPlZWXotmu1FrbJknGmFYlm4MJjDGbJEWttdG0r4s6b0xm84l/I1cLULYy5R289TJZBwCAWekfGFRDU9j1XENNQF8M3MlzReWpaBqDLG2W5HZ14LSSzcRJabS5WJ3huU5ba086C477clolUCSyzjsg6wAAMAv14VrdiI9oRahqwrkb8REtCdf6UFX5KbXGoEVSj8vxqKSXUnecKUbT0Zp2aWuzpCFjjKy1v5hRlUCRuN1/K+NjjDFjsg5S/68ksw4C2nZxUB9dPq7HXhcLAChaW9u2qedYtw6vGjudyFqr3uiItra2+1hd+SiaNQaZGGNCSu4aNORyOuacmzZrbZ+19hepm6RLkq7QFABPkHUAAJiNnbt2625NvXZcHR6znm3H1WHdranXzl27/S6xLJRMYyBpYaYHOM3DjDmJxy2SWpzpSABE1gEAYHaCwaDOnDuvtVvatf16hV58/462X6/Q2i3tOnPuvILBoN8lloVSm0rkKecqAVcKgHFSWQeTzQ2tDi5IrksAAGASwWBQezr3aU/nPr9LKVs0BjnmXElolaRFixbp7Nmz/hZUxhKJBOOfJ8msg1M68t2Jc0OTWQcbpfO/5eeRJ/zu+4vx9w9j7y/Gv/iVUmPgtrZgDGttzOsinIXNvZL0wgsv2PXr13v9kpjE2bNnxfjnx+PVm3Xh5jltuzio9hfmqqEmoBvxEXV/+iCZdbB6s3T+t/w88oTffX8x/v5h7P3F+Be/kmkMrLUxY0xMyfCzK+NOR5RcgAzAC5myDubMkyqrJixOdhNesjSr3ZAAAEBulUxj4Lgk90XIy0UWAeCJrLMOJILQAAAoYKW0K5EknZC0weV4i6Tjea4FKAu3+2/JWjvmdubMmQnHJI0Goa3vP6V31wR0bWNY764J6Hufn9K8934qfXN/6hcDAACeKcbGIKRJMgmc+f0RY0wkdcwY0yhpyFp7Mi/VAZhUehDailCVAhXGCUKbr7VzBlVxmf4dAAC/FM1UImPMASUbglclhYwxJ5RccNxjrU1fU/ADSW8aY24695dba92uIgDIs7mffKD2NZMHoV24+IG4ZgAAgD8mbQyMMX8m6T968JonrbX/Ot0vstZ2OH9sy/C4mKSOqR4DwB/JILSw67lkENqXea4IADAbiURCXYcO6mjPEfUPDKo+XKutbdu0c9duQsmK0FRXDIYkeRFXmnFbUQCliSA0ACgdiURCzeuaVBfvV/fKgBqawroRH1HPsW41n/oHEouL0KSNgbX2K0nv5LEWACUu2yA0AEDh6zp0UHXxfnWtqh79O31FqEqHVwW042q/ug4dJMW4yBTj4mMARerx6s268E2ttl28p+uxh3r42Op67KG2Xbz3JAjNyTvIdHuu/nm/vx0AKGtHe46oNRJwXTfWGgnoaO/bPlWGmSqaxccASkA2QWiPHpJ3AABFoH9gUA1Nk68b+2LgTp4rwmx53hg4i5jtTBYcAygd0wpCc/IOmuYOqn3NXDXUJOet/urTU7pw85zuv/rrfJcPABinPlw75bqxJeFaH6rCbHgylcgY831jzF9KkrX2D5KeMcZ834vXAlAc3ILQ3G4SeQcAUAy2tm1TT3TkSYilw1qr3uiItra+7lNlmKmcNwbGmP9DyayB3xhjvjTGHJf0Z5IiU38lACTN/eQDtX9r8ryDudc+8KkyAEDKzl27dbemXjuuDo9ZN7bj6rDu1tRr567dfpeIafLiisHvrbWvW2tfstY+K+k9Sf9BbFMKIEvJvAP3mY7JvAMvdlIGAExHMBjUmXPntXZLu7Zfr9CL7w9o+/UKrd3SzlalRcrzNQbW2r+X9Pdevw6A0kHeAQAUh2AwqD2d+7Snc5/Onj2r9evX+10SZsGLKwZ9xpifG2Oe9uC5AZSBZN7BA9d5q8m8gx/5VBkAAKXLi8bggKTlkj4zxtwwxhwxxvyNMabGg9cCUIIy5h1c+1CSyDoAACCHvJhKdNla+44kGWMWSGpRco3BZucGAFPLlHfw1stkHQAAkGM5bwyste84W5NGnewC1hgAyFrWeQdkHQAAkFM5nUpkjFlljPlPki4RaAZgJrLJO5DIOgAAINdm1BgYYxYYY/6TMebfpx+31l6VdEJSmzHmb3JQHwC4IusAAIDcmukVg4ikXkk3UyFmxpitxph/b639ylp7UASaAfAQWQcAMDOJREL7O/dq2eI6VVZUaNniOu3v3KtEIuF3afDZTNcYLJTUIemKkguLfyDpHUnWGBOTFHVuAOAJsg4AYPoSiYSa1zWpLt6v7pUBNTQl12f1HOtW86l/IJiszM30ikGLtfagtfZ31toOJ+W4QtJ3JP1G0p+stexABMAzZB0AwPR1HTqouni/ulZVj1mfdXhVtRbF+9V16KDfJcJHOV18bK29Yq3tkNRjjPnLXD43AKTLmHWwms8mAGC8oz1H1BoJuK7Pao0EdLT3bZ8qQyGYaWNgnC1JXVlr/17J/AIA8EYq62DpRv344oi+fWpAP744oo+WbkxuVTpnnlRZlTEEjSA0AOWkf2BwyvVZXwwM5rkiFJIZrTGw1v6tMeaSMeaipAPW2s9yXBcATCrrrAOJIDQASFMfrp1yfdaScK0PVaFQzGYqUYukBklRY8wNY8wRZwvTvzHG/FzS8tyUCABjZZN1MLr2wAlCW99/Su+uCejaxrDeXRPQ9z4/pXnv/VT65v7ULwYAJWRr2zb1REdc12f1Rke0tfV1nypDIZhxY2CtjVlrN0jaLOlfnf/2Sjqp5Falf5uLAgFgNghCA4Andu7arbs19dpxdXjM+qwdV4d1t6ZeO3ft9rtE+GjWi4+ttSettRustQslPWOtrbDWbrbWsok4AN8RhAYATwSDQZ05d15rt7Rr+/UKvfj+HW2/XqG1W9rZqhQzzjFwRTMAoNAkg9DCrueSQWhf5rkiAPBXMBjUns592tO5z+9SUGByul0pABSaVBCam1QQGgAAmOKKgTHmzyT9Rw9e86S19l89eF4AmCAZhHZKR747dt/uJ0FoG8fuYgQAQJmaairRkCQvpgYNefCcAODq8erNunDznLZdHFT7C3PVUBPQjfiIuj99kAxCu/ahJE1YgzBeeMlS3e6/lY+SAQDwxaSNgbNe4J081gIAuZcKQrt8XBcufqDhxJeqDi7Qg7/YmMwxeOtlsg4AAFCOFx8DQCHJOgjNyTpomjuo9jVz1VAT1o34iH716SlduHkumaQMAECJY/ExgJKVTRCaRNYBAAASjQEAkHUAAIBoDADAyTpwn1mZzDogogWAvxKJhPZ37tWyxXWqrKjQssV12t+5V4lEwu/SUEJoDACUPbIOABSyRCKh5nVN+vhYt7pXWn2yMazulVYXjnWreV0TzQFyhsYAQNlLZh08GF1zkPIk6+BHPlUGAFLXoYOqi/era1X1mHVQh1dVa1G8X12HDvpdIkoEjQGAsvd49WZd+KZW2y7e0/XYQz18bHU99lDbLt5LZh2s3ixVVskYk/H2XP3zfn87AErM0Z4jao0EXNdBtUYCOtr7tk+VodSwXSkAZMo6mDNPevSQvAMAvugfGFRDU9j1XENNQF8M3MlzRShVNAYAylrWWQcSeQcAfFEfrtWN+IhWhKomnLsRH9GScK0PVaEUMZUIQFnLJuuAvAMAftratk090RHXdVC90RFtbX3dp8pQamgMACBL5B0A8MPOXbt1t6ZeO64Oj1kHtePqsO7W1Gvnrt1+l4gSQWMAAFki7wCAH4LBoM6cO6+1W9q1/XqFXnz/jrZfr9DaLe06c+68gsGg3yWiRLDGAACylMo7mGyeb3VwQXJtAgDkWDAY1J7OfdrTuc/vUlDCuGIAAFki7wAAUMpoDAAgSxnzDq59KElkHQAAihJTiQAgW5nyDt56mawDAEDRojEAgCxknXdA1gEAoEgxlQgAspBN3oFE1gEAoHjRGABADpF1AAAoVjQGAJBDZB0ASEkkEtrfuVfLFtepsqJCyxbXaX/nXiUSCb9LA1zRGABADqWyDtyksg4AlL5EIqHmdU36+Fi3uldafbIxrO6VVheOdat5XRPNAQoSjQEA5BBZBwAkqevQQdXF+9W1qnrMeqPDq6q1KN6vrkMH/S4RmIDGAAByKGPWwerNfpcIIA+O9hxRayTgut6oNRLQ0d63faoMmBzblQJALmXKOpgzT6qsmvBmwU14yVLd7r+Vh6IB5Fr/wKAamsKu5xpqAvpi4E6eKwIyozEAgBzJOutAIggNKHH14VrdiI9oRahqwrkb8REtCdf6UBUwNaYSAUCOZJN1MLr2wAlCW99/Su+uCejaxrDeXRPQ9z4/pXnv/VT65v7ULwagoG1t26ae6IjreqPe6Ii2tr7uU2XA5GgMAMAHBKEBpW3nrt26W1OvHVeHx6w32nF1WHdr6rVz126/SwQmoDEAAB8QhAaUtmAwqDPnzmvtlnZtv16hF9+/o+3XK7R2S7vOnDuvYDDod4nABKwxAAAfJIPQJl+YOJz4Ms8VAci1YDCoPZ37tKdzn9+lAFnhigEA+IAgNABAoaExAAAfEIQGACg0NAZZMsZEjDFvGGM2GWNa/a4HQHHLKgjNyTvIdHuu/nm/vx0AQAlgjUEWjDERST3W2g3O/RPGmCFr7UmfSwNQrLIJQnv0kLwDAEDe0Bhkp0dSR9r916y1MZ9qAVDkphWE5uQdNM0dVPuauWqoCetGfES/+vSULtw8p/uv/jrf5QMAShRTiTIwxoQktVhrrxhjGo0xEZoCALMxnSA08g4AAPlSdI2BM8f/xBTnQ8aYA8aYVud2YJYvGZEUM8ZskhSVFDHG9MzyOQEgK+QdAADypWgaA2NMj/OGvE3JN+uTOaHkeoBea22vpB5jzOlZvHREUkhS1Fobs9b2SVrIAmQA+ZDMO3Cf9ZnMO/gqzxUBpS2RSGh/514tW1ynyooKLVtcp/2de5VIJPwuDfBc0TQG1to2a22bkm/8XaU+1bfWRtO+Lpp2bjavfyXtblTSK7N5PgDIBnkHQP4kEgk1r2vSx8e61b3S6pONYXWvtLpwrFvN65poDlDyiqYxyNJmSZddjp9W8kqDJMmZYtST4ZZqJKKSYuOe76akhV58AwCQjrwDIH+6Dh1UXbxfXauqx6zpObyqWovi/eo6dNDvEgFPldquRC1K7iA0XlTSS6k7zhSjrDiLjkPjDockXZpBfQAwLY9Xb9aFm+e07eKg2l+Yq4aagG7ER9T96YNk3sG1DyVpwhqE8cJLlup2/618lAwUraM9R9S9MuC6pqc1EtD23re1p3OfT9UB3iuZxsB58x6SNORyOuacm6lfGGM2peUWbNDY7UsBwBuZ8g7eepmsAyBH+gcG1dAUdj3XUBPQFwN38lwRkF8l0xgoi6k9xpjQTLYatdZ2OKnHrZKWSzowbs0BAORc1nkHZB0AOVEfrtWN+IhWhKomnLsRH9GScK0PVQH5Y8bPWy10zpvzNmvt6nHHI0rO/V89/k27MaZFyXUGz3idQeDU1ypJixYtWv3ee+95+XKYQiKRUDAY9LuMssTY509zc7Mqmv4Xre8/pb/7zvwxUyCstdp28Z4+WrpRj8//VmfOnPGx0vLB779/Zjv2x/6v/1N/PPOP+tXqif8vbb98T/+u+X/Ulv/1f8tFqSWJ331/NTc3X7bWvpT5kZOjMfDQCy+8YD/99NN8vRzGOXv2rNavX+93GWWJsc8fY4zmBUN6d03A9VPO67GH+vHFEd3/OjZhATO8we+/f2Y79qldiRbF+9UaCYyu6emNjuhuTb3OnDvPG98p8LvvL2PMrBuDUtqVyG1twRgkFgMoRWQdALkRDAZ15tx5rd3Sru3XK/Ti+3e0/XqF1m5ppylAWSiZNQbW2pgxJqZkINn4+f8RTdxyFABKQirrYLJ50dXBBcl1CQAyCgaD2tO5j92HUJZK6YqBlNxC1G0R8nJJfXmuBQDygqwDAEAulFpjcELJrUTHa5F0PM+1AEBePF69WRe+qdW2i/d0PfZQDx9bXY891LaL95JZB6s3S5VVMsZkvD1X/7zf3w4AwCfFOJUopEkyCay1vcaYNmNMxFoblSRjTKOkobQMAgAoLZmyDubMkx49JO8AADClomkMjDEHlGwIXpUUMsacUHLBcc+4XYh+IOlNY8xN5/5ya63bVQQAKHpZZx1I5B0AAKZUNI2BtTaVNNyW4XExkUoMoEzc7r814ZjbloHGGFVcPq6muYNj8g5WhKp05LsBbbs4qI8uH9fjfBQNAChIpbbGAAAwibmffKD2b80dE9wkJZuG9hfmau61D3yqDABQCGgMAKBMkHeAUpRIJLS/c6+WLa7TD77/fS1bXKf9nXuVSCT8Lg0oOjQGAFAmUnkHblJ5B0AxSSUVf3ysW90rrT7ZGFb3SqsLx7rVvK6J5gCYJhoDACgT5B2g1HQdOqi6eL+6VlVrRahKgQqjFaEqHV5VrUXxfnUdOuh3iUBRoTEAgDKRVd4BUESO9hxRayTgum6mNRLQ0d63faoMKE5FsysRAGCWMuUdHNsqSRPeZI0XXrLUdTckIN/6BwbV0BR2PddQE9AXA3fyXBFQ3GgMAKAMZJ13QAgaikh9uFY34iNaEaqacO5GfERLwrU+VAUUL6YSAUAZuN1/S9baKW+SRkPQ1vef0rtrArq2Max31wT0vc9Pad57P5W+uT/1CwF5tLVtm3qiI67rZnqjI9ra+rpPlQHFicYAADAqPQQtfTHnke/O19o5g6q4fNzvEoFRO3ft1t2aeu24Ojxm3cyOq8O6W1Ovnbt2+10iUFRoDAAAowhBQzEJBoM6c+681m5p1/brFXrx/QFtv16htVvadebceQWDQb9LBIoKjQEAYBQhaCg2wWBQezr36bM/Dqjvd/+PPvvjgPZ07qMpAGaAxgAAMIoQNAAoXzQGAIBRhKABQPmiMQAAjMoqBK2ySsaYjLfn6p/3+9sBAEwDOQYAgCcyhaDNmSc9ekjeAQCUIBoDAICkaYSgSaN5B01zB9W+Zq4aasK6ER/Rrz49pQs3z+n+q7/Od/kAgFliKhEAQFJ2IWiptQfkHQBA6aExAABMG3kHAFB6aAwAANNG3gFmI5FIaH/nXi1bXKfKigotW1yn/Z17lUgk/C4NKGs0BgCAaSPvADOVSCTUvK5JHx/rVvdKq082htW90urCsW41r2uiOQB8RGMAAJg28g4wU12HDqou3q+uVdVj1qccXlWtRfF+dR066HeJQNmiMQAATFvGvINrH0oSWQeY4GjPEbVGAq7rU1ojAR3tfdunygCwXSkAYPoy5R289TJZB3DVPzCohqaw67mGmoC+GLiT54oApNAYAACmJeu8A7IO4KI+XKsb8RGtCFVNOHcjPqIl4VofqgIgMZUIADBN2eQdSGQdwN3Wtm3qiY64rk/pjY5oa+vrPlUGgMYAAOAJsg7gZueu3bpbU68dV4fHrE/ZcXVYd2vqtXPXbr9LBMoWjQEAwBNkHcBNMBjUmXPntXZLu7Zfr9CL79/R9usVWrulXWfOnVcwGPS7RKBsscYAAOCJVNbBZHPJq4MLkusSUHaCwaD2dO7Tns59fpcCIA1XDAAAniDrAACKC40BAMATGbMOVm+WKqsyZh2QdwAA+cFUIgCANzJlHcyZJz16SN4BABQIGgMAQM5lnXUgkXcAAAWCqUQAgJzLJuuAvAMAKCw0BgAAX5F3AACFgcYAAOAr8g6KQyKR0P7OvVq2uE6VFRVatrhO+zv3KpFI+F0agByhMQAA+CqVd+AmlXcAfyUSCTWva9LHx7rVvdLqk41hda+0unCsW83rmmgOgBJBYwAA8BV5B4Wv69BB1cX71bWqesw6kMOrqrUo3q+uQwf9LhFADtAYAAB8lVXeAXx1tOeIWiMB13UgrZGAjva+7VNlAHKJ7UoBAP7KJu/ACUKTpAojza2Qhh9J1ZXSg8fSY+diQ3jJUt3uv+XjN1Oa+gcG1dAUdj3XUBPQFwN38lwRAC/QGAAAfDOtvIP/8k9P8g6+NVcNNQEn7+CBLnxTq/uv/loDb72c72+hLNSHa3UjPqIVoaoJ527ER7QkXOtDVQByjalEAADfkHdQHLa2bVNPdMR1HUhvdERbW1/3qTIAuURjAAAoCuQd+Gfnrt26W1OvHVeHx6wD2XF1WHdr6rVz126/SwSQAzQGAICiQN6Bf4LBoM6cO6+1W9q1/XqFXnz/jrZfr9DaLe06c+68gsGg3yUCyAHWGAAAikIq72Cyee7VwQXJtQnwRDAY1J7OfdrTuc/vUgB4hCsGAICiQN4BAHiLxgAAUBTIOwAAbzGVCABQHDLlHRzbKkkTFiePR9YBALijMQAAFI858/S46Se63/QTSU9yDyRJ8QFp5++kb+6r4vJxzf3kAw3f+0rV8xfowbd/NBqWlspNAACMRWMAACh46UFoU/rm/pMQtDVz1VATdkLQTunCzXO6/+qvvS8WAIoUawwAAAUvmyA0iRA0AJgNGgMAQMkgBG2sRCKh/Z17tWxxnSorKrRscZ32d+5VIpHwuzQABYipRACAkpEMQQu7nkuGoH2Z54r8k0gk1LyuSXXxfnWvDKihKTmtqudYt5pP/QPBZAAm4IoBAKBkpELQ3KRC0MpF16GDqov3q2tV9ZhpVYdXVWtRvF9dhw76XSKAAkNjAAAoGYSgPXG054haIwHXaVWtkYCO9r7tU2UAChWNAQCgZGQVglZZJWOMjDGqrDB6KmBUYZL/rawwo+eeq3/e729nVvoHBtVQ4z5juKEmoC8GBvNcEYBCxxoDAEDpyBSCNmee9Oih9F/+6cm2pt+aq4aagLOt6QNd+KZW91/9tQbeetnv72ZW6sO1uhEf0YpQ1YRzN+IjWhKu9aEqAIWMxgAAUBLSsw4e60n42f2vY9L53yZvjvRtTVNTbZLbmga07eKgPrp8XI/zW37ObW3bpp5j3Tq8aux0ImuteqMj2tra7mN1AAoRU4kAACUhm6yD1NqDctjWdOeu3bpbU68dV4fHTKvacXVYd2vqtXPXbr9LBFBgaAwAAGUnua3p5PPvhxNf5bmi3AsGgzpz7rzWbmnX9usVevH9O9p+vUJrt7SzVSkAV0wlyoIxJiSpVVLMOTRkrT3pW0EAgFlJbWs62fz76uCC5BSkIhcMBrWnc5/2dO7zuxQARYArBtlptdb+wlrba63tlRQxxkT8LgoAMDNsawoAE9EYZGfDuPtXJDX6UQgAYPay2tYUAMoMU4mys9AY02OtbXPuvyKpw8+CAACzkGlb02NbJWnC4uTxwkuW6nb/rXxUDACeK7rGwBizSdJma+0rk5wPSXpT0k3n0HJr7WzfxL8m6XfGmBZJPZJ6rLWxWT4nAMBPc+bpcdNPdL/pJ5KebG8qSYoPSDt/l/EpUtuj5lIikVDXoYM62nNE/QODqg/XamvbNu3ctZsFwwA8VTSNgTGmx/ljRNLCKR56QlKbtTbqfF3EGHPaWjt+OlDWrLVXjDG9klokHVDyisGVmT4fAMA/6XkHhSaRSKh5XZPq4v3qXhlQQ1NYN+Ij6jnWreZT/8BuQgA8VTRrDKy1bc5UnhOTPca5mhBNNQXO10XTzs2I05T0WGtXS2qTdGI2zwcA8E82eQfTlUgktL9zr5YtrtMPvv99LVtcp/2de5VIJKb1PF2HDqou3q+uVdVaEapSoMJoRahKh1dVa1G8X12HDk67NgDIVtFcMcjSZkmnXY6fVvIN/UlJMsa0Slqd4blOW2tPGmMaJd1MNRjW2l5jTFTJNQZsWQoAZS6Xn/If7Tmi7pUB1+C11khA23vfZutRAJ4ptcYgtQZgvKikl1J3nC1Hs7VQT/ILUl/fZ4xpc384AKBkfHNfFZePa+4nH2j43leqnr9AD779o+QC5TnzJI39lD/1hj75KX9AO64mP+XP9s18/8CgGprCrucaagL6YuBObr4vAHBRNFOJMnEWHYckDbmcjjnnZuKSxm1X6lxFOD7D5wMAFINv7mveez/V+v5TendNQNc2hvXumoC+9/kpzXvvp9I3yeXKR3uOqDUy+af8R3vfzvol68O1uhEfcT13Iz6iJeHamX8/AJBBKV0xmGpBsqRk8zDd3YSstTFjzM+NMQf0ZKcjT5KPHz9+rD/96U9KJBIaHh7W48ePc/0SZWXBggX6l3/5l6wfX1FRoerqagWDQT3zzDOqqCiZvhnADFRcPq6muYP6u+/MH3Ml4Mh3A9p2cVAfXT6ux8rtp/xb27ap51i3Dq8a22hYa9UbHdHW1vZZfU8AMJVSagw8Y629oix3IXLWL7RK0qJFi3T27NmsX2f+/PkKhUKqqanRggULVFFRkXEPbUzu0aNHqqyszOqx1lo9fvxYDx480NDQkKLRqO7du+dxhaUrkUhM63cfucX458bcTz5Q+5q5rlcC2l+YqwsXP9B9SXULF+hGfEQrQlUTnuNGfER1zyzI+ufxnTVr9X//t99q++VBvf7nc9VQE9CN+Ije/v8e6FZVrX66Zi0/2ynwu+8vxr/40RjkmLN+oVeSXnjhBbt+/fqsvu7OnTsaGRnR4sWLaQZy5Ouvv9bTTz897a+rq6vTv/3bvykQCKiurs6Dykrf2bNnle3vPnKP8c+N4XtfqaFm8isBw4kvpcoq3RmK6Vf/fa6O/A/PTPiUv/u/f62BoQf6n/7nLVkHof3V1f9XXYcOanvv2/pi4I6WhGu1tXWHTpJjkBG/+/5i/ItfKc2VcFtbMEYhh5J99dVXevbZZ2kKCoAxRs8++6y++uorv0sB4JPwkqWqrrBTzvevrpD06KEe/+d/0oXH/07bLt7T9dhDPXxsdT32UNsu3tOFx/9Oj//zP2ngi8+zfu1gMKg9nfv02R8HNPLokT7744D2dO6jKQDguZJpDJw3/TElA9DGi2jczkKFZmRkRHPmzPG7DDjmzJmjkRH3NwQASt/t/lv62/+6Rz3RkQm5Bqn5/n/7X//35IE583T/1V/ro6Ub9eOLI/r2qQH9+OKIPlq6Ufdf/fXo7kUAUOhKpjFwXJL7IuTlkvryXMu0cbWgcPCzALBz127dranXjqvDY64E7Lg6rLs19dq5a/eTB8+Zp8dNP9H91/5edkef7r/293rc9BOaAgBFpdQagxMat7Woo0VsLwoAmIZgMKgz585r7ZZ2bb9eoRffv6Pt1yu0dkv7tELLAKBYFOPi45AmySRwUonbjDGRVFKxkzngyfaiAIDSlprvT9owgHJQNI2BkyMQkvSqpJAx5oSSC457nO1EU34g6U1jTCpzYLm11u0qAgAAAABH0TQG1toO549tGR4Xk9Qx1WMAAAAAjFU0jQEgSdFoVB0dHerr61MsFlMkElFjY+Po+VgspqGhIb355pv64Q9/6GOlAOD4zY8lZd7UILxkadZZBwDgBRoDFJVIJKITJ06ora1Nvb29unz5skKh0JjHXLlyRatXr9bPfvYzHT582J9CASAlPiDt/F3Ghw10/SAPxQDA5GgMysxz9c9nFbRT6J9c9fX1qbGxcUJTIEmNjY1qaWnRL3/5SxoDAJ4KL1nKG3oAJYPGoMwMfPF50X9yFYvFFI1G9cYbb2R8bDQaVSTilnkHALPn9gHK2bNntX79+tH75KIAKBallmOAMtDXl8yq27Bh8s2mLl26pFAoRFMAAACQJa4YoOicPn1akvTSSy+5nu/o6FAsFtOxY8fyWRYAAEBRozFA0ZlqfUFfX596e3t14sQJdiUCAACYBhoDFJXU+oJIJKKOjo4J50KhkP7whz8oFArp66+/VjQa1YEDB7RhwwZt2rTJp6oBAAAKH40BikpqfcGBAwcyvtE/c+aMnnrqKUWjUQ0NDeWjPACYucqqrBYqF/qucQCKF40BikpqfUFLS0vGxzY3N+vpp59WT0+P12UBwOw9elj0u8YBKG40BigqfX19ikQirusLAKAQkXUAoFiwXSmKRmp9QTZXCwCgUNzuvyVrbcYbAPiNKwZlJttPrsJLluahmunJJr8AAAAAM0NjUGaKecHa8ePHJWW3vgAAAADTQ2OAgtfR0aFoNKqTJ09Kkl555RU1NjbqwIEDPlcGAABQOmgMUPBoAAAAALzH4mMAAAAAXDFA6frnf/5nnT9/Xn19fYpGo4rFYmptbWWrUwDFjSA0AB6hMUDJ+su//Ev91V/9ld544w2/SwGAjKaVd0AQGgAP0BgAAFAAsv10P5urBQAwE6wxAAAAAEBjAAAAAIDGAAAAAIBoDAAAAACIxgAAAACAaAwAACg9v/mxpOQORlPdnqt/3udCARQStisFAKDUxAfIOgAwbTQGAAAUkWkFoQHANNAYAABQRLIJQiMEDcBMsMYAAAAAAI0BAAAAAKYSochEo1F1dHSor69PsVhMkUhEjY2No+djsZiGhob05ptv6oc//KGPlQIAABQXGgMUlUgkohMnTqitrU29vb26fPmyQqHQmMdcuXJFq1ev1s9+9jMdPnzYn0IBAACKDFOJUJT6+vrU2Ng4oSmQpMbGRrW0tOiXv/xl3usCgKJSWZUx64C8A6B8cMWgTCUSCXUdOqijPUfUPzCo+nCttrZt085duxUMBv0ub0qxWEzRaFRvvPFGxsdGo1FFIpE8VAUARejRQ/IOAIyiMShDiURCzeuaVBfvV/fKgBqawroRH1HPsW41n/oHnTl3vqCbg76+PknShg0bJn3MpUuXFAqFaAoAlCWyDgDMBFOJylDXoYOqi/era1W1VoSqFKgwWhGq0uFV1VoU71fXoYN+lzil06dPS5Jeeukl1/MdHR2KxWLq7u7OZ1kAUDBu99+StTbjDQDS0RiUoaM9R9QaCUwIwDHGqDUS0NHet32qLDtTrS/o6+tTb2+vTpw4oY0bN+a9NgAAgGLFVKIy1D8wqIamsOu5hpqAvhi4k+eKspdaXxCJRNTR0THhXCgU0h/+8AeFQiF9/vnnevfddxWLxfT73/9ebW1tamlp8alyAACAwkZjUIbqw7W6ER/RilDVhHM34iNaEq71oarspNYXHDhwQJs2bZrysV1dXaPblcZiMT3zzDO6efMm6w4AAABcMJWoDG1t26ae6MiE+aXWWvVGR7S19XWfKssstb4g0yf/0WhUX3311ej9UCikTZs26cCBA57WBwAAUKxoDMrQzl27dbemXjuuDut67KEePra6HnuoHVeHdbemXjt37fa7xEn19fUpEom4ri8Y7x//8R8Vi8XGHBsaGvKmMAAoVb/5sSSRdQCUAaYSlaFgMKgz586r69BBbe99W18M3NGScK22trYXdI5Ban1Ba2trxsdGIhHdunVLTz/99OixK1euqK2tzcsSAaD0xAfIOgDKBI1BmQoGg9rTuU97Ovf5XUrWsskvmMyVK1ckKatQNAAoF+QdAEhHY4Cicfz4cUmZ1xe4ee2110bXJwAAkm7338r4mPFbWwMoXTQGKHgdHR2KRqM6efKkJOmVV15RY2Nj1guJOzo69M4777AbEQAAwBRoDFDwZrOTUG9vrzZv3qzGxkZJySlFqT8DAADgCXYlQsk6c+aMFi5cqEgkolgsplgspkuXLvldFgAAQEHiigFKUjQa1V//9V9PON7T0+NDNQAAAIWPxgAlKRKJKB6Pj9muFAAAAJNjKhEAAJi9yqqMIWgEoQGFjSsGAABgUtPKOiAIDShqNAYAAGBS2WQdSOQdAKWAqUQAAAAAaAwKibXW7xLg4GcBAADKDY1BgaisrNSjR4/8LgOOR48eqbKy0u8yAAAA8obGoEA89dRTSiQSfpcBRyKR0FNPPeV3GQAAAHlDY1AgampqNDQ0xFWDAvDo0SMNDQ2ppqbG71IAAADyhsagQDz99NOaP3++PvvsM8ViMY2MjDDPPY+stRoZGVEsFtNnn32m+fPnE44GALn2mx9LElkHQIFiu1KHMSYiqUPSaWvtSZfzb0iKSlooKWqt7cvx66uurk5ff/214vG47ty5w9WDWRoeHlZ1dXXWj6+srNRTTz2l2tpaPf3002y9BwC5Fh8g6wAoYDQGkowxLc4fI0q+8R9/vkfSiVQzYIw5bYy5ZK2N5bgO1dTUMIUlR86ePatVq1b5XQYAlIVpBaEBKEhMJZJkre1z3vTHJnnIq+OuEJyW9KrnhQEAUCRu99+StXbKG4DCRmOQgTGmUdLQuMMxSRvyXw0AAADgjYKbSmSM2SRps7X2lUnOhyS9Kemmc2i5tbbDw5IimnglYcg5DgAAAJSEgmkMnHn80iTz/NOckNRmrY06Xxcxxpy21ub7E/xQnl8PAAAA8EzBTCWy1rZZa9uUfOPvyrmaEE01Bc7XRdPOeSU07v5CTb4eAQAAACg6BXPFIEublVz4O95pSW2STkqSMaZV0uoMz+W6LamLmCZewQgpuXUpAADItcqqrLaMDi9Zqtv9t/JQEFAeiq0xaJHU43I8Kuml1B1rbW+uXtBa2+fyl9Ozcm9QAADAbD16SN4B4IOiaQycRcchTdwhSEp+qh/y8OX7jDGN1torzv1GST/38PUAACg5ZB0Aha1oGgNNvSBZUrJ5mEnomLMlaYtzizhNSG/ac70m6U1jzEIlm4IDuQ43AwCg1GU77YfkecAfxdQYeMa5EnBF0i8mOR+TlNoStc/tMSnO+oZW5+4DY8y1HJWJ6auVNOh3EWWKsfcX4+8vxn/2Vuu/bcvqgcaYy2l3GXt/Mf7+emG2T0BjkGPO+oZeSTLGXLLWvpThS+ARxt8/jL2/GH9/Mf7+Yez9xfj7yxhzabbPUTDblWbBbW3BGEzvAQAAAGamaBoD501/TO6Jw27pxAAAAACyVDSNgeOS3BchL1eGuf8+ydm2qZgRxt8/jL2/GH9/Mf7+Yez9xfj7a9bjb6y1uSgkZ5zFu23W2gkBZc65DdbaV8Ydvyzp51kGlgEAAAAYpxCvGIQ0SSaBs7A3YowZnU7kbDU6RFMAAACAUpP+vtdrBbMrkTHmgJINwauSQsaYE0ouOO5JCxaTpB8omSlw07m/3Fq7Ia/FTsHJQHhTUnp9HZN/BWbDGLNJ0ubxV5HSzofEzyPn0rI/nlUy2yMqqWP8BgCMvzecfyTa0g5FlBz/6LjHhcT4e865mh211vaNOx4S459Tzr/9ByS95xxqkbRB4/7+Yey9kza2Xyr5b4CUnLURc3kM458bbc6/uz1KrqmdsCFP6r3yrMfeWssthzdJpyVF0u5HJJ32u65Suzn/c/Q4432Zn0dex75RUuu4Y29Isuljzfh7Nv4RSW+MO7ZJ0p8Yf19+HiHnd7/F5Rzjn/vxtuNuN8f/3jP2no5/xGVsNykZ/Mr4ezfuJ1x+99NvPbka+0KcSlS0nE+vozbtU7vUn51zyBFrbZu1tk3J/1lc8fPwTItNTusbZa39hZIhgT2pY4y/Z9rGH7DJqZSh9HOMf968Kpdd8Rh/z3RIWq3kVYLV1trlduKVMsbeOyc08erk5vQHMP6eiFprjdtNyZ9Hm5SbsacxyK3Nki67HD8tl3/M4Tl+Ht5oM8a0uBzvU/Kyfgrj753NLsdiGrs+i/H3mPP/wWQ74jH+HrHWXrHW9tmx04zTMfYeSL2xHD/u1tpX7NipKox/7t10O+j8HZT+85j12NMY5FaLknOtx4tKIgkw//h5eCebhVCMvwestR123K5tzpzSkMZeQWP8vRcZ/2l1GsbfP4y9N9qU3dbwjH+Ojb9Kn6bRjl3bNOuxL5jFx8Uu7R9mt4TmmCbZaQne4OfhHWvt8klOReR8csH45907knpT/0Aw/t4zxrRO9o814+89ZyHmQuu+4Dskxt4LL0nqcT6ljig5xt9R2sJjxj9/jDFvONN4U/dDysHYc8Ugd9yC18ZwfmjID34eeeSM5SZJP3cOMf4eM8Y0GmPecHZwO56aY+pg/D3k7Ao12ZUCifH30necKS0xa22fMeaAsytUCmPvnZCSjcCQtbbXWdv0c0mX08aU8c8DpzEe/3dQTsaeKwYAcuGApJOWPJG8ceb5XnHepHYYYxZOcbkZuTVhAT7ypif9KoG1tsMYc9MYM2G7WORO2j76kfQ1BtbamDGmT8ntMdmONH/etJNs0z5bXDEAMCvOZeWXvPpLClOz1kadqwUHjDFv+F1PqXM+rX4v4wPhiUne/J9U8sMJeO/3LscuS2p1OQ4PeB12RmOQO25zusaw48Kf4Cl+HnngXJbsUDJ4MB3jn3+9evLmiPH3QOoyfBZjx/jn100l81Ukxt4rqXF1m0I3pCfz1xl/77XJvUHLydjTGOSIM9gxue/WEpHLPtfwDj+PvHlH0ivj/7Jh/L1hjAkZY04780vH+9J5TITx98ybkjY489pHb0q+Kepw7rcw/t5wpgxNuRc7Y++NtHF1m8e+0OVxjL93NsmlQcvV2LPGILcuyf1/muXKbosv5BY/Dw85b4g60psCY0xj2vxTxj/3IkpuRze6A1SaZ53/pj41YvxzbNxe7aOcKVwHxk1zYfxzL6aJv/dSckzTjzP23uhTcgzdpL9RZfw94ly1nOpN/qzHnisGuXVCyTTG8VokHc9zLeDn4RlnF5DjLnu4p++TzPjnmNN0/WKSRd4tkvrSGjXG31+Mf+65/Z0jJT9B7Um7z9h747jGhlimbFBynUcK4++d1NWAyaYNzXrsaQxyyNmlIpK+MMS55D/Ebi2eCWmSvXn5eXjDWWy83Plzo3NrcY6PBm8x/p75/bjtGVM/k4jSki0Z/7wLpd9h/D3R5/K7/4akaPouUYy9N5yxG0qfzuWMa2P61TTG31OpMY25nczF2Btr7WyLRBrnMs+behJfvXyyy8+YubR5va86/z2pZAfdk76VGj+P3DPGTPWXRm/6fvqMvzecv+g3K7mu4Fkl/7F4bfxaD8bfW87fQxE9mfPbp+SUoqhzPiTGP6fSfvel5N/9N9NDntIeFxJj7wnn936U27gy/t5wfv9PSFo92ULi2Y49jQEAAAAAphIBAAAAoDEAAAAAIBoDAAAAAKIxAAAAACAaAwAAAACiMQAAAAAgGgMAAAAAojEAAAAAIBoDAAAAAKIxAAAAACAaAwAAAACiMQAAAAAgGgMAAAAAkgJ+FwAAKF/GmEZJLZI2SGqz1kaNMW84p5dLWmitfcW3AgGgjNAYAAD8tNla22GMWS6pxxhzxVrbkTppjLlpjGm11vb6WCMAlAUaAwCAL5yrBb937kYkvSRp/NWBmJJXDgAAHmONAQDAN9bak84fX5L0c2ttbNxDGiXdzGtRAFCmaAwAAL6w1l6RJGNMRFJIUl/6eeeKgsYfBwB4g8YAAOC3FulJo5Bms6SotTaa/5IAoPzQGAAA/LZB7lcFWiWdlEavKgAAPERjAADwW4uk0+kHnGlEIUk9zqFNea4JAMoOjQEAwDeTrS+QtFBSzMk1iEhiOhEAeMxYa/2uAQBQpowxLZJ6rLUTtiQ1xpxW8kpCjBwDAPAejQEAAAAAphIBAAAAoDEAAAAAIBoDAAAAAKIxAAAAACAaAwAAAACiMQAAAAAgGgMAAAAAojEAAAAAIBoDAAAAAKIxAAAAACDp/wfYRgdRVtR13wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "plt.rcParams['text.usetex'] = True\n", "plt.rcParams.update({'font.size': 22})\n", "\n", "plt.figure(figsize=(12, 8))\n", "\n", "l1 = plt.semilogy(list(range(1,n+1)), P1, \"ks\", label=r\"$P_1$\", markersize=10, markerfacecolor=(0, 0.447, 0.741, 1))\n", "l2 = plt.semilogy(list(range(1,n+1)), P2, \"ko\", label=r\"$P_2$\", markersize=8, markerfacecolor=(0.85, 0.325, 0.098, 1))\n", "\n", "plt.legend(loc=\"lower left\")\n", "\n", "plt.xlabel(r\"$n$\")\n", "plt.ylabel(r\"$\\left| a_n \\right|$\")\n", "plt.xlim([0, 70])\n", "plt.ylim([1e-10, 1e1])\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Numerical instability: simple example\n", "\n", "What's going on? Due to round-off we have solved a different problem:\n", "\n", "P2 : $a_{n+1} = a_{n-1} - a_{n}$ with $a_0 = 1$ and $a_1 = \\phi + \\varepsilon$.\n", "\n", "The solution is\n", "$$ a_n = \\left(1 + \\frac{\\varepsilon}{\\sqrt{5}}\\right)\\,\\phi^n - \\frac{\\varepsilon}{\\sqrt{5}}\\, \\tilde{\\phi}^n,$$\n", "where $\\tilde{\\phi} = (-\\sqrt{5}-1)/2$. \n", "Notice that $\\left|\\tilde{\\phi}\\right| > 1$. This is an unstable iteration." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Slideshow", "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.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }