{ "cells": [ { "cell_type": "markdown", "id": "c7a000e8", "metadata": { "tags": [] }, "source": [ "**ACCUMULATION OF ROUNDOFF ERROR**\n", "\n", "Let's investigate how bad it can get..." ] }, { "cell_type": "code", "execution_count": 1, "id": "06d659dc", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from random import seed\n", "from random import random" ] }, { "cell_type": "markdown", "id": "780a55a8", "metadata": {}, "source": [ "Create arrays of random numbers that will be added to 1, then subtracted from 1. The **roundoff error** should then be the absolute value of the difference between 1 and the actual value." ] }, { "cell_type": "code", "execution_count": 2, "id": "2f3e3490", "metadata": {}, "outputs": [], "source": [ "nrands = 10\n", "boost = 1.\n", "seed(1)\n", "summands = np.zeros(nrands)\n", "for i in range(nrands):\n", " summands[i] = random()*boost" ] }, { "cell_type": "markdown", "id": "6cb1177f", "metadata": {}, "source": [ "Note that the random() command gives a number between 0 and 1. By \"boosting\" it to larger values, we can better reach the regime of catastrophic cancellation.\n", "\n", "Next, do the adding and subtracting. To make sure we're treating the random numbers sufficiently randomly, let's switch up the order in which we're adding and subtracting them..." ] }, { "cell_type": "code", "execution_count": 3, "id": "3cbe3ffb", "metadata": {}, "outputs": [], "source": [ "x = 1.\n", "for i in range(nrands):\n", " x = x + summands[i]\n", "for j in reversed(range(nrands)):\n", " x = x - summands[j]" ] }, { "cell_type": "code", "execution_count": 4, "id": "4a4014ad", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.0000000000000004" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x" ] }, { "cell_type": "markdown", "id": "41ac04d8", "metadata": {}, "source": [ "Try it again with different values of 'boost' and 'nrands'\n", "\n", "Next, for a fixed value of 'boost,' let's build up an array of the errors to see how it varies as a function of 'nrands'" ] }, { "cell_type": "code", "execution_count": 5, "id": "bda94711", "metadata": {}, "outputs": [], "source": [ "ntrial = 500\n", "boost = 1.e7\n", "fN = np.round(np.logspace(0.,5.,num=ntrial))\n", "N = fN.astype(int)\n", "error = np.zeros(ntrial)\n", "\n", "for itrial in range(ntrial):\n", " Nnow = N[itrial]\n", " summand = np.zeros(Nnow)\n", " for i in range(Nnow):\n", " summand[i] = random()*boost\n", " x = 1.\n", " for i in range(Nnow):\n", " x = x + summand[i]\n", " for j in reversed(range(Nnow)):\n", " x = x - summand[j]\n", " error[itrial] = np.abs(x-1.)" ] }, { "cell_type": "code", "execution_count": 6, "id": "7df3c301", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "error increases as: N^ 0.9885075585976549\n" ] } ], "source": [ "xxx = np.log10(np.extract(error > 0,fN))\n", "yyy = np.log10(np.extract(error > 0,error))\n", "coef = np.polyfit(xxx,yyy,1)\n", "error_fit = 10.**(np.polyval(coef,np.log10(fN)))\n", "print()\n", "print(\"error increases as: N^\",coef[0])" ] }, { "cell_type": "code", "execution_count": 7, "id": "fc24a2b0", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABX6klEQVR4nO2dZ5gUVdaA31udJpBhyDmIYsCAKEbMGNA1x881Z9cc1rCuYRfXVVddXNOaVndVUMwYQQUjYhbJCDKkYRiYPNPpfj+6q7u6urqne6Ynct7n8WHq1q2qWzLcUycrrTWCIAiCkA1Gay9AEARBaH+I8BAEQRCyRoSHIAiCkDUiPARBEISsEeEhCIIgZI0ID0EQBCFr3K29gJagV69eeujQoa29DEEQhHbFN998U6q1LnI6t1UIj6FDhzJ//vzWXoYgCEK7Qim1KtU5MVsJgiAIWSPCQxAEQciaDi08lFKTlVKPl5eXt/ZSBEEQOhQdWnhord/UWl/QtWvX1l6KIAhCh6JDCw9BEASheRDhIQiCIGSNCA9BEIR2zOL1lYTCLd9aQ4SHIAhCO2XhugoOe2AOD81a2uLPbrfCQylVqJSar5Q6qrXXIgiC0BqsK68F4MfiLS3+7BYXHkqpp5RSJUqpn23jk5RSi5VSy5RSN2ZwqxuAac2zSkEQhLZPOBz501CqxZ/dGuVJngGmAv8xB5RSLuBh4BCgGPhaKfUG4AKm2K4/BxgL/ALktcB6BUEQ2iThaBtxtTUID631HKXUUNvweGCZ1noFgFLqReAYrfUUIMkspZSaCBQCY4BapdRMrXXYNucC4AKAwYMH5/gtBEEQWh/TT94KsqPNFEYcAKy2HBcDe6SarLW+GUApdRZQahcc0TmPA48DjBs3ruVDEQRBEJqdqObRCk9uK8KjUWitn0l3Xik1GZg8cuTIllmQIAhCC/G/r35jY2U9sPX4PJxYAwyyHA+MjjUJrfWbwJvjxo07v6n3EgRBaEvc9OpPsZ+NVoibbSuhul8Do5RSw5RSXuAU4I2m3lQKIwqC0Jb4csUmht74Nt+v3tKk+wRDiZb61nCYt0ao7gvAF8BopVSxUupcrXUQuAx4D1gITNNaL2jptQmCIDQnHy0uAeDz5aVNuk990CY8mnS3xtHiwkNrfarWup/W2qO1Hqi1fjI6PlNrvY3WeoTW+i85epZU1RUEoc2gotu8bmIIT10glHBs9XksK6lk6I1vs3xjVdMe0gBtxWwlCILQ4cmVdcmueRiW+7714zoAXv+uyW7jtHRo4SE+D0EQ2hLmHq8bqXqs2FhFXSDkIDzi0iPf4wKgLpiUwZBTOrTwELOVIAhtCXOPz1Z21AVC1AVCHHjfJ1zy32+pDyaarawO87yo8Kj1h3hx3m+8+/P6Jq05FR1aeIjmIQhCWyLm88jyum1vfZej/vkpALMXlVAXSG22yvNEtvW6QIjH5qzgzR/XNnq96ejQwkM0D0EQ2hKN1TwAlpXEHeD1AbvmEf/Z545qHoEQ68pr6deleUoAtpUkQUEQhHZL8eYattQE2GFA+g/VmM8jC93D7+C7sPszDKX468yFbKysZ79tegFQUlFPXSBM364iPLJGypMIgtAS7PO3jwBYefeR6Seq7EN1N1XXJ41tKK+z3Vbx+JwVAEwY0ROAFaXVAM0mPMRsJQiC0Ew8NGspw/74duw4rnmk58lPf2XkTTMB2FTlTzpfUpkoPKw+D7MlbWlVROj0E+EhCILQvrj/gyVonVxOpCHV4863fiEY1viDYTZWJWsepTaBYvV52M1cfZrJ5yHCQxAEoZmpqg8C8XyMhjQPryuyNdf6Q46aR1l14phVNtnPdc33ZLnazOjQwkNCdQVBaElMk5GJKQQqaiPCI1201d3vLGLs7e9HrnNHrqsJBClz8HnY/SCmcIK4ucrE42qebb5DCw/xeQiC0JL4g2H+8cESDr7/E97+cR3+qEpQURcA4sLFKdrq0U+WU14b4JnPfo0Jgxp/iBp/KGmuXRuxChhTeCjCTDY+x7vsnRy8WTIdOtpKEAShJakLhHhw1lIALv3ft7HxitpE4fHwR8s5euwAijr72FzjZ0RRp9jcN6O1qSBitrInBEKyz6O00m/5uZ6Jxvdc536J7Y1V8N1y2C6pm3eTEeEhCIKQI+qCyVoCxDWPoMWsddbT8/AHw2yq9ieE+G6uiQuCGn8oqYIuJJutzONd1RJuLZ3Ozt4FrAr35trw5dx7yh2Nf6E0iPAQBEHIEfUOWgLEfR6hcPx8KKzZFHVuWyOkrCapyroAlXVxf4aJ1lDgdcVMWj2qlzPFM41DXN9QGu7KLYGzeSl0AAX5+c3WZlCEhyAIQo5oSPMIhJzjrH4rq4n9XB41cQGc++z8lM8q9LnpEVjHVe6XOdb4jFqjgHv8J/E/dQRbQl6g+Zzl0MGFh2SYC4LQFBatr2Bg9wI6+TLbKp38EwAVdabmERce1tyMb1dtznhNnX1ufPWbuFG/xWTvu4QxeDx0FOu3u4hnvk+MLPW4mq/HoERbCYIgOOAPhpn0wFwuszi+G8LJPwERh/knSzby3JerHM9/uWJTRvfvTA3Xeqfzie9Kfhd8h+mhiexf/w/uDp4KBd0zXmcu6NCahyAIQmMxndDzV2auFdibNFnHz356XsKY6QcB+KIB4eHDz5mu97nE/QbdA1W8Gd6Td4rOZebawticfK8r43XmAhEegiAIDpiO63SbciAUZktN3EeRSvMIhMLY8geptcxdV15H5zx3knPcRYgTXHO40v0K/VQZc8Jjeb/vBTy/qjt7F/QE4kLH7CBopam90tMhwkMQBIFIa1hrRz4z2c5pUza55dWfeWn+6thxKuHhVFbdzpkThvDwR8vN1XC4MY9r3dMYYazj2/BIrgpcwgLPTuzTuRewnnxP4vbtcxsYiiQh1Vx0aJ+HIAgCRL78F6+vTDtn2B9ncvVL38eOzUQ8U3iEw5plJYn3ePundQnHqZ6RifCYMLwXoNnH+Ik3vLfwiPdBgrg43381x/lv58vwGHweI+a8L7BpRB6XgdsWlptN35BsEeEhCEKH56qXvuewB+ZQbjExWTE39xnfrYmNbTI1j+gm/dRnv3Lw/XP4YfWW2Bx7Lat/fbwcJ/z2qroO9NjyI//1/JXnvVPooSq52n8Rh/vv5oPwOMxi7l6XQaFFeLxw/p6x6z1uA5eRGF3VnGardik8lFITlVJzlVKPKqUmtvZ6BEFo27wVLflhzd62Ury5Jmlsk6U67eqyGv4991cAFqytiI2HMtyd02keuxVs4MvhTzPm7WMZbazmtsDvObD+PmaE9yNs26K9boPOeRHhkedxMWFET4b1ijjNfS4Dt9F8obl2WtznoZR6CjgKKNFa72AZnwQ8CLiAf2ut705zGw1UAXlAcTMuVxCEdk61peKsU7Y2wKpokl7n6Fd9jT/I+wvWA5H6Uvve81Fs7iZL1dpwVPNQKv1XvpPw6E8pV7pf4UQ9F1VayMZx1zDx0zFUk5/yPh6L5mFiaj8et8Jly+toTvdHazjMnwGmAv8xB5RSLuBh4BAiwuBrpdQbRATJFNv15wBztdafKKX6APcDp7fAugVBaIdsrIxv9m/+uJbZi0q44uBRCXN+2xQRHn2iXfd+WVvByk01eF0GZTZtxaqRmJpHQwqI1WzVgwoudb/OGa4PAMWmHc6l16Q/UlHjo/rTT9Lex+My8EXLtevoQ2PCw8Hn0Zy0uPDQWs9RSg21DY8HlmmtVwAopV4EjtFaTyGipaRiM+BzOqGUugC4AGDw4MFNXbYgCO0Uayc+s8+3XXiYAsZ0RptlRIb0LGBpSVXCXGu/jEx9Cv5gmEJqOc81k/Pdb5NPPdND+/Ng8Hie3edYehV2xudPNp3Zsfo1TMEVjv7pdTBbbQ2hugOA1ZbjYmCPVJOVUscBhwHdiGgxSWitHwceBxg3blwLBa8JgtDWsGoeVrTW+ENhfG5XLMTWNC+ZX/NdHLrwbaryEwrrJOd0Knz4OaL6VU71TaOnquTt0HjuD57Icj0AILbhmw2g0uF1qVg3QtNXH4yZreKCJd/jiuaRSLRVAlrrGVrrC7XWJ2utP041TzoJCoKQSnic++x8xv9lFhBP2FtbXsvqshqC0eq3XfKSv6+/WLGJETfNTPm8HQZ0ASIJfie6PmaW71ou8z/JL+EhHF1/J5cGrmTszuNj883ihT6Xcz7JpO37ct4+w2JzTQFh+lvMP70uA3fU51Hoa/5s87YiPNYAgyzHA6NjgiAITWJjZX1CEUKAafNXM3tRCeW1AcJhHRMeW2oC7HvPR2k1DxOdwiaU7zY4zPiad7038nfP42zSXTjd/0f+L3ATP+oRAOw0MF5vz9WA5rHToK6xnA63y2BU707Re3QD4uYrq2Axw4u3BrPV18AopdQwIkLjFOC0pt5Ua/0m8Oa4cePOb+q9BEFon2ysrKdXJx91gVAs2uo/X6yMna/2B5P6cJg+j25phIdTefUJxgLu2vwKI7yLWBbuz4X+K3kvvDtmnoaJ1x3XDExtIZXw8LoMAkbkWW5DMW5oD2Zdsz/DoyG6pqDzug08UYd5obf5t/bWCNV9AZgI9FJKFQO3aa2fVEpdBrxHJMLqKa31ghw8S0qyC0IHpS4Q4qc15ew+tEfaeVX+IJ3z3HhdRkx41NTHy4hU1QcT6kxBfEPu1y112GzAEkG1g1rB9e6X2M/1E2Xh3lwXuIAZoX0J4Ww+spZKNzd8l6FwGSop8dDrNmI+GdPNYm1bG4+2UsmaR8rVN53WiLY6NcX4TCC1IbFxzxLNQxA6KLe+9jPTvylm7vUHMKhHQcp5wVAYt6HwWPIjqv3xfI+quiC1/kThYW7WA7unFx7D1Vqudk/nKNdXlOlO3Bk4g8rRZzL9h9K0a7dqGdbcjO4FXpRK9NN4XQZGVCgYdvsbFs3D4vMw/SipTGu5oK34PJoFcZgLQsflpzWRf9dV9c6JfyahsMZtGAywCIIaf4g8T2T7q6wPJnUAvGb6DwAM7O4slPqyCd87V/G+93oOML7nweBx7Ff/AE+GjsCTF7lmRFEhT5w5jgv3H550vdfS4c9jyc2YduGeXDpxRMJctyUE1ynCy5QPVp9Hj4JIJ8GLbffKJW3F59EsiOYhCB0Xc9MMp/i6DobCVNYFCYQ0HpeKlfGASNZ4r04+1gfqOO/Z+ZRVO5ct6VnoTTjuRiUXu9/gLNf7eH9RPB06lIeDx7CJuAPcLKRY4HVzyJg+LF5fkXiPAk9Ce1i3RfMYXtSJH4uTP3ZNjcNR89Bxn4cpZAq8LlbefaTjO+WKDq15CILQcTGFRqoy6He9vZBd7vyA6vogbpfB0WP7x84Fw5puBRFneCrBAXHzTwF1XOZ6lTm+KznPNZM3wxNYe8an3BE8M0FwQLzarel3sAqK+04cy+xrJuKxmK3siX1uW4mRsNaxOQ6ywzHD3Odp/q29QwsPMVsJQsfF1Ddq/c5FB9/5OVIMsXhzLS5DMXZQNx45fdfYeVN4pMOlAzw4bB5fFl7NtZ7pfBHenkn+v3Ft4CJqCwc4XpPvTSyZbvVvbNOnMz0KvQkOc2WTCHZhErYkJDppHiZWn4fPvfXkeTQL0sNcEDoupuZR43f2eXSNhtmur6iLbdbWroDd8r2O1wEYhDnWmEuPp/fimHUP4O07hmPrb+fCwNUs1QOBxGgrKzHNw5MsPMyffWmyyV22+lQhrWMO83RZ7V53cp5Hc9KhfR6CIHRcTFeHPczWpEteXLMwzTnWjdxZ89AcbHzLte5pbGusJpS3I0yeQUXRXnw3ZXbCzKBDngfEN27zT6tz3Hy+1ZRlJ9lsFSlLAuk1D49LxbJJ8kTzaBpithKEtkddIMTB93/CF8s3JYwvK6li1sINjr01nNAxzSOF8LAk+Hli5pz4ludzG7GIK4DxaiGveP/Mv7334SXApf4/UH/ObBh5ED5P8nd2qgZPcYd5as0jrfCwaRe9LE77dOW0XIaK1bnKE59H0xCzlSC0PVZsrGZZSRW3v5mYB3zw/Z9w7rPzOeDejxN6Ztj5tbSaksq6WGFAe44GwNcry5i9qCR2bJpzfAmZ3QYFXjdj1Eqe9vyNab47GaBKuTFwHof67+Ht8J64ovWmnBzQqUKEzV7iBVHfh1VgmVpIuiKI5lp3H9qdx/9vNybt0Df2rk5mq4O36w1EfCemNiRmK0EQOhym5SVViG0gpNlcE6BnJ8duCxxw78e4DUXfaO8NJ7PViY9+kXDsdti0i/zF3K0f5FDfp2zRhTzX+Vw2b38WL34cL/Adi15y2Ox//9Q8x/V5XAZXH7IN+21TBJDQvMl8vjeN5mEanwylOHT7vkA8osruXAd4+PRdqYpmzpt+mJYwW4nwEAShRTE3wnTJz9bw20AozPyVmxnTrwtdo36KYFjHNspUDnMrHrP4oMugN5v5g3sGp/70MX7t5p/B3/FE8EjuPGhvykoTTWbmh77Tpp3yWS6Dyw6M9wvpbPG9+DIwW5lY//eYgtbpMp/bha9TRFgEov9vJVS3iYjPQxDaHqamoHE2OUGiH+OpT3/l1Ce+5B8fLomVH7fOSRWqa8XtMqB2Mz2//Cuf+K7iZNfH/NTnOM7q/Dj3BU/iuL2255idByQ5q61C4+WLJnDNIdtk8KzEe3Syah4ZmK2cMN87ncMcIomREPe7NCcdWniIz0MQ2h6mwFhWUsV2f3qXX9ZWJM+xaB6bawJApBe5tQd5THgE0mseedRzUOl/4cGxdJr/MO+Ex3Og/17mjLqBGm8vIO5QT7c5jxvaI22hRBOPTXhYe4KYIbf2OVZiS7CoHqbMbEh4xMxWLSA8xGwlCEJWVNUHE76ms8Xuo1i5qZpt+3ZOnGPRPMzufmENW2ot/cPDztFW5rVugpzi+ojL3a/SZ8MW2GYS1fv8kav/FWkV5DJUkhnJHulkJ92mH5+T+E3eyaGhVCZmKyum2aphzcOMthLNQxCENsSSDZXscNt7zPi2uMG5q8tqGHrj20khufZyIrX+EBV1gZRzzK/pYDjM2c98nfQcu/DYXF3H0cbnfOi9jrs8T7NK9+HJbR6F017C3W/H2DyPS8U2cfNPowHhkc7RbeK2Jfk5mZDSO8wjaIvqERce6Z8dCIvZShCENsjCdRET00eLNzY495tVmwF4Yd5vCeN24bG+oo45SxNLmNc6CI+yaj8rNlYnPWfVpuiY1rDkfbo/fzAPeadSi4+z/Ndxkv9PrO86FkjctF2GkdTFryHNw52B8LBrJ07O9oaElB0zpaShvulxzaP5t3YxWwmCkDFmhNTGyjpuevUn/jx5+5TOXzNJLpVZyeTv7y1OurbG0WyVHJ41uk9nlpZUsumXOWx582ZG1P4InQbxB/9lzPXtS9fOPthUE9v0rZu2x6Vi/gVvhpqHXTAcsWNfZv60PmEsm8gsJ4ZGq/9aCzma797QvQNitsoNEm0lCLnF3MS+XFHG/776jfd/WU95bYC/vbsoKWTWTJKzj9cGGo6OsmonZiZ3IJgsPI7ut5nH3PfSc9pkOtf8xt/dF/LYji/xlt6L968+gF0GdwecNQqXoWKbsSkUGvZ52PwZjr6fpjVg6tMljyV3Hc4Zew6JjY3p1wWAsQPTB/+IwzxHSD8PQcgt9hap5bUBps9fzSMfLwfghknbxs6ZGkmS5pGiFlXCHAfNw7TnAwxSG7jK/QrHLvqMSiOfewIn83ToMHrkd6fwl1LGD+tBUWdfzMxj90NExuK1oMwS6a4Gvuztpzv5Gq7M2xjs2twB2/ZusGMixEN1Jc9DEIRmobSqnn3+NptlJZVZXWe3HFXUBtlcE4mA+tTmtzBrT32/ekvCc1L137Di5PMIhjS9KOf+Ts8zy3stRxhf8evo89i3/gH+FTqGWvLYXONnfXkd2/aNfKmbioQ99wIiAsXI0mwVFzcROvniX/hH7dQPgN5d8pKu22NYD0YUFSaNZ0NDggPgd7tEysRLYURBEJqFD37ZQPHmWv4999esrrP7HSrqAiwviTisy2sTI6ZClrknWMqFWLWKAdG8iS55bu45YafYeI0/hNaa8toA/lCYztRwWtWzzPFdyTHBd5kWmsh+9Q+wacJNlNMJiJieavwhqv2hpLwNq68ipo24VOx8pg5zbTNJWcNwz9lnGCvvPjKhmq/JSxdOYNY1E9PeOxfcccwO/HDboVknITaGDm22EgTBGXOPTFcixMqCteWM6dclVv7CpKI2wLKNVUCib2NdeS31wbiZyZrcZ9UqhvUqZM2WWgb3LIj13YaIdvLlijLOf2oulxTMYqpvBt3rq3gjPIEP+pzLm8WRr3Brm9j+3fJZtamGUFjHfBOmT8PaI8NtKELRBkumGSrTUF07VrNVJmG8djLJG8kGl6FifUyaGxEegrAVYppf7F/STsxetIFznpnPPcfvRL3N5PRD8RaWlZjCI3KuLhBiwpTZMa0C4pFXEBcex+06gM4+N58uK6XQ605ozVpf78f9/bN84P4n/QJlfBwey/OFv+fDLX3Zv6AIiIQK9yyMF0/s1zWPVZsitalMYWDu51aNwm0o6omYrZStN3hDmof9f1ehxWyVrsGTE1/88cAWMS81F+1SeCilDOBOoAswX2v9bCsvSRDaF7HKtg1PNc1SizdU0t3WQOnnNZG8j1N2H8SLX6/myxWb2BL1gazZUhubV+iNbzW1/hDb9+/C/SftzMMfLQMim73HUCjCHG7M49bfZtAvWMw3ehRX+i/lK70dA8gHahMinLrkx3/ubxFWptnG7qMAi9nKUJaku8wyuM1Kviad85Ir5mZKv64Nlzppy7S48FBKPQUcBZRorXewjE8CHgRcwL+11nenuc0xwEBgE9BwqqsgCAmYm2QmZitTO1FAnUOYrdtQDIvmJpzy+JeO9yiwfKFvqfHHuvgVRcuu1/qD9NzwKW94b2NHYyWr9BBeGnkPN/w8AFPSmY5zqxZjzXuwajpxn0fk2OqrMXM+XBafhylEG9I8hhd14uNrJ3LFi9+xqqyGEUWdYudaws/QlmgNzeMZYCrwH3NAKeUCHgYOISIMvlZKvUFEkEyxXX8OMBr4XGv9mFLqZWBWC6xbEDoMZpXWTMxWsaJ8hnKMlHIZKmFDd8L6RV9W7WeHbpF8haLOPnZWy7i17GVGf/AjqyniKv/FvB7em/DPiZuxKTycakUB9LL0/zBDc5VNOETORZ3ohpHUW6ShDG6IJPG9duneQKL/pjE+j/ZMiwsPrfUcpdRQ2/B4YJnWegWAUupF4Bit9RQiWkoCSqliwKyQ1nDcnyAICdQHo2XRM9E8onMUUBdM/ufmNlQsIdCKK+qYBqisC/DvuSs4afdBbKr2RxzdJQvZ7YvbeM33HmXBrqzb6w4Omj0EP84OX7P0xi6Du9Oz0Etfm9nH2pM8lrdhmBqWRfMwVOycadZK16nPCbuvBLY+zaOtvO0AYLXluDg6looZwGFKqX8Cc5wmKKUuUErNV0rN37ix4To8grA1UW8p+VEXCHHFi9+xusy5d3hMO1FQ72C2SqV5mIJj75E92VBRz11vL+SaaT/QpW4dJ6+5Gx7Zi07rPufewIkc53mYqrHnpBQcEE8SzHNHmi2dsNvAhPPWnuXeNGYrlysevmv6MDpH/SiZCg/ru8eeuZUJj3bpMNda1wDnNjDncaXUOmCy1+vdrWVWJgjtA1N4aA2zF5Xw+vdrqQ+EufHwbWO1lUzMfddQirqgg8/DZVCQpkR7t/xIOG1Pyjlo5X+Z6nsPd6kBe16C3vsqlsxYyd/3G95gmXJT80i1wXezCI9Y6G1UMwhZlu0x4lrJdYeNZrt+XZg4uijtvVNhzUjf2sxWbeVt1wCDLMcDo2NNQppBCYIzpu8iGA5TXR/JwXh3wXom3vsxP69JrAUX72LnnB1uqPQ+jyKvnyvdL/OJ7ypODM1kRmhfPj70PTjsL7g69eLxM8ex+9AejlngVoINdNPr4iA84j4Pq8M8XrIkz+PihN0GWvJBshMe1ryQphZEbG+0Fc3ja2CUUmoYEaFxCnBaU2+qlJoMTB45cmRTbyUIHQpT86gLxIWHyW9lNewwIP7BZW67CmeHuTuF2cqHnzNcH3LN4jcpcJfzVmgP7g+eyArdnxeKhiTNz/TLPVUinzUcOK55RN/BaraKOdMd7r2VCYCm0Bqhui8AE4FeUcf3bVrrJ5VSlwHvEYmwekprvaCpz5LCiIKQzOvfr+HxOSuAiCZRbStcqIhENpkbcMxhntbnEd9KXIQ4zjWXK92vMEBtYlWnPbhsw2R+0sNjc4o6+5Luk2l3vVTFC60NkMxQXfM53SzZ66bDPOiQ5GIN1T1n72EZrWdrpTWirU5NMT4TmJnLZ4nmIQjw1YpNlNcGOHT7vgBMnb0sdq42EKLKpnlMeWcRF//3W5b/9QhchkroJVEfcvJ5qOhmrTnM+Jrr3NMYaazl+/BwrgtcyD47H89P78Z7drgNxZCeyUX+PBk6nFNZlvK88etNQXTmhKF08rk5bte4c900TYXCzoIQIkLnT5PHZLSerZW2YrZqFkTzEAQ4OZq4t/LuI4HEXg+1/hA1DmYriITzFnjdFrNV6g23f9k85va4i0E1C1kaHsCF/qt4LzwOUBxo0yiG9Sp01DKc6jwZKjkLPpXZymr2ipcnUZw4blDCPDMfxMlEZQoPMV41TIcWHqJ5CEIy1hali9ZXsmi9c1n2FRurI5tpTPOIRzyZ7KhWcHvtdIznfmBQl4G8NewWrlg4mhCpaz5t06ez4/M8tp4bewzrwTWHjuakx75IGE/ll1AJORept/+/n7ATr3xbzM6DuiWdy9ZhvjXToPBQSuURSdTbF+gP1AI/A2/nwi/RnIjmIQhxPlpUwgHb9s64y9wdb/3C5mo/k3aImLsMFTdhjVBruNo9nSNd89gS6gKHTYFx57Dko1WEFi5LuI9VyxjUI5/T9xjs+DzDUFx58ChCYc0/Zy+jc56HAd2T6z/ZlZZ/nzmOz5Yn9hJJ5z/pXujlvH2HO54T4ZE5aY2MSqnbgc+AvYCvgMeAaUAQuFsp9YFSaqc0t2hVpA2tsLVjjY46+5mvqawLZFz9ddWmarbUBmICw1DQPVjC3e7Hed97PfsbP/JA8DjO7/ZvmHAJePLwRQWTVbuxJs89cPLO7DWyV8pnXnnwNrGWq4aK1Kt6+qzdOXRMn9gcu+Zx8Jg+3DZ5+4Qxp86BmdBQJ0EhTkOaxzyt9W0pzt2vlOoNOH9GtAFE8xC2djZW1iccaxLrMf065QgmTJnN+oq6pGtLKusp8LgIa+hGJXuveIALql9CuzTPhCbxcPAYyujCWHe8OKDpt/AYBnWEo2PxjdyVwaYe87FE9/EDtu3Nmz+sjZ3PJJw2ndkqHTGfh8iQBkn7N6m1flsp5VJK3ZvifInWen7zLE0QhKaywSYUdBiq6uPCQynFpzcc4Giu0RoIVDNh9ZPM8V3J2OL/8YFrXw6sv487g/9HGRENwRreaob1WhP+rJpHJl/21ox26zpj98jAtJRp2K8dMVtlToM+D611SCm1T0ssRhCEphMIhflm1Wb2HN6TJRuqEs6ZGeWDexTw1Fm7A5HyImZ3PRMvAU51zeYy96sUrang3fDulO1xPVN/crGWRIFk3XDNO1g3b2+C5tHw5nzQdr05dpcBXHfYaMt18fOZaB6NFR7mrZ36gAiJZBpt9V20RPp0oNoc1FrPaJZV5QiJthI6AuW1ATZW1jGydyRK6dvfNtOva17KZkJ/e2cR//70V966fB8+tzmSQ1pTVRdkv216MbJ33NxkdtczCPM741Ou9rzMQFXK56Ex/HfE5TywqCvX5A8jGF6V9DyrPDC1Bo9D2CxkJjzyPC7+cfLOCWPW6zJxZzRWeGTallfIXHjkEWm8dKBlTBOpbttmEZ+H0JYpq/ZT6HPha6AV6QmPfM7SkqpYnsZx//ocl6FY/tcjHOf/si7S3W9LTYBvVm1OOBcOQ3V9kEJbIUOXAYcY87nWPY3RRjE/hYdyY+B8Pg3vwKFGX2ADIa0TakQ5YVbgTWm2aqRZKMFslYnPo4WLFDolPXZ0MhIeWuuzm3shgrC1seudHzBxdBHPnD0+7byl0R7hWuvYJhpK0z/Weq6ksp5CrytWgiQQClPlD8ZKkAOw8lOe51Z28i5hebgfl/j/wDvh8eioS7SyLpJEGA5rx5IeVmI+D4uQsCb/NVZ4uFL4P1LRUJHFVJg1usYP65HxNT/cdmjW/cs7AhkJD6XUQOCfwN7RobnAFVpraQErCI3ALNT38eL0vWZ+LN4S+7nGH8ooR8PUDs548isABnTPj/k+ymsDaB2tQLvuB5h1Byz7kD704IbA+bwc2i8hwQ+goi4ARExe6YSW9b0SfB6WjbWhNq+psF7WnA7zbgVe3rtyv6w0ia75qXuQdGQy/T/8NPAGkSTB/sCb0bE2jeR5CG2Vhr7gTY6e+lns52p/kBp/vJSIvRquif3W/S29vZeVVDFUreOIxTfDY/vBmm/gkDs52fswL4UOSBIcYBEe4Yj2Yf/wt1qyDtouko8xeWz/2Jj1qzxVaZGGsF6XidnKqdRJpozu2znjRMqtmUyFR5HW+mmtdTD63zNAUTOuKydIPw+huSmvCfDNqrKsr6t3aKrUENX1IWotFXBXb3bu/GfXDkzh0ZvNVL9yOR96r6PP+o9gv+vgih9g7z8QdOWlfG55TUR4aB0xW6Uz0WzXrwsr7z6SXSylP6wJe43VPBLNVqnn7TK4W3SOREs1N5k6zDcppc4AXogen0rEgS4IWzVnPvUVPxSXxyrQZkq9Q1+MhrBrGlVRX8SWGj9V9UEGdo+YWrTNqT280M+N7hc4y/UuBmGeDx3MkRfdR1G/eH5vurVXRJ8TCkcc5gUeN3UOpdmtWO+XECnVyE3dSHE/O8+duwfry5MTHoXck6nwOIeIz+MfRKKsPgfEiS5s9fxQHDGJhrXGlUVuwOYaf4Nz6oOJAsYuPMwN/OD7P6G0yh+LxgpFhUc+dZzteo/fz5uJy1XFq+G9+UfwBIp1b87sk1hp1imKqmu+h/LaQOw4lELzcDLAWR3WCT832ueRmfDo5HMnhCALzUcmhRFdwF+11ke3wHoEoV0SCmvSmclvfOVHOue5ufnIMazdUsvB989xnPfmD2t5Yu4Kpl80gbLqRAGzfGN1wkZslhkprUqcp0IBznB9wB/cr9JbbWFTnwM5bcWhLNZxTcPuezArrXfyuWP9PYYXFbJgTQX+aA+PUFijNfg8DVu7rWVIEnM0mu4wF4tU2yDTDPMhSimv1rrhz6U2hCQJCi1FQ/kPL369GoCbjxzD4g2JJdCLN9cwsHsB/mCYy1/4DoDVZTVJpqGbXv0p4bjWbvoKh+HnV3is4hb6e9bzVXhb6o57hrWdd2Lx8i9j0964bG/smH6SGyaN5vQ9hjDuLx/Su7OP2qLCWMn2QFSINJSXAokaRk58Hlk6zIXmJ1OH+QrgM6XUrUqpq83/mnNhuUAc5kJL0VAIq5UKiykIYJ+/fQSQ0NFvzZa6JM3DTp0/xO8e/gzQHGB8h35sH5hxHjXkc5b/ek7234oeOD4hVBZIaBlrYkZ/ed0GhqHYcUBXdhncPaH3hunkzyQBL9Vm39g8D6vZSvqMtw0y9Xksj/5nAM6dXARhK8baYK+yLkAorBP6ZpvU+IOUVNQnjUOiAFq7pTbWk/uViydw/CNfJM3/dVM17uIvmeZ9ifHGYsL1w1DHP8mF7/VgeU0tECn1Yd/snb7+zQ6Bplbx7DmRxMWps5fG5gSijaAyMVtZn+HKQZJggvCQ4oVtgkx9HttorU9vgfUIQrskpDWlVfUs2VDJxc9/S3ltIObADlj6fq/dUsva8lrHe1hNX2s219KvWyR81qmG1bbqN/aZdx83+L6hRHfj5sA5XPH7O3jtpxKWly6Kzctzu5Ic3E4beMiieVjpa3m2WRY9k2xql+HsJG+syckq/6TybdugQ/s8BKGlCIU1Zz41j4XRulJWXor6OwBWb65l3ZbkUFJty95eu6WWbgWRzOVOeW58boP6YJhBagNXu1/mGONzKkP53B08hWdCh1GHj/NCBn+duSjhvj6PkfSl7lS6w3y2XTD07JSsPSX5PBwsdu4cO8xVgtmqUbcQckymZivT5/EGiVV172+WVQlCG6Ss2k+XPDduB5t/WOskwREMhXEZilte+zl+jyo/823FCiHiT7AKj801/lgXwDy3iyHeSs5gGqe6ZhPExaOhyTwaPIoK4mGpNf4gYwd2jYUPQ0QY2H35jpqHNoVHomDoVehLmpscqpssPRJMVTnwUeQiV0TILeLzEATguuk/8EPxFt6/av+Uc3a98wMOGF3E0w6FDG9/c0HS2KZqf1Ldoyc//ZXSqmSfR1V9MEF4hHQkj6ObUYPn47t4PTwVtyvIi6EDeCh4LBvpnnSPF+etZlVZYta5Uoou+Yn/zJ1atKYyWzlpHvY5TrhzoG1Yyba2ldD8ZFpV93YApVSB1tq5JkILopTaFzidyPrHaK33auUlCe2c6d+kr/EZjm6uH6UoZDjzp/VJYxsq6pIK9P2yroL9tynikyWJ96muD8a+/gHmLSlmu+VP8YnnDdSn1czz7c+fKo5hle4bqYjrUNfquS8jvTY657ljlXABenfO44GTd+bKl74H0vs87FpFj8IMzFYO5HqDl2irtkdGobpKqQlKqV+ARdHjsUqpfzXmgUqpp5RSJUqpn23jk5RSi5VSy5RSN6a7h9Z6rtb6IuAt4NnGrEMQMmXtltqEjHCz1lNDlFTUOxYvHDckWWt4/stVhMMaN0FOdc3iY9/V/NHzAj8Zo+HCuTzc/Y+s0n0B+MNBo2LXOW3SvTolm5pG9Uls/GTHVHrsWoVTgcBMoq08mXRsyoJEs1VOby00kkz/hh8ADiNaz0pr/QOwXyOf+QwwyToQjeh6GDgcGAOcqpQao5TaUSn1lu2/3pZLTwP+18h1CEJG7HX3bA57YG7s2Kwy2xA/ry1ni4Ogcdp8/z13Od5Fr/G+93qmeJ6kWBdxUv2t3OC7FfrthMcd2TGvnzSa8/cbHruuyEFQOAmUhES9NBVnM4mkyijPowlVbZ1I1c9caD0y9XmgtV5t+0vLvrJb5D5zlFJDbcPjgWVa6xUASqkXgWO01lOAo5zuo5QaDJRrrStTnL8AuABg8ODBTlMEISXry+vI97joGo14svopqlKUQrfzwIdL+eCXDUnjPreLkb07saykCtDsZ/zI9e6XGPLRShYxiHP91zArvCugGBC9xjR/9bSZkXp19rK+IjF6y0mzsO73Tj4PEyd/xgvn78mpT8Qz1O3Czym5vrGZ5KmQ3I62R6aax2ql1F6AVkp5lFLXAgtzuI4BwGrLcXF0LB3nkqaniNb6ca31OK31uKKiNl89XmgjmL6NPafMYuwd7zvOMU1R9uq1VgZ2j+RHLFgbicC6YdK2sXM+t8H/ztuDyT2KecHzF/7j/RtdqWbxXvdxhH8Ks8K7AYkdA80N325G6u6QiOikWRgZhro6+TMmjOiZUE+qdXweOb2dkAMyFR4XAZcS2dDXADtHj1sNrfVtWuvP082RZlCClWAo3GC5bn8osZ6Uk3/j+pd/5Nxnvk7bk+OsvYZyxI59Y8f7juoV+7lX7Qp6v30O/6y5npFGMbcGzuJA/32sHjiZMEbCV7vpRPdGzVZ2Z7HTJu1y0Cys2kY6s08qf4ZVTmaUJJhj05LUs2p7ZBptVUokuqm5WANYa0QPjI4JQqP4tbSaYb0KY8daa3a760PKawN8duOBDOiWnLUNUBdIbPVavCU5uHBFaTUrSqtjeRgmdx6zPf/6eDnryuso8Lpj/TUACn1uBqqNXOl+hYM+mgu+znw6+CIuWDKe7t26E9hSS120BLvXbRCMNn0yNSHTbGVmq+87qhdf/VrmuKm6DcVH105M0Iwy9V9n4s/IpCR7rs1MYrZqe6T9TVFK3aKUStkJXil1oFLK0SeRJV8Do5RSw5RSXuAUIm1vm4QURtw6mb1oAwfc+zEzf1oXGwuFdaw3xeqy1NHmdm3CyeFtYq96e+RO/endJVJSJN9rxARUT8rp89mfmO29mqOML/ht9DlwxQ98MeBsasiL+Rnqo/ez+h1Ctp7gwWh9qefO3YMldx3uqEW4DcWwXoUML7JGWGUmPTLRKnyt0KJVwnPbHg1pHj8Bbyql6oBvgY1AHjCKiOnqQ+Cv2TxQKfUCMBHopZQqBm7TWj+plLoMeA9wAU9prZOzrrJESrJvXVTVB5m7ZCPLN1YBsGBtOUfs2A9I7OudrlptXSAU+9pvaK7fJmhchsIb9Tfke1x0M+q5yv0y57neJv/7IC+E9uOh4HHcO+4IhhT0IKxLgHi/bVNwWXNDTJ/H9ZNGUx8Mc+RO/WzPTF6Xo88jQ80jk0gmu4C59IARmd28CWSgEAktTFrhobV+HXhdKTUK2BvoB1QAzwMXaK2dK7ylv+epKcZnAjOzvV8Dz3oTeHPcuHHn5/K+Qtvkz28s4OVvimO+BmXp7GctOrgprfAIJ2gfW9J0/LN3+nMbCo/LwIefkcufZejCRznAvZm3Qnty+GVTuem+JUDcr2AuyRQWZr8Oq+nIFGS9O+fxz1N3SVpDpj6PhnwGF+4/nMc+WZF2jolVeJjFH5sb0TzaHpn6PJYCSwGUUgbQqTGCQxByzaaqek7/91c8csZurNkc+ZU0O+tZ9xurw7fM1nnP6huoD4YSfBmb05it7CYuFyEOqn2Pe33P0P+7MsLDD+SohQfysx7OyqJRQFR4RDdfsyaUPQvd6rQONdBkymlT9TSQ5+HEHw/fjj8evl3aObH1ZWDayjUiPNoemWaY/08p1UUpVQj8DPyilLqueZfWdCTaquMza2EJi9ZX8tCspeRFN12zNIe53awsrebX0lg9T8qqE2tLWSOs6gLhhA59ac1Wses0k4x5+B7fm3PL7qdEd2fefs9gnPkqP+vhSdfFQl1TZHUnah4pHw+k0jwab7bKhExCdXON1LNqe2T6KzVGa10B/A54BxgG/F9zLSpXiMO849O7SyTDes3m2liUVCyhL/q1OvHejznioXiGeKlNIFg1iPpgKEF4bE5jtvIHw+xl/Mxr3lt51PsAKIOf9vkXv/PfQf9dDgPg+F0HMryoMOE688s9rJ3rSVmP9x+dPkfJMdrKweeRy823dTSPFn+k0ACZZph7lFIeIsJjqtY6oJTKvO9mKyEO846P+UtYvLmGvl0jkU4bKyPC479fruLE3QYmXWM3W9UHbJqHP7XZylAR5/tOajnbvDeV/3k/p1j34trAhdx7yRR2NFysPDg+/76TxiY9P5XPw8TURIYXFTr6ORLWE91VD96uN1+uKKOqPuhoosql8Mikqm6uEbNV2yPT34LHgJVAITBHKTWEiOO8TSOaR8cnENUaSqv91PgTHdibqv2c/Fhy+9aSysREQavjuy5g83nYtJRd8jfyL88DvOG7lcLNC7k98H8cWH8fL4f2ByMzc45p9jlhXESwHbljYgRVTHj06uRYmNCKqXn0KPRy+h6RMjyO5UlyuPmK8BAgc4f5Q8BDlqFVSqkDmmdJwtZKOKwJhnVWm1MwGo3kUopKh4KFa20Z5T63weqyWkJhHfsaTzRbOfs8+lPKFe4ZnBCeQ63h5f7ACQw+9FqefiOzCCX7GgC27duFlXcfyQ+rtyScN30emdSHMpUMl6Fi5qpM/SCNpTH3Om7XAc7ZhM34TKF5yUh4KKV6ArcB+xD5FfgUuINold22ipit2hfXvfwjr3xbnFX4p5lx7TJUQg+LVIwo6sQv6ypYu6WWQT0K2FRVz1xLb426QCjBbEVNKTe7X+VM1weA5p2CY/hT2WGU0YWbg8kVbTPB7jOw+yhMM1YmlWnNL3JDqViIrpPPI5eVaBuzkd9/0s5NeqbZz11oO2Tq83gRmAMcHz0+HXgJODjlFW0AyfNoX7zybaQhk71EiBOry2oo8LpiiXqGIkFjSMXwokJ+WVfBvvd8xE4Du7J4faWj5lFILee63uF89TYFrjpeCe3HA8Hj2X7E9pSVRSrlVmZYXdeOvY1tKp9HJpqHuZG7DGUJ0W3er3SF4vpJo9lzeM9mfY6VbXpLA9O2RqbCo5/W+k7L8V1KqZObY0GCsL68jqG9CtPO2feejzAU/OXYHYHI5mnP+HbCWrLjR0uvbxN/fS2D177IJ75/0UtV8E5od+4NnsRyHSnyvK+liq095Lex2IWE6RPJxM6foHnENI7mjWUJhTWXTGxYm/e4FIN6FDQ4LxMMQzGgW75jC1+hdchUeLyvlDoFmBY9PoFIKRFByBqtNX9/bzHH7TqQkb3jm7nZPnVdBsIDIlFPVrNVIJQsPHxuI0Gz6NfV2fxhEOZY41NOm3cNXevX8Vl4ey7hNOYFhlHU2QfRCK5uhfGe5M9/+VtmLxzlgNFFjm1sU2kemZCoebSMIzvcQOKiyaI7D8/pc2dfu79j7xChdUgrPJRSlUQ+YxRwJfBc9JQLqAKubc7FNRXxebRNKuqC/Ovj5bww7ze++9OhsfGu+R4q64JsqEhfNt1K3GzlLDx6FHpZZ3Gau5Qiz2NYihpqDjXmc617GtsYa1jn2pZHiq7g0dWD6d81D8rrGNazMBb+69Q/I1OeOHMcgVDy7peUYd4I4WGolnMqh8KZ7eC5Xk9rJCcKqUn7W6q17qy17qK17gz0IuIwPxg4AJjcAutrEhKq2zYxNxV7DkXX/MhXvbUzXiis+WxZKUEHwQDENmOXoWIVZ6108iV+HykVH9vT+IVXvbfxuPcfuAhzsf8KrujyDx5dHQl57RoVFEN6xk0v3Qs8NBa3yyDfm7wB2h3cpuaRydYbM1sZKlZgsTm+zqdfNCGW7JhKexO2LjKNtjoPuIJIn43vgT2Bz4GDmm1lQofgn7OWUh8Mc+1hoxuca26E1oZN7y1YzyX//Zb9tiniP+eMT7rGFCqGUvhDYbxuI8H3Yf/6NZTi2cN9lL52B/u7fmSt7sH1gfN5JbQfIVxsEy3bfuiYPjEh1t/S+6OTr/HCIxV2c1MmPTVMzNdzKdWsPS92H9qDWVfvT30w3GAwg7B1kKnP4wpgd+BLrfUBSqltybIUu7B1ct8HkWKAmQgPM1rKKjxWborUpPpmZZnjNaapSqnIz13yPQnCw+p0HqbWscc3/2Hg2nep8nTmrvrTeS50CPVENIxOPjdV0XDf0/ccwuqyGn4sLmfHAXHNtdCX+43T43YO1c1E9bD6PFRzR1kpJYJDiJGp8KjTWtcppVBK+bTWi5RSDe8GwlZNJrbxZz9fyZE79aNXJ18sv6LMUk9q7ZZIpdyaQIhgKJwU5uqPmqoCoTBhDQUeF1uIm8NchqIPZVzhnsFJro9hgw/2u54bV+7FW0sSm0L17ORlk1mRFzhjzyH8bpcBCaavznmZ/pPJHHs5kWwc5mb+Ri7zOAQhEzL9LS1WSnUDXgM+UEq9DqxqrkUJHYPizak79pnc9sYCLnruG4Cktq4Aa7dEtBCt4bPlyTmppuZhRlRZfQpdqeLsmqf4xHcVJ7g+4bnQIcw69H048GYCnuS8gR6FXqqiuRumxmL3mWzXrwvb9s1tzoEnhc8jE0RkCK1FRr+lWutjtdZbtNZ/Bm4FniRSJLFNIyXZW5elG6ocx+ttQmLtllpWl9Uk1aaau3QjsxeVxPpy/P6peay0lFYHqPFHNntT8BT63ORTxyWu15jf+Rp+VzuDt8N7cqD/Pm4P/p5Afi8gOVEPoGdhPJIq1Yd8gdfN1NPSFyvMFrvW0NjaUZkoH/YKv4LQWLLWwbXWnzTHQpoDyTBvOX7bVMMZT37F9Ism0Cfax3tpSUR4DLYlilXYyoisLa9j33s+Shj7csUm/u/JeQAUdfJREg2VteeAmD3G6wJhPAQ5xv82R/v+S5Eqh+FHcPG6I3inpEdsvqlRmIl5B27bm9mLIu1gu+THneHpNmJ7aG2uySSz3MRpnamirX7686HNvnZh60F+k4Sc8L95v/FbWQ0vf1McG1taUgkk+wnKa1N35zP59rfNsZ/P3zfeUMmeA7KlJoAizDHGp8zyXsO5Ff9ihe7H2a6/wKkv8KsxJGG+uS+bfgarychqokrnfG7uqrLZVMDNJiy3c55HHN5Czsi990/YKjEFhLU44bKo5lFSWU8gFI599WYiPMyIqRV/PYKQ1vxl5kIAvl+9xdJXXDO6/DP+5H2a7YzV/BweysP9r+XvKwbSyxvRfuybq7JpHlbzlVV4pPv4b+6vd1N2NHf0lCA0BdE8hJxgCo+q+rhgWF5SxcDu+WysrOedn9fHxisyFB4el4omvxl8e+sh5HtcPPP5Sv785i/srhYx3Xs7f66+gzz8XOa/nMn+u/i1+wRAxTZgeymNmNkqqnFY+313smhI6aKXrJrH9386pMF3aQkUcee5bubaVoIAIjyEKCc++nmSMzpTZi/aEMvNsGoedcEwR4/tz/BehTw5N973wtQ8RvdJHbVUHwwnJMv1KPQyoHs+26lVPOW5h+m+OxisSrgpcC6H+P/OW+EJaAwKvGZRwch1ycIj8qc7lh+RmeZx+A59mRhtCetNoa20JiIuhJambfzmC63KovUVfL1yM1e+9D2vXbp31tdfN/3HWHZzZV2QN39Yy+UvfAdEvvSP3rk/D3y4NFZq3RQelx80ksv+953jPf3BMD6rfb5sBX8O/IO9vB9TST5TAqfy24jTeWdJYkPLAq879lxINlvFNY9kn4c138KueTxyxm6xn70ug56FXq47bLQ0KRK2Wtql5qGUGqyUek0p9ZRS6sbWXk97J2hJtMv+2jBlNf7Yl3plXYA/vJgoEOwF7SrqAnjdRlL7VSt+U/OoXA9vXQ1Td2f3us95JDSZfesf4LHQZHYe0T/Soc6CqXmYW7pd8zBlgrXr3qFj+rDL4G5YXRnpfNaGofjm1kM4ZfzgnCXnzbhkr0ZdZ318YVQL6pqf+xIqgmCnxTUPpdRTwFFAidZ6B8v4JOBBIhV7/621vjvNbXYEXtZaP6+UeqlZFyykZXNNAK3BlDsbK+sbjACqqA3QNd+DUoojd+rH2z+ui50zTUKGv5xLw/+FB2dCOAC7ncWtJYcxbXFi9vj9J+1M/675TP1oGWARHtFd1Z7kbg/V9bgMHj9zHADT569OmtdS7Dq4e5PvceSO/SitqufU8YNzsCJBSE9rmK2eAaYC/zEHlFIu4GHgEKAY+Fop9QYRQTLFdv05wJfAy0qpc4iXiRdagU3Rhkhmk56VmxrOKi+vDdAl6pwutFWZ7e4JckzVNI7dMI2CcBXseCIccBP0GI7/xe+AtbG5sXLkFtORvWptSod51ERlNTtZBUZ7NEYZhuLsvYe19jKErYQWFx5a6zlKqaG24fHAMq31CgCl1IvAMVrrKUS0lASUUtcCt0Xv9TLwtMOcC4ALAAYPli+x5qK00p80lu9xJbSENb/yA6FIRdaK2mDMtGL6KNwEOcn1CVcxg6LKzfyQvwdPeE5n6vG/j93Hbv6KOb0tm37MYR41QSX7PBKvtQqJVIKkpcnG+S3NkYTWoq04zAcAqy3HxcAeaea/C/xZKXUasNJpgtb6caXUOmCy1+vdzWmO0HQ2ObRi3W1IdzbX+FmwNuLMNgVFeW2AznkeymsD9OoUKQVS4FFMNj7navd0hhkb+EFtyz973MyvBTvF6kyZ5HkSXXSmxmH1Wed7TJ+Habay+zwSHebWs1YNpi3UGcxmDW1gucJWRlsRHlmhtf6ZSCvchuZJeZJmprQqWfNQKiJAFqytIKQ1XaMNlLbUBBjYPSJEhvcqgKUfcvaCmyjyLmZheBDTtrmPGZVjCKPAFqoLJGVHux3MVuac1HkekT+d8vysGkxrCg+znMuYfl0yvkYUEKGlaSvRVmuAQZbjgdGxJiGFEbOnqj7I1NlLU3buszJ9/mrufOsXx3OmA3jJ+kq6WTQPrTW9Nn/PpauugP8eT16oiiv8l3CEfwrD9joutmv7g+GkMiA+m/AwTUtWc1Nc84iQ5DCPmauSpYNVg2nN7O49h/fkrcv34ey9hzY4ty1oSMLWSVsRHl8Do5RSw5RSXuAU4I2m3lTa0GbP58tKuff9JXyzanODc697+ceU50zhsXBdBd2i7VyD6xZQ9eyJzPD9me61q+CIe3lv4lu8Ht4HbftV9AfDST4Ou9nK5WC2GtWnMyN7d+LO30UC+VL5PJywajBGK//L2GFAV+nRIbRpWvyfiFLqBeALYLRSqlgpda7WOghcBrwHLASmaa0X5OBZonlkibnXrt5cGxt75Ztiht74doPXXnHQKCDSBGpQj3xOHjeI+07amR6BddzneYR9PzyGvOIvuCdwEvvV/wPGn09evnM/7PpgCJ9d87AJk7jwiG+yXfLcfHj1/uw7KpINrlP4PGLH1vtZzmXjMD9vn2E8fNquGc8XhI5Aa0RbnZpifCYwM8fPEp9HI1ldFg+5vf6ViIbh1MnPilnfqro+iFKKv03qC3P+hp7/FEca8HjoKNZseyHP/VAREzSFXudfQbMfuZVUmofVbGVfX6pQXSes98nmm/+Wo8ZkMTs1n994IBV1Ddf9EoS2QLt0mGeKUmoyMHnkyJGtvZScsbqshgVry5m0Q+rsbJOHZi2lf7d8TthtIADfrCpDaxg3tEcDV0KxRfOw82tpNbe9kawYdsmL+DbCteUw+y744l8QrEPt+n/s/9k4NtCDA+vzgArGDoqYEu15GSZ+J4d5qlDdNLYoeyvcdGYrq1xpDZNR/2759Cc/q2skVFdoLdqKz6NZ6Ig+j0c+Wc7V037IaO79Hyzh2unxude9/CN/f29xRteuTtNC9vIXvmXOko1J4108Qc5zvc1z1RfAnL/DNofBpfNg8oNsICKwNkWTCc2oqFSaR30wjM+TPtrK1CLsZdatpKptFTtv+dnVxkJ1BaEt06E1j47IkvWVSV/TmbC+vI4VG6sTWq2mo7gstfAwe22YuAhxgmsOB7x/NZM865mrx7LvhQ9B/52TrjVDe01BUODLXPOw9/qO1aeK7vR2Hwmkrm3lRGN9Hq1JO1mm0AHp0JpHR3OYa61ZsqGyUdd+vrwUgK9XRqKolpVUMvTGt/l8Wanj/HUVdUlCYu7SUra5+R1L2XXN4cZXvO+9nr95niBY2I9T/Lfwf/U3OAoOiJcxyW9A83AK1bWbkuKhupFjpw5/qWpbxe5pPddIn0drYpoKu+TJd6DQsnRo4dHRzFYllfVJ/b8z5bNlmxyP312w3mk6WsPaLYl+j3vfX4w/FGZdeR17Gz/xuvdWHvE+SAiD8/1Xs/7EN/ky7Ow8Nk1C9VGBlE7zCIc1wbBOEgZ2q5RZn8qIaR4O98rCYW60Q83jzAlD+PPkMfx+r6GtvRRhK0M+V9oRjdU6tNYxzQOSncipePqzXzlkTN/YcXV9kLFqGde7X2Jv1wKKdS+u8V/Eq+F9CGPwl7zUpcB7FnopqYyXMjEjpwocemr7owmKdmFg38/NXIyY8PAkfwvZ3zRdB9lMS7K3Jdwug7OkGKLQCnRo4dHRoq2WbKhq1HW/llazLtrpD2B9RV2a2XGe/WIVz36xCpehGKHWcHPNyxzo+4pNugt/DpzJ/0IH4ScuMDqnER7dCjwJwsM0W7ldBl63kWAiM39O1jxsPg9bZVwnn4c9z8M0SJmtWhNqW7WR8iSC0B4Qs1U7Ysn6xmkeny1PNFn9lkHZdJP+lHK361He917P+PCP3Bc4gcP0QzwTmpQgOCA5D8NKt/xER701cspelj2V8Nh7ZC9+P2FI7DmmpmD6Kpx8HnaBZkZkmc+3rrmtVNUVhPZAh9Y8OhpLShrpLF9WSv+ueayNah+r00RSmfSggkvcr/N/rg8AxVOhw/lX8Gg204UClwsIJV2TLjfCLI5oYtUSCrxuNtfEk+NMv4gvKdrK4PZjdmDO0lJ+La2O9R8393wnn8fLF03g8he+i1X4NQXEKbsPprSynosnxrVS0TwEIXM6tObRkdBas6yRZqsvVmxir5G9YsfpcjhcgSqucL3CnLyrONv1Lq+F9uGgwP38JXgGm8m8yqud3p19sZ/zPEaCoCm0Oc1jwiOFJmNeaobWpgvVHV7UibMszmRPVCB53QZXHzo6IUmxPTrMBaG16NDCoyOF6q4rr6OyPkgnX/bK4paaAHuP7Bk7/s1B8/AS4Oi619nn7YO5yvMK37l34VD/PdwQvIB19Eqab/J/ew7JaA03HL5t7Gd7sl+BLVzXH4xoNfY8DxNzYzdsZisn4QH28iW5L08iCFsjHVp4dCSfhxlpNapPp0Zdv9cIi+ZhER5Kh+i2ZDqzfddwYc0TfB8YyNH1d3Kr73qW6wGO96rxx01WvTpFNIqTxkVKoNx8xHb89dgdk67pkudh8tj+QNxZbmLXPMxwZLupK7bm6J/mZm9GjzmZrSBRi/CkKZebGG0l4kMQ0iE+j3aCKTy26d2ZX6L2+0wZ2bsTfbrEq9f+VlYLWnOY8TWXLrqJ3nUr+UEP507jYt6r3Q6AfXsUZNSPvKizj5V3Hxk7Pn+/4SnndvKZTur0moeJ2RTJjt2kVB/VVFKZuYwMNQ8lPg9ByBgRHu2EJRuqKOrso1th6nDYVOw9omfC8aiab5n05W383vszG/Vgfjv4UY55qzNWY83A7skbd4HXlaB1QER4ZIqZTZ4sPJI1Breh6NfVuUhgrEtgNLq3PmDmhaQyc1num85sJT4PQcgYER7thKUbKtmmsSarqLN8B7WC690vsZ/rJ7bU9ea6wAUU7nw6JwwbBnyacM3A7skbd8TZ3AThEfXX2ENqnTSPgd3zU1bLNTUEM1fDTCp0CtWFbMxW4vMQhEzp0D6PjkI4rFlaUsWo3p0bdf1eXctg2u95y3cLOxi/cmfgdK7u8xTTQxMJq8SNe3y0XLuTUKiuTy6N0qtTZoUWgZizP2zLcLfneQAMSmGygvjGbub/xTWPhn0eRpqa7IbkeQhCxnRozaMtZZibnfis/gEnfllbwREPzWXC8J68cMGeAKzZUkuNP8Q2fTrz0Kyl1AfDjPnTuwnXmeYk8/5aa/qyiT+4Z9D5yTngzuPB4HE8ETyCKgoYtiWuQViTsAf1KGDeyjLHL++gQ1kT02GeCabmYb9PgUMEWTrh0aeLj1/WxTWNvKjwSaUFpevhkWqeks8qQUhLhxYe7bGT4PRvVgOR3AwT01k+um+nWGmR0/cYnHDdE3N/jR/UlLFsxh184nseRRjGnw/7Xss/7poHRL70reG6a8sjBRAfPm1XZi8qSbhv9wJP2mKMdv9FOsyoqnSaR11Ui+jf1bk9LcA/Tt6Z9xasZ5s+EU3s1N0HEQqFOT1F2HC6ZlEJ85SYrQQhUzq08OgomDWtRlrMVjcfmVi99om5v1JAHcz5O/qzBxlRX8WM8L48EDyeTw8/K2HuoB4FLLKUOvlsWSn5HheHjOmTIDzGDuxK90Ivc5fGiyqOKCpk+cZqIFKvKhvMNrXBcGKpd6vmsTFa/6p359TCo1uBl5N3jwvPhooDZmqCErOVIGSOCI92wNINlfTtkkfX/BSbddDPma73uNz9GswuZ2n3/bis4kiW6EGO0+3C49NlpewxvEdKh7OVbft2iQmPbExWEI+2slu/rJqHKTyKumR373Sk83NYcUmoriBkjFh22wFLSiqdkwPDIfjhJZg6jjs8z7Jc96f4+Dc4suRidtp1Qsr7DbKE4a7dEukwuM/I1FnkVkb3jWs/RdkKj6iGEQjZNA+L8DAjp3pnEcXVEJn7PETzEIRMEeHRxgmFNUs3VMXs+xE0LH4HHt0HXr0A8rpypv8GTvHfwk1f55PndnHDpG1T3nNwj3gYrtnnY+/GCI8sN3hTs/HYyo6Y97GWXklntsoWV8Zmq5w9UhA6PO3yn4tSaoxSappS6hGl1AmtvZ7mZHVZDfXBMKOjwmO8WsjL3tvhhVMgWA8nPA0XfMKc8FhAMWfJRq48ZJu0G7s1kqnGH6JXJ2/s/g0xoFtc8GQrPAIh52S+XQd3Z+Yf9mW7fvHCiz0y7LWeCZmWGpGS7IKQOS0uPJRSTymlSpRSP9vGJymlFiullimlbmzgNocD/9RaXwyc2WyLbQXs3QLN453cq+D5E5jmu5OBaiMc9QBc+hXscFzCJ/Oo3p04c0L6YoX2sh97jeiVsV+g0KIdZOvzMCOp8m15HUopxvRPrNibaYRUJmR6LynJLgiZ0xoO82eAqcB/zAGllAt4GDgEKAa+Vkq9AbiAKbbrzwGeA25TSh0N9KQDsW5LYpe/klULedAzlW1f/xzyuvHXwKk8GzqMxeOOdbz+9qO3TzIL2bGXHsnU3wGJRQyz1TzqA5HckrwUyXzNhfg8BCH3tLjw0FrPUUoNtQ2PB5ZprVcAKKVeBI7RWk8Bjkpxq0ujQmdGsy22FenNZnjrKk6d/yx+lxv2vQb2+gOP3/5Z2uv2ykAQ2L/89x6VhfDwunEbimBYZy08ak3hkabj4NTTdsm4x3qmZBxtlWC2yukSBKHD0VZCdQcAqy3HxcAeqSZHhc9NQCHw9xRzLgAuABg8eLDTlDZJJ13F9e6XONv1Lnwb5m3vJGb1PpMHD5qU0+eYAgAS/RgNUeB1MaZ/F34sLk9Z9TYVZl7IDgNSl8g/aqf+Wd0zEzLO87BmmIvmIQhpaSvCIyu01iuJCoY0cx5XSq0DJnu93t1aZGFNwV8DXz3Kf6vvpdBVw+vhvZh8+VSu/cdizu6fWcOlbBjUo4BfS6uzFhxKKZ47Zw+mzV/NsF6FWT1ztyE9ePGCPdk9Wj+rpchUixCBIQiZ01airdYA1oy2gdGxJtEumkGFAvD1k/DQLjDrdn52jeEI/xSuClzKKt0bfyjMqAwjobJhTDSy6YBtizK+xnSWdy3wpO3bkY49h/fMqTM8E8R/IQi5p60Ij6+BUUqpYUopL3AK8EZTb9qm29CGw/DTyzB1d3j7augxDM5+l5vzb2GRjpjZlpoNoBpZij0d+4+OCI0Dt+2d8TVO1W/bAyI8BCH3tEao7gvAF8BopVSxUupcrXUQuAx4D1gITNNaL2jqs9qk5qE1LP0AHtsPXjkXvIVw2nQ4+x0YkpgVHq9plXvhYZqcvK5sChu2Sytni2s6grA10BrRVqemGJ8JzMzls9pSSXaA3dRieHoq/PY5dB8Kx/0bdjg+ZWrz4g2VDOqRn7JNa1MYN6Q7D5+2K3sOz9z/0F6Fh8gOQcg97XM3yJA2U5J9wwKe8NzLIa5voawPHHkf7HImuNNnUS/dUJlx5ne2KKU4cqd+WV3Tqb0KD5EegpBz2orPo1lodZ/H5pUw4wJ4ZG/2MBZxT+Bk+MN3sPt5DQoOgF9Lq5vFWd5Y2q/mIcJDEHJN+9wNMqTVNI+qEpjzd5j/NBgu2PsK9p01hnI6cb038/DWQEg3i7O8sbRXh3mmhREFQcicDi08WtznUVcOnz0EXz4CwTrY9UzY/3ro0p/yWW836paN7VveHLRXzcOUHWK9EoTc0T53gwxpMc0jUAvznoBP74fazREn+AE3Q88RTbqtoZon0qqxtFfhYUZbiflKEHJH+9wN2gqhIHz/PHz8N6hcCyMPhoP+BP3GNvqW1qpOQ3oWZtUjvCHm33Iwugllozr52qfZyhQa4jgXhNwhwqMxhMOw8HWYfRdsWgYDx8PxT8DQfXL6mFE51jqyLaFupzlChlsCMxJaZIcg5I72uRtkSM59HlrD8tkw6w5Y9z0UbQenvACjD2+WBhDbtKFIK2jHobrRvxtxnAtC7ujQobo5zTAvng/PTobnj4PaMjj2Mbj4M9j2iGbrHOTYt7wVaa8+DzFbCULuaZ+7QUtSsghm3wmL3oKCXnD4PbDbWeBumgkoE6z9wtsChe3U52GKDHGYC0LuEOGRjlWfwzNHgqcQDrgF9rwYfC2nDWRb8ry5KWynPo9wNEpAFA9ByB3tczfIkCb7PAbtAQfeArueBYUt3+3W18LtWhuivZqtPO6IdTZdEypBELKjfe4GGdLkPA/DFWn/KgDt12HeJc/DtAsnsF2/tmUGFIT2TPvcDTowddE+322R9urzABg/rGW7FwpCR6dDR1u1RzZU1LX2EpIYHvW95OcwYVEQhPaNCI82RhMSwJuNo3bqD0i0kiAIcTq08Gj1kuyCIAgdlA4tPNpkG1pBEIQOQIcWHoIgCELzIMKjjdGUqreCIAgthQgPQRAEIWtEeAiCIAhZI8JDEARByJo2LzyUUsOVUk8qpV62jBUqpZ5VSj2hlDq9NdcnCIKwNdKswkMp9ZRSqkQp9bNtfJJSarFSaplS6sZ099Bar9Ban2sbPg54WWt9PnB0jpctbMVcd9hoTh0/uLWXIQhtnuaubfUMMBX4jzmglHIBDwOHAMXA10qpNwAXMMV2/Tla6xKH+w4Efor+3HaLQQntjksPyFHXSUHo4DSr8NBaz1FKDbUNjweWaa1XACilXgSO0VpPAY7K8NbFRATI96TQnpRSFwAXAAweLF+SgiAIuaQ1fB4DgNWW4+LomCNKqZ5KqUeBXZRSf4wOzwCOV0o9ArzpdJ3W+nGt9Tit9biioqIcLV0QBEGAdlCSXWu9CbjINlYNnN3QtU1uBiUIgiA40hqaxxpgkOV4YHRMEARBaCe0hvD4GhillBqmlPICpwBvNMeDpDCiIAhC89DcobovAF8Ao5VSxUqpc7XWQeAy4D1gITBNa72gmZ4vJdkFQRCageaOtjo1xfhMYGZzPjv6nKb1MBcEQRAcafMZ5oIgCELbQ+kOXAPcjLYCTgaWAl0Bqw0r3bH1515AaQ6WZH9eY+emOuc03ph3ztX7plpTY+bl6p1Tndta3rk5fq9Trakx89rLO+fq33K68w3927WPNcc7D9FaO+c6aK23mv+AxzM9tv08vzme39i5qc45jTfmnXP1vtm8c0PzcvXOqc5tLe/cHL/XW+M75+rfcjbv3ND/g5b4e7b+t7WZrewJhemOHZMPc/z8xs5Ndc5pvL28c0PzcvXODf3/yAVt+Z2b432zuW9Heedc/VtOdz6T39WW/nuO0aHNVrlCKTVfaz2utdfRUmxt7wvyzlsL8s65Y2vTPBrL4629gBZma3tfkHfeWpB3zhGieQiCIAhZI5qHIAiCkDUiPARBEISsEeEhCIIgZI0IjyzZGvunO/WR7+gopX4X/Tt+SSl1aGuvpyVQSm2nlHpUKfWyUuri1l5PSxH9Nz1fKZVpM7p2jVJqolJqbvTvemJj7yPCg6x7rXeI/unZvLN27iPf7sjynV+L/h1fRKRCQbsky3deqLW+CDgJ2Ls11psLsvz3DHADMK1lV5lbsnxnDVQBeUSa8TWO5sg8bG//AfsBuwI/W8ZcwHJgOOAFfgDGAH8Edo7O+V9rr70l3tly/uXWXncrvPN9wK6tvfaWemciH0TvAKe19tpb4p2BQ4i0hTgLOKq1195C72xEz/cB/tvYZ4rmQaTXOlBmG471Wtda+4EXgWOI90+Hdqy5ZfnOHYJs3llF+Bvwjtb625Zea67I9u9Za/2G1vpwoN2aZLN854nAnsBpwPlKqXb5bzqbd9Zah6PnNwO+xj6zzbehbUWceq3vATwETFVKHUkLlABoYRzfWSnVE/gL0T7yWusprbK65iHV3/PlwMFAV6XUSK31o62xuGYi1d/zRCJmWR8t0DKhhXF8Z631ZQBKqbOAUsvG2hFI9fd8HHAY0A2Y2tibi/DIEp1h//SOhHboI9/R0Vo/RORDYatBa/0x8HErL6NV0Fo/09praCm01jOAGU29T7tU0VqIrbHXuryzvHNHRd45x+8swiM1LdZrvQ0h7yzv3FGRd87xO4vwoPV7rbcG8s7yzsg7yzs35ZnRkC1BEARByBjRPARBEISsEeEhCIIgZI0ID0EQBCFrRHgIgiAIWSPCQxAEQcgaER6CIAhC1ojwEIRWQCmllVL3WY6vVUr9uRWXJAhZIcJDEFqHeuA4pVSv1l6IIDQGER6C0DoEgceBq1p7IYLQGER4CELr8TBwulKqa2svRBCyRYSHILQSWusK4D/AH1p7LYKQLSI8BKF1eQA4Fyhs5XUIQlaI8BCEVkRrXQZMIyJABKHdIMJDEFqf+wCJuhLaFVKSXRAEQcga0TwEQRCErBHhIQiCIGSNCA9BEAQha0R4CIIgCFkjwkMQBEHIGhEegiAIQtaI8BAEQRCyRoSHIAiCkDX/D/uV56Ks+jWOAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.xscale(\"log\")\n", "plt.yscale(\"log\")\n", "plt.xlabel(\"N\")\n", "plt.ylabel(\"abs(error)\")\n", "plt.plot(fN,error)\n", "plt.plot(fN,error_fit);" ] }, { "cell_type": "markdown", "id": "9dbec280", "metadata": {}, "source": [ "Cautionary tale: When I first did this in IDL, there were some better features, and some worse ones.\n", "\n", "**Better:** For boost=1, x always stayed equal to 1 at machine precision. In fact, even for boost=1e9, x only began to deviate from 1 for $N>9$.\n", "\n", "**Worse:** For boost=1e9 and larger values of $N$, the error increased as $N^{1.3}$." ] }, { "cell_type": "code", "execution_count": null, "id": "4d0aec4e", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.2" } }, "nbformat": 4, "nbformat_minor": 5 }