{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "****Floating Point Error in Horner's Rule Polynomial Evaluation****\n",
        "\n",
        "The following example is taken from *Applied Numerical Linear Algebra* by James Demmel, SIAM 1997.\n",
        "\n",
        "Here are three methods for evaluating the function\n",
        "\n",
        "$$f(x)=(x-2)^9 = -512 + 2304x -4608x^2 +5476x^3 -4032x^4 +2016x^5-672x^6+144x^7-18x^8+x^9$$\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEDCAYAAAA849PJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VGX2+PHPmUlCAoGQ0KSHDqFDaAKC9CZYQMFeWHRd\n3a/6E0XdVVRU7AXddbFhwYpgARWkCagIAULvECDUECAQICQz9/n9McOQMiEJmWQmyXm/Xnllbj8z\nmdxz73Ofe64YY1BKKVX22PwdgFJKKf/QBKCUUmWUJgCllCqjNAEopVQZpQlAKaXKKE0ASilVRgV8\nAhCRD0XkiIhs8NH6fhGREyIyO9v4BiLyl4jsEJGvRCTEF9tTSqlAFfAJAJgGDPLh+l4GbvEy/kXg\ndWNMY+A4cJcPt6mUUgEn4BOAMWYJcCzzOBFp5D6SXyUiS0WkeQHWtwA4lW19AvQBZrhHfQxcXbjI\nlVIqsAX5O4BLNBW4xxizXUS6AP/BtQO/VFWAE8YYh3s4EahdyBiVUiqglbgEICLhwOXAN64DdwDK\nuaddCzzjZbH9xpiBxROhUkqVDCUuAeBqtjphjGmXfYIxZiYw8xLWmQxUFpEg91lAHWB/4cJUSqnA\nFvDXALIzxpwEdovIKHC134tI20Ku0wCLgJHuUbcB3xcqUKWUCnAS6NVAReQLoDdQFTgMPAUsBP4L\n1ASCgS+NMd6afrytbynQHAjHdeR/lzFmrog0BL4EooA1wM3GmHO+fTdKKRU4Aj4BKKWUKholrglI\nKaWUbwT0ReCqVaua6Ohof4ehlFIlxqpVq44aY6rlZ96ATgDR0dHExcX5OwyllCoxRGRPfuf1SRNQ\nXvV6RKS3iKSISLz750lfbFcppdSl89UZwDTgbeCTi8yz1BgzzEfbU0opVUg+OQPwVq9HKaVUYCvO\nawDdRGQtcAB42Biz0dtMIjIOGAdQr169HNMzMjJITEwkLS2tKGNVqtBCQ0OpU6cOwcHB/g5FKa+K\nKwGsBuobY1JFZAjwHdDE24zGmKm4ir0RGxub4yaFxMREKlasSHR0NJlqASkVUIwxJCcnk5iYSIMG\nDfwdjlJeFct9AMaYk8aYVPfrn4BgEal6KetKS0ujSpUquvNXAU1EqFKlip6pqoBWLAlARC5z19xH\nRDq7t5tciPX5KjSliox+T1Wg80kTUOZ6PSKSiKteTzCAMeZdXEXW/i4iDuAsMNpoDQqllMph+a5k\nqoaXo3H18CLflq96AY0xxtQ0xgQbY+oYYz4wxrzr3vljjHnbGNPSGNPWGNPVGPOHL7brL3a7nXbt\n2nl+Jk+efEnr6d279yXf6LZ48WL++OPCx/juu+/yyScX64Wb//UOG3bpvXWjo6M5evSoz9Z3qW6/\n/XYWL17s8/Vmf39K+droqcvp99pvxbKtgL4TOFCFhYURHx/v1xgWL15MeHg4l19+OQD33HOPX+Px\nFYfDQVBQ4b6WTqfTr9tXqqTQYnA+8ssvvzBq1CjPcOYj37///e/ExsbSsmVLnnrqKa/Lh4dfON2b\nMWMGt99+OwA//vgjXbp0oX379vTr14/Dhw+TkJDAu+++y+uvv067du1YunQpEydO5JVXXgEgPj6e\nrl270qZNG6655hqOHz8OuM44Hn30UTp37kzTpk1ZunSp11hOnjzJ0KFDadasGffccw+WZfHhhx/y\nwAMPeOZ57733ePDBBwv0GR07doyrr76aNm3a0LVrV9atWwfAxIkTueWWW+jevTu33HIL06ZN49pr\nr2XQoEE0adKERx55xLOOefPm0a1bNzp06MCoUaNITU0FXEfmjz76KB06dOCbb74hIiKCkJAQACZM\nmEBMTAxt2rTh4YcfzhGXt+3fd999nunDhg3zejbx2Wef0blzZ9q1a8fdd99dqMSjlD+U6EOdp3/c\nyKYDJ326zphalXjqqpYXnefs2bO0a3fhgWSPPfYY1113HePGjeP06dNUqFCBr776itGjRwPw3HPP\nERUVhdPppG/fvqxbt442bdrkK54ePXqwfPlyRIT333+fl156iVdffZV77rmH8PBwzw5twYIFnmVu\nvfVWpkyZQq9evXjyySd5+umneeONNwDXEe6KFSv46aefePrpp5k/f36Oba5YsYJNmzZRv359Bg0a\nxMyZM7n++ut57rnnePnllwkODuajjz7if//7n9eYr7zySux2OwCpqak0b94cgKeeeor27dvz3Xff\nsXDhQm699VbPmdSmTZtYtmwZYWFhTJs2jfj4eNasWUO5cuVo1qwZ999/P2FhYUyaNIn58+dToUIF\nXnzxRV577TWefNJVWaRKlSqsXr0awPPZJycnM2vWLLZs2YKIcOLECa8xZ99+XjZv3sxXX33F77//\nTnBwMPfeey/Tp0/n1ltvzXNZpQJFiU4A/pJbE9CgQYP48ccfGTlyJHPmzOGll14C4Ouvv2bq1Kk4\nHA4OHjzIpk2b8p0AEhMTueGGGzh48CDp6el59ilPSUnhxIkT9OrVC4Dbbrsty5nJtddeC0DHjh1J\nSEjwuo7OnTvTsGFDAMaMGcOyZcsYOXIkffr0Yfbs2bRo0YKMjAxat27tdflFixZRtaqrl+/ixYs9\nZybLli3j22+/BaBPnz4kJydz8qQrgQ8fPpywsDDPOvr27UtERAQAMTEx7NmzhxMnTrBp0ya6d+8O\nQHp6Ot26dfMsc8MNN+SIJSIigtDQUO666y6GDRuW6/WI7NvPy4IFC1i1ahWdOnUCXAcF1atXz/fy\nSgWCEp0A8jpSL26jR4/m7bffJioqitjYWCpWrMju3bt55ZVXWLlyJZGRkdx+++1e+4Zn7jKYefr9\n99/PQw89xPDhw1m8eDETJ04sVIzlypUDXBeyHQ6H13myd188Pzx27Fief/55mjdvzh133FGoOLKr\nUKGC1zgzx2qMoX///nzxxRf5WgdAUFAQK1asYMGCBcyYMYO3336bhQsXXnTZoKAgLMvyDHv7exlj\nuO2223jhhRfyfnNK5dPBlLPcaF9Ae9kOZggUcVdivQbgQ7169WL16tW89957niaIkydPUqFCBSIi\nIjh8+DA///yz12Vr1KjB5s2bsSyLWbNmecanpKRQu3ZtAD7++GPP+IoVK3Lq1Kkc64mIiCAyMtLT\nvv/pp596zgbya8WKFezevRvLsvjqq6/o0aMHAF26dGHfvn18/vnnjBkzpkDrBOjZsyfTp08HXGcG\nVatWpVKlSvlevmvXrvz+++/s2LEDgNOnT7Nt27aLLpOamkpKSgpDhgzh9ddfZ+3atXluJzo6mvj4\neCzLYt++faxYsSLHPH379mXGjBkcOXIEcF3f2LMn31V4lfKq2wsL6W2Lp4Nte5Hv/KGEnwH4S/Zr\nAIMGDWLy5MnY7XaGDRvGtGnTPDvrtm3b0r59e5o3b07dunU9zRfZTZ48mWHDhlGtWjViY2M9Fzcn\nTpzIqFGjiIyMpE+fPuzevRuAq666ipEjR/L9998zZcqULOv6+OOPueeeezhz5gwNGzbko48+KtD7\n69SpE/fddx87duzgyiuv5JprrvFMu/7664mPjycyMrJA6zz/Xu68807atGlD+fLlsyS0/KhWrRrT\npk1jzJgxnDvnelzzpEmTaNq0aa7LnDp1ihEjRpCWloYxhtdeey3P7XTv3p0GDRoQExNDixYt6NCh\nQ455YmJimDRpEgMGDMCyLIKDg3nnnXeoX79+gd6TUtnF2PawxmpMo2LYVkA/Ezg2NtZk7ye/efNm\nWrRo4aeI1LBhw3jwwQfp27evv0MpEfT7qgqi3YQviQ+9m8kZo5nwnPdOFnkRkVXGmNj8zKtNQCpf\nTpw4QdOmTQkLC9Odv1JFIOVsBi1sewHYaKKLZZvaBKTypXLlynm2tyulLl3fVxdztSQAsNGKLpZt\n6hmAUkoFgKOp6bS0JXDQRHGM/HeOKAxNAEopFSBaSgIbrfpc3a5WsWxPE4BSSvmZMYZQztFIDrDR\nRBfbPU6aAJRSys9+XHeQ5rIPuxg2WdFEVggplu1qArgEpbkc9IEDBxg5cmSBllm6dCktW7akXbt2\nnD17tkDLfvfdd2zatKlAyxREbuWbv/nmG1q0aMGVV15ZZNtWKr92J52mpS0BgI2m+O4l0V5Al6A0\nl4OuVasWM2bMKNAy06dP57HHHuPmm28u8Pa+++47hg0bRkxMTL6X8UXJ5g8++ID33nvPc5dzcWxT\nqdyEBNloKQmkmPIkmmrFtl09A/CR0lIOOiEhgVatWgFctCzzee+//z5ff/01//73v7nppptITU2l\nb9++dOjQgdatW/P999975v3kk09o06YNbdu25ZZbbuGPP/7ghx9+YPz48bRr146dO3deNPYHHniA\n2NhY3nzzTZKSkrjuuuvo1KkTnTp14vfffwdc1T8HDBhAy5YtGTt2LN5udHzmmWdYtmwZd911F+PH\njyctLY077riD1q1b0759exYtWuR5/8OHD6dPnz707duXxYsX06tXL0aMGEHDhg2ZMGEC06dPp3Pn\nzrRu3ZqdO3d6/dsqdTEHU87y4i9baGXb7e7+WXyPEi3ZhzQ/T4BD6327zstaw+CLN+mU9nLQmXkr\ny1y3bl3P9LFjx7Js2TKGDRvGyJEjcTgczJo1i0qVKnH06FG6du3K8OHD2bRpE5MmTeKPP/6gatWq\nHDt2jKioKIYPH+5ZFqBNmza5xp6enu5pMrvxxht58MEH6dGjB3v37mXgwIFs3ryZp59+mh49evDk\nk08yZ84cPvjggxzv6cknn2ThwoW88sorxMbG8uqrryIirF+/ni1btjBgwADPPQ+rV69m3bp1REVF\nsXjxYtauXcvmzZuJioqiYcOGjB07lhUrVvDmm28yZcoUT6xK5deN7/1FOdJpIXt53xpSrNsu2QnA\nT0p7OejMvJVlzpwAsjPG8Pjjj7NkyRJsNhv79+/n8OHDLFy4kFGjRnnKREdFRRU49szlnufPn5/l\n2sHJkydJTU1lyZIlzJw5E4ChQ4fmq2bRsmXLuP/++wFo3rw59evX9ySA/v37Z4m1U6dO1KxZE4BG\njRoxYMAAAFq3bu05c1CqIHYfPU172UOwOIm3GrP9ucHFtu2SnQDyOFIvbqWlHLS3+fO7zPTp00lK\nSmLVqlUEBwcTHR3t9f1eiswlmy3LYvny5YSGhvpk3fnZJmT9PGw2m2fYZrPl6/NUypt2NleF23ir\nEcH24muZ12sAPlRaykEXRkpKCtWrVyc4OJhFixZ5SiT36dOHb775huTkZMBVPhmyvo+CxD5gwIAs\nVVDPn5FdccUVfP755wD8/PPPnmsIF5O5TPW2bdvYu3cvzZo1K/B7V6qgzl+jamvbyUETxWFynhkX\nJU0Al+D8NYDzPxMmTADwlIP++eefPReAM5eDvvHGG/MsB3355Zd7mhjgQjnojh07eppPwFUOetas\nWZ6LwJl9/PHHjB8/njZt2hAfH+95ZGJxuOmmm4iLi6N169Z88sknnsdBtmzZkieeeIJevXrRtm1b\nHnroIcB11vTyyy/Tvn17du7cme/Y33rrLeLi4mjTpg0xMTG8++67gOuxk0uWLKFly5bMnDmTevXq\n5Rnzvffei2VZtG7dmhtuuIFp06ZlOdJXqqjMWJUIQFvZyVqrOApAZ6XloJUqQvp9VRcTPWEOlTlF\nfOjdvJgxmv86h5MweWih1qnloJVSqoRoa9sFQLwp/jMATQBKKeUHGU7Xc6fbyQ4sI6y3GrB10qBi\njaFEJoBAbrZS6jz9nqqLmfzzFsDVA2iHqUUq5SkXZC/WGHySAETkQxE5IiIbcpkuIvKWiOwQkXUi\nkvMhq/kUGhpKcnKy/nOpgGaMITk5uci7qaqS67dtSQgWHW3bWWXl/lzrouSr+wCmAW8DuVUjGww0\ncf90Af7r/l1gderUITExkaSkpEtZXKliExoaSp06dfwdhgpQR06m0UwSqSRnWGk1Y9a9lxd7DD5J\nAMaYJSISfZFZRgCfGNdh+3IRqSwiNY0xBwu6reDg4DzvhlVKqUB3Ms3BCLurGWiFac5r9fK+a93X\niusaQG1gX6bhRPe4HERknIjEiUicHuUrpUqzTratHDRRxVoBNLOAuwhsjJlqjIk1xsRWq+afD0Up\npYrSp38mAIZOtq2stJpRnBVAMyuuBLAfyFxBrI57nFJKlTmfr9hHHTlKTTnmTgD+UVwJ4AfgVndv\noK5AyqW0/yulVGlgjKGTuNr/V1rN/RaHTy4Ci8gXQG+gqogkAk8BwQDGmHeBn4AhwA7gDHCHL7ar\nlFIl0YkzGXSybSHFlGer8V9PMV/1AhqTx3QD/MMX21JKqZJsXeIJDp1Mo0vIFuKsZhhs9G6mF4GV\nUqrUG/7279TiKI1sB/ndcj1+9e0bL/ne2EIp2Q+EUUqpEqi73VU0YZnVqtDVPwtDzwCUUqqYdbdt\nIMlEsM2P7f+gCUAppYrNVVOWAYbuto38brXEX/3/z9MEoJRSxWT9/hSaSiLVJMXT/u9PmgCUUqoY\n9bC52v9/d2oCUEqpUm/fsTNET5gDuNr/d1mXcYCqeSxV9DQBKKVUEXtrwXYAypHO5baNLLHaALBk\n/JX+DEsTgFJKFbVDJ9MA6GbbSJiks9Bqzwe3xVKvSnm/xqUJQCmlitjS7UcB6Gtbw2lTjr+sFjSu\nHu7nqDQBKKVUMTH0sa9hmdWac4RQq3KYvwPSBKCUUsWhueyjtiSzwGoPQLDd/7tf/0eglFJlQB/b\nagAWOdvRomYlP0fjoglAKaWKUMLR0wD0s69mvRVNEpFUCLH7OSoXTQBKKVWEer+ymJok08G2g1+c\nnQF4ekRLP0floglAKaWKyPmbvwbbVwDwk9UFgJa1IvwWU2aaAJRSqgicTMvwvB5s/4vNVj12m5os\nfcS/N39lpglAKaWKQJuJ8wCowTE62bYxx+k6+q8b5d+bvzLTB8IopZQPnU138vai7Z7h7M0/gUTP\nAJRSyof+s3gH7yza6Rkeal/OFqsuu0wtP0blnSYApZTyobPpTs/r+nKITrZtfO/s7seIcqcJQCml\nfMgyF15fa1+KZYRZ7gRwQ2xdP0XlnSYApZTyIYMrAwgW19mXssxqxSGqADDpGv8/BCYzvQislFI+\nsGF/CmsTT/DR7wkAdLFtoY4c5SXnDQBc3qhKQNT/yUwTgFJK+cCwKcuyDF9nW8IpE8Y8KxaA6WO1\nF5BSSpV6lUhlmH05Pzq7kUY5Vj7RDxHxd1g5+CQBiMggEdkqIjtEZIKX6beLSJKIxLt/xvpiu0op\nFYhG2ZcQJul86uwPQEiANf2cV+gmIBGxA+8A/YFEYKWI/GCM2ZRt1q+MMfcVdntKKRVozjkudP0U\nLG62/8pKqymbTX0AIsoH+yu0i/JFWuoM7DDG7DLGpANfAiN8sF6llCoRhr51of2/p209DWyH+dTR\n348R5Y8vEkBtYF+m4UT3uOyuE5F1IjJDRHLtDCsi40QkTkTikpKSfBCeUkoVrR1HUj2vb7XPI8lU\n4hersx8jyp/iapj6EYg2xrQBfgU+zm1GY8xUY0ysMSa2WrVqxRSeUkpdmpQzF6p+NpFE+tnXMN3Z\nj3RczT5/TOjjr9Dy5IsEsB/IfERfxz3OwxiTbIw55x58H+jog+0qpZTfzVid6Hl9T9CPnDHlmOYY\n6BkXCA9/z40vEsBKoImINBCREGA08EPmGUSkZqbB4cBmH2xXKaX86pzDybOzXf1dapPEcNsffOHs\nwwkq+jmy/Cl0LyBjjENE7gPmAnbgQ2PMRhF5BogzxvwA/FNEhgMO4Bhwe2G3q5RS/tbsX794Xo8N\n+gmA9x1D/BVOgfnkTmBjzE/AT9nGPZnp9WPAY77YllJKBYJTmZ74VZNkbrQvZKazJwfddX9KgsC8\nO0EppQLYgRNnae1+4hfA/wV9CxjeclyTZb5P7wrsnkCaAJRSqoAun7zQ87qR7GeU/TemO/uxnws9\nF8uH2OnZJLB7MmoxOKWUKoDk1HNZhh8O+pqzlOMdx4X7Xydd3YorAnznD5oAlFIq386kO+g4ab5n\nuLttPYPtK3k1YyTJRAAQWz+Sm7vW91eIBaJNQEoplQ9/7Uom5sm5nuFgHDwd9DF7rOpMdQ7zjB/S\nuqa3xQOSngEopVQ+rEw4lmX4DvvPNLYd4I708Zwj5ML47tHFHNml0wSglFIX4XBazNt0mB/XHvSM\ni5aDPBj0Lb86O7LIap9l/kCs+58bTQBKKXURjZ/4OcuwDYtXg98lnSCeyLgzy7QqFUIoSTQBKKVU\nAYy1z6GjbTv/l34vR4jMMu0fVzb2U1SXRhOAUkplcyotgwohQdhsWZtz2st2Hg76mp+dnfje6u4Z\n37h6OPMf6lXcYRaaJgCllHLLcFrMWrOfR2aso0XNSmw+eNIzLZKTvB3yFodMFI9m/A24kBy++0d3\nL2sLfJoAlFLK7a0F25mycAdAlp2/HSdvBr9DVU5ybcZEThLumXZDbF3Cy5XMXWnJjFoppQphwebD\nvPvbTr4a143Dp9I4fc5J4+rh7D9x1svchklBH3KFfT3jM8ax0TTwTNnx3GDstpLT6yc7TQBKqTLn\n79NXk+6wOOew6PaCq67PP65sxLrElBzz3m+fxZigRbzluJpvnL0942/qUo8ge8m+l1YTgFKq7DGu\nXze+v9wz6p1FO3PMNtY+h/8XPINvnT15zTEqy7Snh7cs0hCLgyYApVSZYlmGdKcFwJq9J3Kd7277\njzwW/AWznV1zXPQFSvzRP2gCUEqVMd1fXHjR6YLFo0FfcU/Qj3zvvJyHMv6OE3uWebY8O6goQyw2\nJT+FKaXKFKdlMMZkGZfusLj+f3+yeu/xLOO3Hz7FjiOnPMMZTouDKWm5rjuUc/wn+E3uCfqRTx39\neDDj3hw7/90vDCE02J7LGkoWPQNQSpUYxhgaPf4TI9rV4s3Rrho8GU6LjQdSWLH7GI99u57/3tyB\np37YyL5jZ0hIPgNAwuShRE+Yc9F1N5V9TAmeQmPZz8SMW5nmHEj2Zh8oWbV+8iLZM2kgiY2NNXFx\ncf4OQ6mAtuNIKmv2HmdUbN1i3u4pIsuHUCW8XLFt8+5P45i78TAAu54fgs0mjHr3D1YmHM9jydwJ\nFjfb5/NE0HROUZ4HM+5lmdU6x3wz770cyzLERkdd8raKg4isMsbE5mdebQIqo2avO8CfO5OzjDuZ\n6SHXgeTzv/ay8YCre166wyItw+l5/cSs9SSdOnexxUuFI6fSPO87u8FvLmH8jHXFHBH0e20JvV9Z\nzJ87kzlyKvdmldx8E7eP6Alz+PiPBCzLdSC679gZ/thxNMe8s9cdYPfR056dP0DDx39i1prEQu38\nm8tevg2ZyLPB0/jLasHgc5O97vy3TRpMh3qRAb/zLyhtAipjTqZlkJrm4L7P1wCuU2OATQdOMuSt\npbw5uh0j2tW+6DrSHRbBdvF6Kny+ffZiPSSOnEqjesXQHOPnbjxEusPiqra1sox/fNZ6T6wD31jC\n7qOnAWhaI5xth1M5meZgypj2Oda3aMsROjeIokI+79J8bd5WOtSPpGHVcLYdPkW/mBr5Wu67Nftp\nWqMiMbUq8fLcLTichseGtMjXsrn5csVemtQIp22dygTZbXR+bgFdG0bx5bhuOebNcLp2nlsOnaRG\nxVAyLCvL5/vOoh30bVEdQZi/+TDjrmjI2n0nPDuzCd+6yh70bFKVFbuPcVlEKKv3HOfB/k2z/I1P\nn3Mw9K2lxNSqxBNDYwA4leZgzHuurpRvjWlPt4ZVqFbRdUYw8PUlHD6VRvyTAzjncNL1+QX8a2gM\nczce4s+dyZw65wDgqR82EhEWzNXta9PzpUUAPH9Naw6cOMvV7Wvx0/pDvPbrNq+f04Nfrb2kz/cy\nkrk/6DtusC/iBOE8mP53Zlk98Nbk07F+JCFBpfNYWRNAKTJzdSIPfb2WL8d1JbJ8CM0uq0iHZ3/l\nzu7RNKlRkTZ1Ijw3vWRmWYa5Gw8B8Nu2pBwJIOVsBoPfWMJ9fZowol0tWj41l7/1bODZCWR2xUuL\nSEo9x7ZJgzlw4izHz6TTslYEJ86k0+6ZXxnTuR5frNjLa9e35doOdbIse/enqwAIC7Z73fku3Z7k\n2fkDbDucCrjahds9M49xVzTk+Ol03lu6m+lju3DHtJW0rh3ByI51aFmrEk1qVCQiLDjXz+8tdwmA\nYLuQ4TQkTB7Kv75bz2fL93Jth9rE1KxE3xY1aFC1AuBqex761lJPHJufGeTpS35V21ps2J/ChJnr\nuT62Di+NbAu4Hin4t0/ieHZEK6IqhCAI+46fYfmuZH7ecIjE42eY90AvJsxc74mraQ1X2YHlu1wP\nJDnncHL3p6toXC2cHk2qeuYb9MZSz+tbutbn2atbse/YGV6eu5WX5271TFuXeIK5Gw9zT69GvPtb\nzr7v5y3ceoSoCuX45M7OAMxYlUhCsqtdfVfS6Rzz//ML10FFlwZR/LX7wsNTPlu+hw71Ijl+JoP/\n9433HfYDX8Xzx84LR/7nk/7bi3bkGt+lqC+HuN0+lxvtCxEspjv78ppjFCmZSjtkNuefPWhZK8Kn\nMQQSvQYQwM45nATbbDkqEp4XPWEO/VrUYOotHVmRcIzRU5dnmd6kejjbj6RedBtLxl/JFS8v8gyX\nD7Hzx4Q+VC7vqmu+/8RZrpqyjGOn03Msu/uFIVmOED9dvod/f7cBgA1PD6TVU67H5yVMHsqG/SkM\nm7Isy/If3h5L6jkn0VXKE1k+xHP0l9nV7WrxXfwBAOpXKc8e90W9zKqGh3A0NWd83jwyqBn39m6M\nZRn+2n2MCuXsvP7rNsJDg/lx7YEs87atE8FaL3eGfnJnZzpFR9H+2XmkZVheY+3bvDoLthzxTJt5\n7+VM+HYdt1/egMdnrad/TA1+3XQ4x7ov9j4BBrW8jG1HTnndARe16CrlPRdVC6ogfyNfC8ZBT9s6\nbrIv4EpbPE5szHT25C3HNezn4g9uX/5YXy6LyHm2GsgKcg1AE0AAOpWWwTM/buKbVYme3g6WZRAB\nY1yn3RHlgz29Gh4d1JwXf9nis+13jo7iuo612X44lfeX7b7ovN/+vRt/7T7G9bF1ic30sOzsIsKC\nSTkbGNe7v42dAAAc2klEQVQYvhrXlZmr9/NV3D5/h6KKSDnS6WLbzGDbCgbbV1BZTnPEVOZzZx+m\nO/qSlK2Of2Yf3d6JO6at5J99GvPQgGbFGLVvaAIowZJTz/HKvK18seLCzumlkW14ZMY62tSJoFmN\ninyzKtGPESoVeMI5QytbAh1kO91tG4i1baOcZJBqQplnxfKjsxvLrNZk5KPVO2HyUDKcFsEl9E7f\ngiQAn1wDEJFBwJuAHXjfGDM52/RywCdARyAZuMEYk+CLbZd07yzawctzt7Lx6YE4LENHL0fRj7h7\neKxLTPFarEqpsqISp4mWQ0TLIRrIIRraDtJKdtPIduF5vZutenzi7M8yqzXLrRZZHtieXYOqFfjl\ngZ4AvP7rds81kZK68y+oQp8BiIgd2Ab0BxKBlcAYY8ymTPPcC7QxxtwjIqOBa4wxN+S17rJwBnC+\nGWdkxzrM0CN7VYoJFiE4CMbh+R0m5wjnLBXlLBU5QzhnCZezVOQsVeQk1eQE1eUE1XD9DpcL3U0t\nIxygChutaNZZDdlgGrDeasAxKuU7pkUP9/Zc1C8tivsMoDOwwxizy73xL4ERwKZM84wAJrpfzwDe\nFhExgdz+VMQSjp7O0iPClzv/cqQTwWkqyhkqkOb5h6rAWcrLOcqRceGfUFz/iJnH2cWJDeP+sRDP\n6+zDFgJZxuVF8jVP3kR8s62LzZ95WHIdn/c8+Gg9Rb9+75/XpazHJibr9wwHweL9PobcpJpQjpjK\nJFGZjSaaxVZlDplIdpua7DaXsc9Uv+jR/cV0qFeZe3o1KnU7/4LyRQKoDWS+mpYIdMltHmOMQ0RS\ngCpAjjs+RGQcMA6gXr16PggvMPV+ZfElLVeBs9STI9STI9R1/64uJ6gqKVQlhaqSkuUoKS9OI6QT\nTAZBnCOIDIJwGDuWZ7cuWO5du+u3YHl+Mg+75jH52H3nZ578yGs9BsGY/MST+zpzm5bbtnObJz/r\nyT0O36yn6Nd/4bXTsrm/U67vVgZ20k3W71kGQaSZEE5RnlOmPKmEcYowUk0YqYSRTu5ddgvjkUHN\nuLN7g1JTz6cwAu4+AGPMVGAquJqA/ByOTyUcPc3vO4+yYf/JvGcGojhJrG0rMbY9xMgeYmx7qCNZ\nc2aKKc8hE8VRE8E6GnLUiuCoieAE4e5/pFBSTRin3f9cZ0wo6QS5f4Kx9GZwVYbc3ash9/Zu7O8w\nAoYvEsB+IHMRkjrucd7mSRSRICAC18XgUivlTAbly9lZuj2J8d+s4+M7O+foB59dOdLpbtvAFbZ1\ndLVtprnNdWLlNMIuU4tVVlM+t/qSYGqw11Rnr6me5dmkSqkLRnasQ9eGVXj4m7W0q1uZ3s2qMe6K\nhv4OK6D4IgGsBJqISANcO/rRwI3Z5vkBuA34ExgJLCzt7f9tn5mX5Y7I3Hb+QTjoa1vNVfbl9LbF\nEy5pnDbliLOa8UPG5Sy3WrDJ1CeN4iu4pVRJtvuFIWw7nEqzyyoCrkSgvCt0AnC36d8HzMXVDfRD\nY8xGEXkGiDPG/AB8AHwqIjuAY7iSRKmX+Xb47GpxlJuD5jPK/hvVJIUkE8H3zu7MtWL502qZr/7K\nvhYabCPDaXBapTo3F4vwckGkumvdqLwNiKnBvFzujM7Lx3d25rYPVwDwt54NEBHPzl9dnE8agI0x\nPxljmhpjGhljnnOPe9K988cYk2aMGWWMaWyM6Xy+x1Bpte9Y7rfLR8tBXgyaym/lHmScfTbxVmPu\nSB9Pl3Pv8ITjLpZYbQu983/35o45xq19ckCey62fOJCdzw/hhWtzVkMEV2mJyPLBfD62C7ueH3LJ\n8dWvUp6P7ujE7hfyt46EyUP56/G+vHerq2db4+r5b/bq3rhKluE3R7fzOt+oTEeJQ9vUzDF9QD4L\nw51X0KPOqAqX1pvFm0/v6lyo5dc+lfd3BaBuVJjX8c297HznPXgFVcO9v8eEyUMZ2PKyLONu6Vo/\ny/AbN+T8uz0yqBlfjetKr6bVXOVGnh7IY4MLV4SvrAm4i8AlncNpea1pU5lTPBg0g5vt88kgiM+c\n/ZjqGMZBqnhZi3ef3NmZelHlCQ6y8flfe3I8xPqVUW09O56q4SE4LcMt3aLpFB1JRPkLPSr+erwv\nNSqFsmF/Cuv3p/CYu/DY+ZtfRneqS2T5YHo3q05ahpMZqxKZNGczVzStxr+HXSgAt/mZQcTtOcb6\n/Sm89Iur2Nhjg5tzMCWNmhGhvPDzFk9cD3+zlh/u605MzUqICHYv9Y2mj+3CTe//Ra2IUN67LZZP\n/tjD5e4deI1KofSPCfVULz2ZlsE/v1jD4q1JOdYzfmAzmtWoyNhP4gjL1tNjRLva/N+X8VnG/bNP\nY+7v24Q7ujfgbIaDjvWjmLPuwsNDbu1W31OKuUqFEJIz1UVaP3EArSfOA1w7smXbj9KqdiUsAwdO\nnKVWZddOctOBk7xzUwcqhgbR/N+/eJa3Cex6wfWeZq87QMtaEVQOC6b9s7965unTvDqVQoM8dYYA\nOjeIYoWXM8wNTw8kPFv102dGtOTWbtEXPoN3fmftPu/Pwp19fw+vBfPGD2zGyoRjWT7vfw+NoU2d\nyvR4cSG1Koex133g07BaBa6Prcszsy/0BG9YtQIrHu9Hw8d/8ox76qoYhrsrv0Zn64757NWtuKlr\nPQa9sZSmNcK5un1tHvjqwt/NW9Xa7O9b5U0/MR+wLEPK2QwclmHoW0uzTTWMsS/k0aAvCecsnzn7\nMcVxLUfJvcLgb+N7c/8Xa7Lc9du5QRRXNL1QuGr8wOZ0blAFyxh2JZ2mXd0IOtS7UN8k7l/9c6z3\nqra1+HHtAU/3t1a1I2hVO8KTAM4TEQa1ch0FhwbbuaptLT5ctpubsx2VhYXY6dmkGj2bVKNLgyhP\n6eLz7u7VyPO6d7NqVPXy4JCf/tmTcw4ntSqHUaNSKIsf7k1khRAiwoJ5cWSbXD+jSqHBdGtYhcVb\nk7jvysaeqpG9mlbjH1e6ir092K8pN3WthzGwdt8Jr2cOvzzQk+aXuW4ciql14QaihtUqMKjlZdhE\nuLtXQya7k9mEwc1zPHhlSOvLqOZ+b5mrc069Nee9OA7nheJxyx69MkuCGtamVo75q1QI4T83dSA0\n2M5LI9ty7HS6pziZtydcedsJdm6QtYb99//oziMz1vJ1XNZ7T0KDbbSq7fpe9mxSlaXbXT3Ovv17\nNzrWv7COjs/+SvLpdCwDl0WEsuP5ITz0VbwnAXRtWIVbu0Xzddw+thw6xfiBzbyWB68YGux5mEzH\n+pHMf+gKwkKCPGcK9aLKA/BAv6aA67tiGeOJURWe1gLygX9/t4FPl+/JMb46x3k5+H/0sq/jT2cM\nTzluY5vJ+6lNCZOHknrOwbLtR7nnM1eJ5Ju61OO5a7w3zeRXWoaTfcfO0KRG1lP07+P3U7tyWIl7\n2IUxhgMpadSuHMbR1HMs236U3s2qeSqZ5mbcJ3Ge9ubzZxR5STmbwZQF23lkUPNC1YY3xtB64jwm\nDG6eI6FmdvtHK+jTvHqWI/fs3l64nVfmbWNM57o4nIYDKWeZPrYrABN/2MjljaowIFvTSmY7jqSS\nes5BapqDb1cnMqZzPU+yMMaw8cBJrzvb7+P3839fxme5i/aHtQc85aDPV4mdu/EQd3+6io1PD/Q8\nk2H5rmQiy4fwvyU7ef6a1toXvwhoMbhidOx0Oh0yna6fd4VtLW8Gv0Mo6TznuInPnP3I6x7XJ4fF\ncGePBp7hM+kOYp6cS9XwEBY81DtLM44qnPNHz/lNACorb8XSPl2+h6oVQhjcOuc1FFV8ir0YXFn2\nwk+bs40xjLPP5tGgL9lm6nBvxgPsNvn7h8i88wcoHxKkO6gi8u7NHTh8svQ/SrKoeCuWlv3CrQp8\nmgAKIfWcI0tpZhsWzwe9z+igxcx2dmF8xt2cJfeHSXw1riuRFUJIPH6GiDDf9QJReTt/jUOpskwT\nwCVauOUw9aIu9FwIxsHrwe8wzP4XUxxX86pjFJmbfERcvWY6PTefU2kOXrquDV0aunq4NK2hfZaV\nUsVPE8AlSD3n4M5pcdSJdHXxC8LBf4LfoL99NZMybuJ9Z85mm/gnBxAabGf9xIHFHa5SSnmllcAK\nKPH4Gc5fOE88fhbB4qXgqfS3r+ZfGXfk2PmLwH9v6nDRh5ErpZQ/6BlAASzccpg7p8Xx7IiWnnGP\nB33OtfZlvJxxPZ85c/a93z5psNc+0Eop5W+aAArg/F2Q//5+IwCj7Iv5W9BPfOQYyDvOEVnm3fn8\nECxjdOevlApYmgAKYMvBU57X7WQHk4I+ZKmzFZMcN5O9j7/dJth99OATpZQqCpoA8mniDxtZkeCq\nvVKZU/w35A0Om0juz7gfJxfuZpz/0BVsO5zqrzCVUirfNAHk07Q/EtyvDC8Ev08VUrgm41lOkLUL\nZ+PqFWlcXbt1KqUCnzZQF9Ao+28Mtq/kFcf1bDTRAMTUdBUR+3JcVz9GppRSBaNnAPnw1PcbAKhN\nEhODPuZPZ0yW7p6z7++BzUt5Y6WUCmR6BpCHhVsO8/GfrkqfE4M/QYCHM+7O8jB13fkrpUoiTQB5\nePFn14NO+tvi6G9fxeuO69iPqy5/y1qVWPrIlf4MTymlLpk2AV2EMYath09RnjQmBn/MZqsuHzkH\neabP+WdPP0anlFKFowngIvq99hsAf7PPobYk88/0+3C4P7LbL4/2Y2RKKVV42gSUi9V7j7Mz6TRV\nSWFc0Gx+cnZmlWnmmV5ZH86ilCrhNAHk4tr//AHAP4NmUo4MXnbckGV614b5f5i7UkoFIm0C8mLt\nvhMARMtBbrQv4HNnX3abmowf2Ix/XNmYtAynPstUKVXi6RmAFyPe+R2Af9i/J4Mg3nJcC8DljVxH\n/brzV0qVBpoAsjl/9F+bJK62/86Xzis5SgT1q5Snfb1IP0enlFK+owkgm/NH/3cHzcYA/3MM49Zu\n9Vn8cG+/xqWUUr5WqAQgIlEi8quIbHf/9nqILCJOEYl3//xQmG0Wh2oc5wb7Yr51XsEhqvDMiFaI\n6N2+SqnSpbBnABOABcaYJsAC97A3Z40x7dw/wwu5zSJ3a9CvBOPgXedV/g5FKaWKTGETwAjgY/fr\nj4GrC7k+v/rH9NWUI50x9oUssNqzx1zm75CUUqrIFDYB1DDGHHS/PgTUyGW+UBGJE5HlInLRJCEi\n49zzxiUlJRUyvPzbeugUc9YfZKhtOVXlJB87BxbbtpVSyh/yvA9AROYD3g6Fn8g8YIwxImJyWU19\nY8x+EWkILBSR9caYnd5mNMZMBaYCxMbG5rY+nxv4xhLAcFvQPHZYtVhmtQLg//o2Ka4QlFKqWOWZ\nAIwx/XKbJiKHRaSmMeagiNQEjuSyjv3u37tEZDHQHvCaAPypneykrW0X/8q4g/PP+H2wf1P/BqWU\nUkWksE1APwC3uV/fBnyffQYRiRSRcu7XVYHuwKZCbtenLMt1ojHavpDTphyznD0A+NfQFv4MSyml\nilRhE8BkoL+IbAf6uYcRkVgRed89TwsgTkTWAouAycaYgEoAL87dQhhpDLMvZ46zK6cJA+CuHg38\nHJlSShWdQtUCMsYkA329jI8Dxrpf/wG0Lsx2itr/ftvFtbYVhEsa3zh7AdC4erj2/VdKlWp6J7Db\nKPsSEqwarHSXfA4N1o9GKVW6lflqoGv2HqeOHKGbfROvZIwChHt7N+LGLvX8HZpSShWpMp8ArvnP\nH/zTvgzLCDOdrkc8PjKouZ+jUkqpolem2znOpjsBGGpfzkrTjANU5f4+jf0clVJKFY8ynQD+Pn0V\nTSSRZrZEZju7AnB/H73xSylVNpTpBLB4axJD7cuxjPCLszOvjmpLSFCZ/kiUUmVImd3bvTF/G2AY\navuLv6wWJFGZazvU9ndYSilVbMpwAthOU0mkiW0/c6wuANrvXylVppTZBAAw1P4XTnfzT7u6lf0d\njlJKFasy2Q009ZwDgL621awyTTlKBHH/6O7nqJRSqniVyTOA2z5cwWUk08qWwEJne3+Ho5RSflEm\nE8CqPcfpY48HYL7Vge3PDfZzREopVfzKZAIAV/PPXqsaO0xtgu1l9mNQSpVhZXLPF8o5uts2sMDq\nwPkHvyilVFlT5hKAMYbutg2ESgYLrA48pE/8UkqVUWUuAbR/9lf62OJJNaH8ZbWgZ5Oq/g5JKaX8\noswlgBNnMuhhW8+fVgwZBNG+XqS/Q1JKKb8ocwmgrhymvu0Iy6zW7H5hiL/DUUopvylTCWDm6kS6\n2zYCsMxqpaUflFJlWplKAM/M3kQP2wYOmUh2mlr+DkcppfyqTCWAlDPnuNy2gWVWa/rHXObvcJRS\nyq/KVAKIkb1ESSrLnK1479ZYf4ejlFJ+VWYSwO6jp+luWw/A71ZLP0ejlFL+V2YSwNLtSXS3bWSr\nVYcktOunUkqVmQSQfPIMHWzbWWE193coSikVEMpEAvhuzX7mL15IRTnLSk0ASikFFDIBiMgoEdko\nIpaI5HpVVUQGichWEdkhIhMKs81L8eHvu+ls2wLASqsZa/7dv7hDUEqpgFPYM4ANwLXAktxmEBE7\n8A4wGIgBxohITCG3WyC7j54m1raVRFOVg1QhskJIcW5eKaUCUqESgDFmszFmax6zdQZ2GGN2GWPS\ngS+BEYXZbkGdSsugs22rtv8rpVQmxXENoDawL9NwonucVyIyTkTiRCQuKSnJJwFEyyGqSQorrWY+\nWZ9SSpUGeT4UXkTmA95um33CGPO9rwMyxkwFpgLExsYaX6yzk811krLSaqbln5VSyi3PBGCM6VfI\nbewH6mYaruMeVyyclqGzbOGYCWeHqc0vt3cqrk0rpVRAK44moJVAExFpICIhwGjgh2LYLgAxT/5C\nB9t2VllNASFIn/+rlFJA4buBXiMiiUA3YI6IzHWPryUiPwEYYxzAfcBcYDPwtTFmY+HCzr9yjpM0\nsh1kjdWY+lXKF9dmlVIq4OXZBHQxxphZwCwv4w8AQzIN/wT8VJhtXao2tt0ArDWNeF8LwCmllEep\nbg85fc5BW9kJwHqrIU1qVPRzREopFThKdQL4Y2cy7Ww72WHV4iQV/B2OUkoFlFKdAP72yUra2Xaw\n1jTydyhKKRVwSnUCqEUy1SSFeKsR17TP9d4zpZQqk0ptAliz9zhtba72/3irMdUqlvNzREopFVhK\nbQJ4ZvYm2tp2cs4EscXUo25kmL9DUkqpgFJqE4BloL1tB5tMNBkEcVOX+v4OSSmlAkqpTQDr9h2j\npSSw1moIgM0mfo5IKaUCS6lNAPXlMOGSxkYTzf/r39Tf4SilVMAptQkgRvYAsMmqz/19m/g5GqWU\nCjylMgGkOyxa2hLIMHa2mzr+DkcppQJSqUwAr/26jRjZww5Tm3SC/R2OUkoFpFKZAN79bSctbXvY\naKLp16K6v8NRSqmAVCoTQDVOUF1OsNGqz3taAVQppbwqlQkgxnb+AnA0Itr9UymlvCmVCaClJACw\nyejNX0oplZtSmQBibHvYa1XjsWu6+DsUpZQKWKUzAUgCG000A1rW8HcoSikVsEpfAnBm4MTOeqsB\n5UPs/o5GKaUCVqGeCRyQ7MH0T38ZMDwSUvrenlJK+UrpOwPw0N4/Sil1MaU4ASillLqYUtlGMmVM\neyqFaQkIpZS6mFKZAK5qW8vfISilVMDTJiCllCqjNAEopVQZVagEICKjRGSjiFgikmvVNRFJEJH1\nIhIvInGF2aZSSinfKOw1gA3AtcD/8jHvlcaYo4XcnlJKKR8pVAIwxmwGtOKmUkqVQMV1DcAA80Rk\nlYiMu9iMIjJOROJEJC4pKamYwlNKqbInzzMAEZkPXOZl0hPGmO/zuZ0expj9IlId+FVEthhjlnib\n0RgzFZgKEBsba/K5fqWUUgWUZwIwxvQr7EaMMfvdv4+IyCygM+A1ASillCoeRX4jmIhUAGzGmFPu\n1wOAZ/Kz7KpVq46KyJ5L3HRVIBAvOmtcBaNxFYzGVTClMa58PwlLjLn0VhYRuQaYAlQDTgDxxpiB\nIlILeN8YM0REGgKz3IsEAZ8bY5675I3mP7Y4Y0zAPRBY4yoYjatgNK6CKetxFbYX0Cwu7Nwzjz8A\nDHG/3gW0Lcx2lFJK+Z7eCayUUmVUaU4AU/0dQC40roLRuApG4yqYMh1Xoa4BKKWUKrlK8xmAUkqp\ni9AEoJRSZVSJSwAi8qGIHBGRDblMjxSRWSKyTkRWiEgr9/i6IrJIRDa5K5j+XyDElWm6XUTWiMjs\nQIlLRCqLyAwR2SIim0WkW4DE9aD7b7hBRL4QkVAfxpXn90Rc3hKRHe74OmSadpuIbHf/3BYIcYlI\nOxH5073cOhG5IRDiyjS9kogkisjbgRKXiNQTkXnu7/0mEYkOkLheci+32T1P4QqxGWNK1A9wBdAB\n2JDL9JeBp9yvmwML3K9rAh3crysC24AYf8eVafpDwOfA7ED4vNzDHwNj3a9DgMr+jguoDewGwtzD\nXwO3+zCuPL8nuLo4/wwI0BX4yz0+Ctjl/h3pfh0ZAHE1BZq4X9cCDvrqb1mYuDJNf9P93X87EP6O\n7mmLgf7u1+FAeX/HBVwO/A7Y3T9/Ar0LE0+JOwMwrhpCxy4ySwyw0D3vFiBaRGoYYw4aY1a7x58C\nNuPamfg1LgARqQMMBd73VTyFjUtEInDtpD9wT0s3xpzwd1zuaUFAmIgEAeWBAz6MKz/fkxHAJ8Zl\nOVBZRGoCA4FfjTHHjDHHgV+BQf6OyxizzRiz3b3sAeAIrps3/RoXgIh0BGoA83wRjy/iEpEYIMgY\n86t7+VRjzBl/x4WrqGYoroOxckAwcLgw8ZS4BJAPa3E9owAR6Yzrtug6mWdwn861B/4KkLjeAB4B\nrGKMJ6+4GgBJwEfiapp6X1ylPPwal3HVlXoF2IvrSDbFGOPTncd5F/me1Ab2ZRpOdI/Lbby/48q8\nbGdcO5Cd/o5LRGzAq8DDvo6lMHHhOmM6ISIz3d/9l0XE7u+4jDF/Aotwfe8PAnONuyT/pSqNCWAy\nrowZD9wPrAGc5yeKSDjwLfCAMeakv+MSkWHAEWPMqmKMJc+4cB1ldwD+a4xpD5wGJvg7LhGJxHWE\n1ABXc0YFEbnZ1xv34/fkogoTl/so8lPgDmOMTw82LjGue4GfjDGJvozFB3EFAT1xJaZOQEPgdn/H\nJSKNgRa4DtBqA31EpGdh4ijyYnDFzf1h3gGuiym42ot3uYeDcX3o040xMwMkrhuA4SIyBNfpXSUR\n+cwY4/OdWgHjKg8kGmPOH53MoBgTwEXiGgjsNsYkuafNxNU2+pmvtp2P78l+oG6m4TrucfuB3tnG\nLw6AuBCRSsAcXGXcl/sqpkLG1Q3oKSL34mpnDxGRVGOMT75nhYgrCFdds/P7je9wtcV/4Oe4bgaW\nG2NS3ev5GddnuPRSYyl1ZwDi6rkS4h4cCywxxpx070Q+ADYbY14LlLiMMY8ZY+oYY6KB0cDC4tr5\n5xHXIWCfiDRzT+sLbPJ3XLiafrqKSHn337QvrnZUX203P9+TH4Bb3b01uuJqhjoIzAUGiKsHUySu\nyrdz/R2X+3OchatdeYYv4vFFXMaYm4wx9dzf/Yfd8flq51+Yv+NKXGef56+T9MFH3/1CxrUX6CUi\nQe4k0ovCfvfze7U4UH6AL3C1f2Xgahu7C7gHuMc9vRuuK+tbgZm4e2EAPXBdRFkHxLt/hvg7rmzr\n6I3vewFdclxAOyDO/Zl95y1mP8X1NLAF1zOpPwXK+TAur9+TbLEJ8A6udvT1QGym5e8Edrh/7giE\nuHAdOWZkWi4eaOfvuLKt53Z82wuosH/H/u5l1wPTgBB/x4Wr58//cO30NwGvFTYeLQWhlFJlVKlr\nAlJKKZU/mgCUUqqM0gSglFJllCYApZQqozQBKKVUGaUJQCmlyihNAEopVUb9f2laX3OOqy6jAAAA\nAElFTkSuQmCC\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7f5638e32828>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "import numpy as np\n",
        "\n",
        "#Evalue  polynomial in factored form\n",
        "def f(x):\n",
        "    return (x-2.)**9\n",
        "\n",
        "#coefficients for expanded form\n",
        "coeffs = np.asarray([-512., 2304., -4608., 5376., -4032., 2016., -672., 144., -18., 1.])\n",
        "\n",
        "#Evaluate polynomial using coefficients\n",
        "def p(x):\n",
        "    return np.inner(coeffs, np.asarray([x**i for i in range(10)]))\n",
        "\n",
        "#Evaluate Horner's rule for polynomial\n",
        "def h(x):\n",
        "    y = 0.\n",
        "    #[::-1] looks at all elements with stride -1, reversing the order\n",
        "    for c in coeffs[::-1]:\n",
        "        y = x*y+c\n",
        "    return y\n",
        "\n",
        "#Define 8000 points between 1.92 and 2.08\n",
        "xpts = 1.92+np.arange(8000.)/50000.\n",
        "\n",
        "import matplotlib\n",
        "import matplotlib.pyplot as pt\n",
        "\n",
        "#plot functions evaluated at each point using Horder's rule and using factored form\n",
        "pt.plot(xpts,[h(x) for x in xpts],label='Evaluation by Horner\\'s rule')\n",
        "pt.plot(xpts,[f(x) for x in xpts],label='Evaluation in factored form')\n",
        "pt.legend()\n",
        "pt.show()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "It seems Horner's rule is inaccurate when $f(x)\\approx 0$, lets try to understand why.\n",
        "\n",
        "The first method uses $(x-2)^9$ directly, encurring a backward error of $\\epsilon$ (machine epsilon), since given $z=\\textit{fl}(x+2)$, we can obtain $z^9$ to the same precision, therby solving the problem for $\\hat{x}=\\textit{fl}(x+2)-2=x+ \\Delta x$, $|\\Delta x|\\leq \\epsilon$. \n",
        "\n",
        "The second uses the inner product formula\n",
        "$$f(x)=p(x)= \\sum_{i=0}^9c_ix^i =\\begin{bmatrix} -512 & 2304 & -4608 & 5476 & -4032 & 2016 & -672 & 144 & 18& 1 \\end{bmatrix}\\begin{bmatrix}  1\\\\ x \\\\ x^2 \\\\ x^3 \\\\ x^4 \\\\ x^5 \\\\ x^6 \\\\ x^7 \\\\ x^8 \\\\ x^9 \\end{bmatrix}$$\n",
        "\n",
        "The third uses Horner's rule, which requires fewer operations\n",
        "\n",
        "$$f(x)=h(x) = c_0 + (c_1 + \\ldots (c_8 + c_9x)x \\ldots )x$$\n",
        "\n",
        "Each addition in the last two methods incurs a relative error of at most $\\epsilon$. An error in the innermost parenthesis, would correspond to evaluating the function at a slightly perturbed $x$. However, the error in the summation done last, contributes directly to the result. When $\\left|f(x)\\right|<|c_0|\\epsilon$, the result will contain no accurate significant digits.\n",
        "\n",
        "In terms of backward stability, extrapolating the above argument implies that the backward absolute error bound has the bound\n",
        "$$\\textit{fl}(h(x))-f(x)=\\Delta x \\leq \\epsilon \\left(1+\\left|\\frac{df^{-1}}{dx}(x)\\right|\\right).$$\n",
        "\n",
        "More generally, given their factorized form, we can evaluate any function with unconditional backward stability. But factorization is hard!"
      ]
    },
    {
      "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.4.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}