{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ "ESE - Numerical Methods I: Basics of Computer Arithmetic" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Floating point (FP) arithemtic $\\neq$ arithmetic in ${\\rm I\\!R}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The motivation for this session can be summed up by the following problem:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "1.2-1. == 10.2-10." ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 1, "text": [ "False" ] } ], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "or, actually the same, but expressed in a more challenging manner:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "1.2-1. != 10.2-10." ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 2, "text": [ "True" ] } ], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Is that a bug in Python? " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That would be a lot of people using a buggy language. Let's try this in Matlab:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bummer (remember the binary convention `0 == False` and `1 == True`). Apparently it has nothing to do with whether we use Python, Matlab, C++, etc. Let's try some reverse engineering to get closer to the root of the problem:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "0.2-0." ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 3, "text": [ "0.2" ] } ], "prompt_number": 3 }, { "cell_type": "code", "collapsed": false, "input": [ "1.2-1." ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 4, "text": [ "0.19999999999999996" ] } ], "prompt_number": 4 }, { "cell_type": "code", "collapsed": false, "input": [ "10.2-10." ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 5, "text": [ "0.1999999999999993" ] } ], "prompt_number": 5 }, { "cell_type": "code", "collapsed": false, "input": [ "100.2-100." ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 6, "text": [ "0.20000000000000284" ] } ], "prompt_number": 6 }, { "cell_type": "code", "collapsed": false, "input": [ "1000.2-1000." ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 7, "text": [ "0.20000000000004547" ] } ], "prompt_number": 7 }, { "cell_type": "code", "collapsed": false, "input": [ "10000.2-10000." ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 8, "text": [ "0.2000000000007276" ] } ], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The problem seems even trickier than we had thought. First, the answer is never $0.2$ exactly (besides the trivial case $0.2-0.$, which really triggers not much at all), but it seems to deviate more and more with larger numbers. It must be concluded that numerical calculations on binary systems like computers contain errors, errors that do not exist for calculus in ${\\rm I\\!R}$. Before we explore the nature and source of these errors, let's quantify them in terms of a 'correct' solution:\n", "\n", "First, the absolute error:\n", "\n", "\\begin{equation}\n", "\\epsilon = |\\tilde{x}-x|\n", "\\end{equation}\n", "\n", "where $\\tilde{x}$ is the computed solution and $x$ is the exact solution.\n", "\n", "For theabove problem, let's see how the absolute error scales with the magnitude of the numerical values involved." ] }, { "cell_type": "code", "collapsed": false, "input": [ "%pylab inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "prompt_number": 9 }, { "cell_type": "code", "collapsed": false, "input": [ "# remember that in numpy, unless otherwise specified by the user, numbers are treated as floats\n", "# it is good practice though to indicate what numbers you consider as floats by using the decimal point\n", "ls = np.logspace(1.,25,25)\n", "x = 0.2\n", "# operators act elementwise on array\n", "x_tilde = (ls+x)-ls\n", "absolute_error = np.abs(x_tilde-x)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 10 }, { "cell_type": "code", "collapsed": false, "input": [ "plt.plot(ls, absolute_error)\n", "plt.xscale('log')\n", "plt.yscale('log')\n", "plt.xlabel('$O(x)$')\n", "plt.ylabel('$O(\\epsilon)$')\n", "plt.title('Absolute Error')\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEkCAYAAAD5BJxYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1cVHXe//EXmnlTrUYoWblSXqSQmGs/gjBsLCKTzbvN\nkowMMlFaJd20m6urW9tVa1eKvMFNSc1qddfUslDIRjQDcVNTktBy3F2xlNzwJm9Sz++Pc+1cEqgg\nc+bMDO/n48GjnTMz5/OZ9sSb853v+Z4gwzAMRERELNTE7gZERCTwKWxERMRyChsREbGcwkZERCyn\nsBEREcspbERExHIKG2m0HnzwQf7nf/7Ho/t88803iY+P9+g+RQKBwkYCnsPhIDg4mOPHj1fbHhQU\nRFBQkE1dNTyY3nzzTZo2bcoll1zi/vnFL37Bt99+68EuRTxDYSMBzeVysX79etq1a8eyZcvsbsfj\nevbsycGDB90/Bw4c4PLLL6/xuhMnTtRp29nU9/Uip1PYSECbN28eCQkJpKSkMHfu3BrPHzx4kEGD\nBnH55Zfz+OOP8/3337ufmzRpEt27d6d169Z069aN0tJSAH788UdmzJhBVFQUiYmJvP/++7XWdrlc\nNGnShFOnTrm3ORwOZs+eTVlZGSNHjuSzzz7jkksuITg4GDB/oS9cuJBbb72V7t27M3v27BpnZKc7\n2wIgYWFhTJ8+nbi4ONq0acPXX39NkyZNWLRoEV27duX2228HYNmyZdx+++1ERUUxc+ZMfvzxx2r9\n//z1IufFEAlgnTp1Mt566y2jvLzcaNasmfHdd9+5nxs2bJhx0UUXGXPnzjV2795tDBkyxBgyZIhh\nGIaxdetWo3Pnzsbu3bsNwzCMsrIyY8+ePYZhGMYzzzxj9O7d2/jqq6+Mjz/+2AgLCzM++eQTwzAM\nIzc317j55psNwzCMnTt3GkFBQcbJkyfdNR0OhzF79mzDMAzjzTffdL/2P1599VXj1ltvNbZu3Wrs\n2LHDcDgcxqxZs2r9bKfXqk1YWJgRGRlpFBYWGkePHnX3M2DAAOPrr782jh49aqxatcr45S9/aeTn\n5xvl5eXGbbfdZjz77LPV+j/99SLnS2c2ErDWrl3L7t276devH+Hh4URGRvL2229Xe02PHj144IEH\nuOKKK3j++edZsWIFp06d4uTJkxw9epTt27dz6tQpOnfu7B6eWrp0KU888QTXXnstt956K0OHDuW9\n996rd39GLWclCxcu5MUXX+S6666jU6dOZGZmsmTJkjPuo6ioiEsvvdT9Ex4eXu35IUOGEB8fT/Pm\nzd3bxo0bxzXXXEPz5s1ZsmQJQ4cOJSEhgfDwcJ544okan+X014ucL4WNBKy5c+eSmJjIJZdcAsDg\nwYOrDaUFBQVx/fXXux9fe+21/PTTT2zbto1u3brx0ksv8cQTT3DllVfyzDPP8OOPP3Lw4EG++OIL\nbrjhBvf7brjhBtasWdPgfg8fPsy6detISkpyh8eDDz7IunXrzvie2NhY/v3vf7t/tm/fXu35mJiY\nGu85fdu6detqfJYtW7Zw8ODBs+5DpL4UNhKQjhw5wsKFC1m1ahXt27enffv2/PGPf2Tz5s188cUX\ngHlmsWnTJvd7vvrqK5o1a0ZERAQAQ4cO5bPPPqOoqIiVK1eSm5vLJZdcQrdu3diwYYP7fRs2bKBX\nr141emjbti3NmjVzzw47ceIEW7ZscT/ftGnTamc3F110ETExMaxYscIdHj/88AP//ve/z/vfwwUX\nXHDWbT179qzxWaKiotwBfaZ9iNSXwkYC0pIlS7jgggvYtm0bmzdvZvPmzWzbto34+HjmzZvnft3G\njRtZsGABFRUVvPDCC/Tp04cmTZqwYcMGiouL+emnn2jZsiUXXHCB+xdw//79efnllykvL8fpdPLO\nO+8wYMCAGj1cdNFFxMbGMmvWLPbv388f/vCHamcMN9xwA9u3b+fQoUPubSkpKTzzzDN8/vnnnDp1\nit27d7Ny5UrL/j3179+fd955h1WrVrFjxw5efvllBg4caFk9abwUNhKQ5s2bR1paGldddRXt2rWj\nXbt2hIaG8tvf/pa3336bkydPEhQUxIgRI/jb3/5Gjx49uPLKK3nttdcAOHDgACNGjCA4OJjevXtz\n4403cv/99wMwYcIEBgwYwKBBg3jppZf405/+xC233ALUvHZn0qRJfPbZZ0RFRXHq1Cl69uzpfi4y\nMpIBAwZw3XXX0a5dOwAefvhh0tLSeOaZZwgODub222+nvLy81s8YFBTkns12+s/f//73M/57+fl1\nRQ6Hg6lTp/L73/+eAQMG0L9/f8aPH3/G14ucryCjtm8p/cixY8d48sknOXLkCP3796dPnz52tyQi\nIj/j92c2n376KdHR0cyYMYPFixfb3Y6IiNTCJ8MmLS2N0NBQoqKiqm0vLCwkIiKC8PBwsrOzAdiy\nZQudOnUCzC+FRUTE9/hk2KSmppKXl1dje2ZmJjk5ORQUFDBt2jQqKyvp1q0b33zzDQCtWrXydqsi\nIlIHPhk28fHxXHrppdW2VVVVAdCrVy86duxIYmIixcXFxMXFsWHDBkaPHs2gQYPsaFdERM7BbybQ\nl5SU0KVLF/fjyMhIioqKSEpK4pVXXjnn+0NCQqqteyUiIud22WWXUVlZ2eD9+OSZjRW+//57DMPw\n2s8tt9yien5aL5A/m+qpXn1/PPVHut+ETXR0NGVlZe7HpaWlxMbG2tjR2YWFhamen9YL5M+meqpn\nF78Jm9atWwPmjDSXy0V+fr5Pr9kU6AdgINcL5M+meqpnF58Mm+TkZOLi4igvL6dDhw7k5uYCkJWV\nRXp6OgkJCWRkZBASEmJzp2fmcDhUz0/rBfJnUz3Vs4vfryBQV0FBQTSSjyoi4jGe+t3pk2c2IiIS\nWBQ2IiJiOYWNiIhYzm8u6hQRe+zaBWPHwnnc+VrETWc2IlKrY8fg97+HG26AHj3gyBEwDP00th9P\n0ZmNiNSwYgWMHg2RkbBhA/jppR3iQxQ2IuL2j3/AuHGwcSO89hokJdndkQQKDaOJCMePw6RJ5nBZ\nVBSUlipoxLN0ZiPSyBUUwG9/C//1X7B+PVxzjd0dSSBS2Ig0Uv/6F/zud2bAvPYa3HWX3R1JINMw\nmkgjc/w4vPwydO8OnTubQ2YKGrGazmxEGpGyMhg0CDp2hKIic+hMxBu0EKdII/H99xATA489Bunp\nEBRkd0fiDzz1u1NhI9II/PQTJCbC//t/5hCaSF0pbOpJYSONlWHAqFGwezcsWQJNm9rdkfgTT/3u\n1Hc2IgHu9ddh7VpYt05BI/bRmY1IAFu5Eh54wAwaXT8j50NnNqdZunQpy5cv58SJE4wcOZIbb7zR\n7pZEbPfVV5CSAosWKWjEfgF1ZrN3716effZZZsyYUeM5ndlIY7J/P8TGwuOPw0MP2d2N+LOAvC10\nWloaoaGhREVFVdteWFhIREQE4eHhZGdnn/H9kydPJj093eo2RXzaTz/BPfeYa5spaMRX+FTYpKam\nkpeXV2N7ZmYmOTk5FBQUMG3aNCorK5k/fz5jx46loqICwzCYMGECffv2pXv37jZ0LuI7xo6FZs00\nxVl8i099ZxMfH4/L5aq2raqqCoBevXoBkJiYSHFxMSkpKaSkpADw2muvsWrVKg4ePMiOHTt0diON\n1owZ8PHH5uoAF/jUf93S2Pn84VhSUkKXLl3cjyMjIykqKiLptPXPx4wZw5gxY865L4fDQVhYGGFh\nYTgcDhwOhxUti9hi1Sp4/nlzmnPr1nZ3I/7K6XTidDpxuVw1/vhvCJ8PG09yOp12tyBiie3bITkZ\n3n1X651Jw/z8D/EgD61r5FPf2dQmOjqasrIy9+PS0lJiY2Nt7EjEt/zwA/TrZ57V9O5tdzcitfP5\nsGn9v+MBhYWFuFwu8vPziYmJsbkrEd9w4gQMGQIJCTBypN3diJyZT4VNcnIycXFxlJeX06FDB3Jz\ncwHIysoiPT2dhIQEMjIyCAkJsblTEd8wfjycPAlTp9rdicjZBdRFnWejizol0Lzxhjm9uagILr3U\n7m4kUGm5GpFG6sQJmD4dJk40Z54paMQfKGxE/Minn8Ijj0BwMKxeDddea3dHInWjsBHxA3v3muuc\n5efDK6/AvffqTpviX3xqgoCIVHfyJEybBtddB5ddBtu2mbPPFDTib3RmI+KjioogIwN+8Qv45BPo\n2tXujkTOn8JGxMfs2wdPPgkffmjONrvvPp3JiP/TMJqIjzh5EmbONIfMLrnEHDIbOlRBI4FBZzYi\nPmD9enOWWYsWUFAA3brZ3ZGIZ+nMRsRGP/4II0ZA//4wZgwUFipoJDDpzEbEJqdOwbBh5v/etg3a\ntLG3HxErKWxEbPLCC7B7t3kfmhYt7O5GxFoKGxEbLFwIubnmdzUKGmkMtBCniJdt2AB33mmuBtC9\nu93diJydp353aoKAiBdVVMDAgTBrloJGGheFjYiXHDkCAwaYNzkbONDubkS8S8NoIl5gGP+3EsCC\nBbpQU/yHhtF+5vDhw0RHR7N8+XK7WxGp4aWX4OuvYfZsBY00TgEzG23KlCnce++9drchUsPixZCT\nY848a9nS7m5E7OFTZzZpaWmEhoYSFRVVbXthYSERERGEh4eTnZ1d4335+flERkbStm1bb7UqUicb\nN0J6OixZAu3b292NiH186jubNWvWcPHFF/PAAw+wZcsW9/Zf/epXvPrqq3Ts2JE77riDtWvX8tFH\nH/H5558zfvx4pk+fzuHDh/nyyy9p2bIl7733HkE/G6vQdzbibd9+CzEx5s3OBg+2uxuR8+Op350+\nNYwWHx+Py+Wqtq2qqgqAXr16AZCYmEhxcTEpKSmkpKQAMHHiRADmzp1L27ZtawSNiLcdPWrOOEtN\nVdCIgI8No9WmpKSELl26uB9HRkZSVFRU62uHDRtG3759vdWaSK0MAx5+GDp0gGeesbsbEd/gU2c2\nVnM4HISFhREWFobD4cDhcNjdkgSgKVPgyy9hzRpo4vN/zolU53Q6cTqduFyuGiNNDeHzYRMdHc34\n8ePdj0tLS+nTp8957cvpdHqoK5HaLVsG2dnmLZ1btbK7G5H6+/kf4p76WsLn/+5q3bo1YM5Ic7lc\n5OfnExMTY3NXIjV98QU89JA51fmqq+zuRsS3+FTYJCcnExcXR3l5OR06dCA3NxeArKws0tPTSUhI\nICMjg5CQEJs7Falu717zBmivvgo33mh3NyK+x6emPltJU5/FKseOwW23wS23mCsFiAQST/3uVNiI\nNIBhmENnP/wAf/2rJgRI4AnI62xE/M2f/gSffw5r1ypoRM5GYSNynj78EP74R3Pm2cUX292NiG9T\n2Iich9JSePBBWLoUfvlLu7sR8X068Repp8pK6NfPXPPsppvs7kbEP2iCgEg9HD8OiYnmApuTJ9vd\njYj1NButnhQ20lCGYd4u4Ntv4b33oGlTuzsSsZ5mo4l4WXY2fPYZrFunoBGpL4WNSB2sWAF/+IMZ\nNpdcYnc3Iv5HYSNyDmVlkJICf/sbhIXZ3Y2If9JsNJGz2L8f7roLJk2C+Hi7uxHxX5ogIHIGP/0E\nffpA9+7mxZsijZFmo9WTwkbqKyMDXC54/31NCJDGS7PRRCw0fTo4neaEAAWNSMPpzEbkZwoK4P77\n4dNPoVMnu7sRsZenfndqgoDI/zIMePdduO8+858KGhHP0TCaCPDll/Db35qzz5Ysgbg4uzsSCSwB\nc2bz+uuvM27cOObNm2d3K+JHDh2CCRPMu2wOHAgbNihoRKwQEGGzceNGVqxYQdOmTYmIiLC7HfED\nhgELF0JEBHz3HWzdCqNHwwU61xexhE+FTVpaGqGhoURFRVXbXlhYSEREBOHh4WRnZ9d439q1a+nd\nuzdTpkxh+vTp3mpX/FRZmbly88SJ8PbbMHcuhIba3ZVIYPOpsElNTSUvL6/G9szMTHJycigoKGDa\ntGlUVlYyf/58xo4dS0VFBd26dSM4OJigoCBOnjxpQ+fiDw4fhiefhJtvhqQk83bOWhVAxDt8atAg\nPj4el8tVbVtVVRUAvXr1AiAxMZHi4mJSUlJISUkBICQkhPz8fMaNG0dSUpJXexbfZxiweDGMHQu9\nesGWLdC+vd1diTQuPhU2tSkpKaFLly7ux5GRkRQVFVULlQsvvJCJEyeec18Oh4OwsDDCwsJwOBw4\nHA4rWhYfUl5ufhezezfMn29OBBCRM3M6nTidTlwuV40//hvC58PGk5xOp90tiJecOgXPPWeuBPDU\nU2bgNGtmd1civu/nf4gHBQV5ZL8+9Z1NbaKjoykrK3M/Li0tJTY21saOxB88/TR8/DFs3gzjxilo\nROzm82HTunVrwJyR5nK5yM/PJyYmxuauxJe99Za5AsDSpXDllXZ3IyLgY2GTnJxMXFwc5eXldOjQ\ngdzcXACysrJIT08nISGBjIwMQkJCbO5UfFVRkXkms2wZ6DAR8R1aiFMCxj//CbGxkJMDv/613d2I\nBAYtxClymsOHoX9/ePRRBY2IL9KZjfi9U6fgnnvg4oshNxc8NHlGRNDN00TcnnsO9uyBVasUNCK+\nSmEjfu0vf4F586C4GJo3t7sbETkTDaOJ3yopgb59zTtrXn+93d2IBCZNEJBGbfdu8/4zf/6zgkbE\nHyhsxO/8+CMMGACPPGL+U0R8n4bRxK8YBiQnmzc5mz9fEwJErKbZaNIoTZwIO3fC6tUKGhF/orAR\nv/G3v5nf0RQXQ4sWdncjIvWhYTTxCxs3mrdyXrECevSwuxuRxkOz0aTR2LPHXIpmxgwFjYi/UtiI\nTzt61JziPHw43H233d2IyPnSMJr4LMOAlBT46Sfz/jSaECDifZqNJgFv0iT46ivNPBMJBHUOm+PH\nj/P+++/z/vvvc+LECZo2bcrBgwcJDg4mMTGRu+++myZNNConnrFkCUybZs48a9XK7m5EpKHqNIz2\n6aefsnz5cu6//37Cw8NpdtoN3Y8cOcLWrVt56623GDZsGD189BtcDaP5j82bISEBPvwQoqPt7kak\ncfPU785zhs2xY8fYvn07Xbt2PefOSkpKiLbht0NVVRWPPvoobdq0ITIykocffrjGaxQ2/uG77yAm\nxhxCGzLE7m5ExGtTn5s3b14taD788EMAdu3aBUB5ebn7OTuCBqC4uJibbrqJqVOnUlBQYEsP0nDH\njsGgQfDAAwoakUBT7y9Z1q9fD4DT6QTgX//6l8eaSUtLIzQ0lKioqGrbCwsLiYiIIDw8nOzs7Brv\ni42NZcGCBdx2223ceeedHutHvMcwID0d2rc3b4YmIoGl3mFz/fXXc+utt7Jy5UqWLVvGli1bPNZM\namoqeXl5NbZnZmaSk5NDQUEB06ZNo7Kykvnz5zN27FgqKir4y1/+whNPPMHHH3/MBx984LF+xHte\neQW++ALmzgXNMxEJPOd1nU15eTl5eXm0adOGe++9l+YevEWiy+XirrvucodYVVUVDoeDjRs3AjBm\nzBjuuOMOkpKS3O/ZtWsXzz//PJdddhlt27ZlwoQJNfar72x81wcfmGc1RUXQoYPd3YjI6bx2nc2x\nY8c4ePAgISEh7m3XXnst1157bY3Xfv3113Tq1KnBTZ2upKSELl26uB9HRkZSVFRULWw6duzInDlz\nzrkvh8NBWFgYYWFhOBwOHA6HR3uV+tu6FdLSYNkyBY2IL3A6nTidTlwuFy6Xy2P7PWfYNG/enPz8\nfA4cOMDAgQNp2bJljdfs27ePrKwsHA6Hx8PGk/7zPZP4hn37oF8/mDoVYmPt7kZEgBp/iAd56Irq\nOl3U+etf/5o9e/YwdepU9u7dy9GjRzl69CgHDhygRYsWdO/enQkTJtC6dWuPNHW66Ohoxo8f735c\nWlpKnz59PF5HvOv4cXOtsyFDYOhQu7sREavVeQWB5s2bc/3119OyZUt69OhBmzZtrOzL7T8BVlhY\nyC9/+Uvy8/N59tlnvVJbrGEYkJEBl15q3gxNRAJfneb9LF++nBEjRlBYWEhubi4xMTGMGTOGqqoq\njzaTnJxMXFwc5eXldOjQgdzcXACysrJIT08nISGBjIyMat8fif959VUoKYG33tLMM5HGok6z0V55\n5RUee+wx92PDMFi5ciUzZ85kwYIFtPKDxas0G803fPQRPPQQfPYZdOxodzcici5evXlaZGRkjeJ3\n3HEHs2bNYubMmQ1uQhqHbdtg2DBYtEhBI9LY1ClsDh06xNGjR2tsb9u2Le3atfN4UxJ4Pv4Y7rwT\nXn4Zeva0uxsR8bY6hc3NN99MRkYGmzZtqvHcF1984fGmJHDs3m3OOHvoIfO7mmHD7O5IROxQ5xUE\nqqqqGDp0KLt37yY2NpZWrVpRUlJCeno6Q/1g7qq+s/Gun34yw2XSJBg1Cp58UvelEfFHXrvFwM9V\nVFTwySefcOjQIRwOB507d25wE96gsPEepxMeecRcESA7G8LD7e5IRM6XbWHjrxQ21tuzBx57DNau\nNVcFGDhQt3MW8XdenY0mcjYnTpjhEhVlzjL78kvzvjQKGhH5jzqvICBSmzVrzCGzyy+HTz8FPxlV\nFREvU9jIefn2W5gwwfx+5k9/gt/8RmcyInJmGkaTeps50xwyu+IKc8js7rsVNCJydjqzkXqZPRv+\n+EcoLISICLu7ERF/odloUmdr1pjDZYWFcNr97EQkgGk2mnjVzp1wzz3mSs0KGhGpL4WNnNPBg+Yd\nNZ98EhIT7e5GRPyRhtHkrE6eNC/OvPxyyMnRRACRxsZTvzs1QUDO6r//Gw4cgL/+VUEjIudPYSNn\nNG+eee+Z4mK48EK7uxERf+Z339ns3LmT4cOHM3jwYACOHj3KuHHjGDVqFHl5eTZ3FzjWrTPXOVu2\nDHQXbhFpKL8Lm6uvvpo33njD/XjdunVER0czY8YMFi9ebGNngeMf/zAv1HzzTbjuOru7EZFAYFvY\npKWlERoaSlRUVLXthYWFREREEB4eTnZ29jn3s2XLFjp16gTAkSNHLOm1MTl0yJx59rvfQd++dncj\nIoHCtrBJTU2tddgrMzOTnJwcCgoKmDZtGpWVlcyfP5+xY8dSUVFR4/XdunXjm2++AaCV7s7VIKdO\nwQMPQI8eMG6c3d2ISCCxLWzi4+O59NJLq22rqqoCoFevXnTs2JHExESKi4tJSUlh6tSpXHHFFezf\nv5+RI0eyceNGJk+eTFxcHBs2bGD06NEMGjTIjo8SMJ55BvbtgxkzNPNMRDzLp2ajlZSU0OW0y9Mj\nIyMpKioiKSnJvS04OJiZM2dWe98rr7xSp/07HA7CwsIICwvD4XDgcDg80ncgePttWLDAnHnWvLnd\n3YiIXZxOJ06nE5fLhcvl8th+fSpsrOZ0Ou1uwSetXw+ZmbBqFbRrZ3c3ImKnn/8hHuShYQ6fmo0W\nHR1NWVmZ+3FpaSmxsbE2dhT4/vUv866as2ebtw0QEbGCT4VN69atAXNGmsvlIj8/n5iYGJu7Clw/\n/ggDBsDo0eYMNBERq9gWNsnJycTFxVFeXk6HDh3Izc0FICsri/T0dBISEsjIyCBEVxRa4tQpePBB\niIw077gpImIlLcTZSD3/POTlwSefQIsWdncjIr5KC3HKeVu0CObMMWeeKWhExBt0ZtPI/P3v0KcP\n5OdD9+52dyMivk536pR627PHvDdNTo6CRkS8S2HTSBw5Av37Q3q6OdVZRMSbNIzWCBgGDB1q/u8F\nC7QUjYjUnSYISJ39/vewYwesXq2gERF7KGwC3HvvwcyZ5syzli3t7kZEGisNowWwTZsgMRE++ghu\nuMHubkTEH2k2mpzVt9+aEwKmTVPQiIj9FDYB6OhRc4pzaioMHmx3NyIiGkYLOIYBw4aZgfPuu9BE\nf06ISANoNprU6uWXobQU1qxR0IiI71DYBJBly+C116CoCFq1srsbEZH/o7AJEFu2wPDh8MEHcNVV\ndncjIlKdBloCwL595s3PsrLgxhvt7kZEpCaFjZ/btAnuuAPuu8/8ERHxRQobP/XDDzBmjBk0o0bB\niy/a3ZGIyJn5Zdjs3LmT4cOHM/h/LyJZunQpI0aMIC0tjfXr19vcnbUMA+bNg4gIOH4cvvwSHn5Y\nM89ExLf59XU2gwcPZtGiRe7He/fu5dlnn2XGjBk1XhsI19l88QU88oh5Dc306RAdbXdHIhLoAmK5\nmrS0NEJDQ4mKiqq2vbCwkIiICMLDw8nOzq7z/iZPnkx6erqn27RdVRU8+igkJMD995tTmxU0IuJP\nbA2b1NRU8vLyamzPzMwkJyeHgoICpk2bRmVlJfPnz2fs2LFUVFTUeL1hGEyYMIG+ffvSPYBuQWkY\n8NZb5pDZ4cPmkFl6OjRtandnIiL1Y+t1NvHx8bhcrmrbqqqqAOjVqxcAiYmJFBcXk5KSQkpKCgD7\n9+/nqaeeYtOmTUyaNImLLrqIVatWcfDgQXbs2BEQZzdbt5pDZocOweLFEBtrd0ciIufP5y7qLCkp\noUuXLu7HkZGRFBUVkZSU5N4WHBzMzJkzq71v9OjR59y3w+EgLCyMsLAwHA4HDofDY317yoED8Nxz\nMH8+PP+8zmRExLucTidOpxOXy1XjZKAhfC5srOR0Ou1u4ayWLoWMDHM6c2kptGtnd0ci0tj8/A/x\nIA/d3tfnJsxGR0dTVlbmflxaWkpsIxhD+uQTGDECFi6EOXMUNCISWHwubFq3bg2YM9JcLhf5+fnE\nxMTY3JW1duyA5GR45x3o2dPubkREPM/WsElOTiYuLo7y8nI6dOhAbm4uAFlZWaSnp5OQkEBGRgYh\nISF2tmmpqipzXbPnnoNbb7W7GxERa/j1RZ314YsXdZ44AXfdBZ06weuv292NiEhNAXFRZ2M3YYIZ\nOFlZdnciImKtRjUbzZfMnm3ee6a4GC7Q/wsiEuA0jGaDwkIYPNj8Z+fOdncjInJmGkbzUzt3wj33\nmMvQKGhEpLFQ2HjRgQPmhID//m+4/Xa7uxER8R4No3nJyZPQvz9cdRXMmAEeuihXRMRSGkbzM08+\naa7cnJ2toBGRxkfzoLzgzTfNlZuLi6FZM7u7ERHxPg2jWezTT2HgQFi92rwvjYiIP9Ewmh/YtQvu\nvhvmzlVw5UaHAAAL9ElEQVTQiEjjprCxyKFD5ppnEybAnXfa3Y2IiL00jGaBU6dg0CAICYE//1kT\nAkTEf3nqd6cmCFjg6adh/37z3jQKGhERhY3HLVgA775rzjy78EK7uxER8Q0aRvOg4mJzhYBVq6Br\nV0tLiYh4hWaj+Zh//tP8nmb2bAWNiMjPKWw84PBhcymazEzzzEZERKrzu7DZuXMnw4cPZ/Dgwe5t\nhw8fJjo6muXLl3u9n1On4MEHISoKxo/3enkREb/gd2Fz9dVX88Ybb1TbNmXKFO69915b+nnhBdi9\nG3JyNPNMRORMbAubtLQ0QkNDiYqKqra9sLCQiIgIwsPDyc7OPud+8vPziYyMpG3btla1ekYLF0Ju\nrrnuWYsWXi8vIuI3bAub1NRU8vLyamzPzMwkJyeHgoICpk2bRmVlJfPnz2fs2LFUVFTUeP3q1asp\nKiri7bff5s9//rPXLtz8+9/hkUdgyRK4/HKvlBQR8Vu2XWcTHx+Py+Wqtq2qqgqAXr16AZCYmEhx\ncTEpKSmkpKQAsH//fp566ik2bdrE5MmTmThxIgBz586lbdu2BHlhLKuiAgYMMIfOfvUry8uJiPg9\nn7qos6SkhC5durgfR0ZGUlRURFJSkntbcHAwM2fOrPHeYcOGnXP/DoeDsLAwwsLCcDgcOByOevd4\n5IgZNOnp5lRnEZFA4nQ6cTqduFyuGicEDeFTYWM1p9PZoPcbBgwfDp06mbd2FhEJND//Q9xTo0U+\nNRstOjqasrIy9+PS0lJiY2Nt7Ki6P/wBysthzhzNPBMRqQ+fCpvWrVsD5ow0l8tFfn4+MTExNndl\neu89mD4dli6Fli3t7kZExL/YFjbJycnExcVRXl5Ohw4dyM3NBSArK4v09HQSEhLIyMggJCTErhbd\nNm+GESPMwLniCru7ERHxP1qI8xy++w5iYmDSJBgyxILGRER8mBbi9IJjx8wZZw88oKAREWkIndmc\nwa5dMGqU+f3MokXQRLEsIo2QzmwscvQoTJwIPXpAbKx5MzQFjYhIwzSq62zO5YMPzNsEXH+9uRxN\nWJjdHYmIBAaFDbBjBzz6KGzfbk5vvuMOuzsSEQksjXqA6PBhePppc7isVy/YskVBIyJihUYZNoYB\nf/0rREbCN9+Y19FMmAAXXmh3ZyIiganRDaNt2wajR5vXz8ydC+exFqeIiNRTozqzeewxc7isXz/Y\nuFFBIyLiLY3qzGb/fti6FUJD7e5ERKRx0UWdIiJyRrqoU0RE/IbCRkRELKewERERyylsRETEcgob\nERGxnMJGREQs53dhs3PnToYPH87gwYPd215//XXGjRvHvHnzbOxMRETOxO/C5uqrr+aNN95wP964\ncSMrVqygadOmRERE2NhZdU6nU/X8tF4gfzbVUz272BY2aWlphIaGEhUVVW17YWEhERERhIeHk52d\nfc79rFmzht69ezNlyhSmT59uVbv1FugHYCDXC+TPpnqqZxfbwiY1NZW8vLwa2zMzM8nJyaGgoIBp\n06ZRWVnJ/PnzGTt2LBUVFTVef/311xMcHExQUBAnT570Rut14nK5VM9P6wXyZ1M91bOLbWETHx/P\npZdeWm1bVVUVAL169aJjx44kJiZSXFxMSkoKU6dO5YorrmD//v2MHDmSjRs3MnnyZG666SZ27NjB\nuHHjSEpKsuOj1CrQD8BArhfIn031VM8uPrUQZ0lJCV26dHE/joyMpKioqFqIBAcHM3PmzGrvmzhx\n4jn3fdlllxEUFOS5ZutA9fy3XiB/NtVTvfq47LLLPLIfnwobK1VWVtrdgohIo+VTs9Gio6MpKytz\nPy4tLSU2NtbGjkRExBN8Kmxat24NmDPSXC4X+fn5xMTE2NyViIg0lG1hk5ycTFxcHOXl5XTo0IHc\n3FwAsrKySE9PJyEhgYyMDEJCQuxqUUREPKTR3DxNRETs41PDaN5U27I3VikrK2PUqFE89NBDLF68\n2PJ6TqeT+Ph4Ro0axerVqy2vt3btWkaNGsXDDz9Mz549La/3z3/+k5EjR/LCCy9QUFBgWZ3ajhEr\nj5va9m3lsVNbPSuPndrqWXXs1FbLyuNm6dKljBgxgrS0NNavX3/GHqysZ+WxUlu9eh8rRiN39913\ne63WsWPHjCFDhlheZ/Xq1UafPn2MzMxM4x//+Ifl9f5jyZIlxqxZsyyvM2fOHCMvL88wDMMYNmyY\n5fVqO0asPG5q27eVx87p9bxx7NT2+aw6dk6v5Y3j5rvvvjNGjhx5xh68Uc/KY+X0evU9VgLqzMZT\nS+BYUWvZsmX07t2be+65x/J68fHxfPTRRzz66KO88sorltf7j7fffpv77rvP8nq/+c1vWLt2LY8/\n/jg7duywrI4neKJefY6dhtar77HjqX+fdTl2GlqrvsfN+dSbPHky6enp59y3VfWsPlZOr1fv3zOW\nxJ9NCgsLjc8//9zo2rVrte3du3c3Vq9ebbhcLqNz587Gvn373M+d718d51PLMAzjrrvu8lq9qqoq\n46GHHvJKvV27dhkPP/zwedU6n3qGYRjHjx83UlNTLa/TkDMbTx6TdTl2PFWvrseOJ+rV9djx1Ger\n63FT13qVlZXGqVOnjPHjxxsFBQU19uPpY+Vc9QzDs8fKuerV9VgJqIs64+PjayzlcPoSOIB7CZyb\nbrqJp556ik2bNjF58mQef/xxy2pdfPHFLF68GMMwznv8tj71jh8/zooVKzhx4gSjRo2yvF5SUhJz\n5swhLS3tvGrVt17Xrl156aWXCAoKIjMz07I6tR0j+/fvr9dxcz71/rMU0+OPP87q1avrdew0tN57\n771Xr2PHE//N1fXYaWitXbt21eu4qWu9oqIivvnmG1atWsXBgwfZsWMH6enplh0rZ6pn1bFypnr1\nPVYCKmxqc7YlcH6+7I1VtV588UVuueUWj9Y6V72BAwd6rV5SUhLPPfecV+vNmjXLK3V+fozUtlyS\nlfVuueWWBh879ak3cODABh879f1vriHHTn1qdezYscHHzdn+mxs9enS111p5rNRWz8pjpbZ69T1W\nAuo7GxER8U0BHzbeXALH28vtqJ5/1VE9/66leg0T8GHjzSVwvL3cjur5Vx3V8+9aqtdA55xC4EeG\nDBlitG/f3rjwwguNq666ypgzZ45hGIbhdDqNLl26GJ06dTJeffVVv6ulep6rF6ifqzHUC+TP1hjq\nabkaERGxXMAPo4mIiP0UNiIiYjmFjYiIWE5hIyIillPYiIiI5RQ2IiJiOYWNiIhYTmEjIiKWU9iI\niIjlFDYiNjt27Nh5PSfiTwL+fjYidquoqGDt2rVcfPHFfPnll0RERJCUlATABx98QGxsLM2bN6/1\nvdu3b6eqqoqePXt6s2URj9OZjYiFqqqqePrpp+nXrx99+/blscceIzc3l/3797Nnzx4OHDhASEjI\nGd/ftWtXli9fzvHjx73YtYjnKWxELDR06FDGjBlDixYt3Ns6derEZ599Rm5ubp3udHj//ffz/vvv\nW9mmiOUUNiIWqaioYPfu3XTv3r3G9pYtW7J3715atmwJwMmTJ1m0aBEvvfQSs2fPZuTIkXzzzTcA\nhIeHK2zE7ylsRCzyySef1LjLoWEYrF+/nh49enD06FH39s2bN9OvXz86duxIkyZNGDJkCO3btweg\nWbNmnDhxwqu9i3iawkbEIocOHaJVq1bVtq1cuZI77riDNm3aVAubHj160Lx5c4qLi3E4HDgcDvdZ\nD0DTpk291reIFRQ2IhZxOByUlJS4H+/bt4+ZM2fy4osvAnDgwAH3cyUlJVRWVrJ161auvvpq1q5d\nW21fBw8e9E7TIhbR1GcRi3Tu3Jn09HQmTJhAt27d2Lt3LwsWLHCf7Zw+aSAvL482bdrQs2dPlixZ\nwpVXXul+7siRIwQHB3u9fxFP0m2hRWwyZcoUUlNTadu27VlfV1JSws6dO7nnnnu81JmI52kYTcQm\n6enpZGVlnfN1b731FnfffbcXOhKxjs5sRGyUn5/PNddcQ6dOnWp9vqSkhKZNm9KjRw8vdybiWQob\nERGxnIbRRETEcgobERGxnMJGREQsp7ARERHLKWxERMRyChsREbGcwkZERCynsBEREcv9fye2FJiC\nzwdUAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 11 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the absolute error, that is the absolute difference from $x = 0.2$ compared to $\\tilde{x}$, increases with increasing magnitude of the numbers that we perform the arithemtic operation on. Above a certain magnitude, the error appears consant ... we'll come back to this.\n", "\n", "Second, the relative error:\n", "\n", "\\begin{equation}\n", "\\epsilon = |\\frac{\\tilde{x}-x}{x}|\n", "\\end{equation}\n", "\n", "Let's see if this provides more insight:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "relative_error = np.abs((x_tilde-x)/x)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 12 }, { "cell_type": "code", "collapsed": false, "input": [ "plt.plot(ls, relative_error)\n", "plt.xscale('log')\n", "plt.yscale('log')\n", "plt.xlabel('$O(x)$')\n", "plt.ylabel('$O(\\epsilon)$')\n", "plt.title('Relative Error')\n", "plt.ylim(None, 1.0E1)\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEkCAYAAAD5BJxYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9clfX9//EHmjNNhhLlJ5eCOQpIyjQmYejRITZJi1IX\nFhmUolSSK2vrs60+9XXLZYkRTU2lZmblstSaJGRH8gdIM8soQs2jTbYpaohOJPT6/nFtZyKoIOdw\nnR/P++12bu26zjnX63Xcm/M67/f1vt5XgGEYBiIiIm7UzuoERETE96nYiIiI26nYiIiI26nYiIiI\n26nYiIiI26nYiIiI26nYiJyFzWZj4cKF5/XePXv2EBgYiK4uEFGxET8QFhZG586dCQwM5Prrr+fx\nxx+ntra2We8NCAggICCg2XHWrl3r3O7Vqxc1NTXNfn9zvfLKK7Rv357AwEDn44c//CH/+Mc/XBpH\nxJVUbMTnBQQE8N5771FTU8Mrr7zC0qVL+fOf/+yWOG3Vixk0aBA1NTXOx+HDh/mf//mfRq+rr69v\n1r6zaenrRZqiYiN+pW/fvowYMYJVq1Y59+3cuZNHH32U0NBQJk6cyJdfftnke3fu3MmwYcMICQnh\nmmuuYebMmRw5cgSA1NRU9uzZw6hRowgMDGTWrFk4HA7atWvHyZMnefPNN4mJiWlwvNmzZ3PLLbcA\n5hf6W2+9xbBhw+jXrx8LFy6krq7ujJ/jbEUtLCyMl156ibi4OLp27crOnTtp164dy5Yto2/fvgwf\nPhyAlStXMnz4cKKjo5k7dy7/+te/AJx5n/56kVYxRHxcWFiYUVhYaBiGYXz66adGaGio8cILLxiG\nYRj19fXGpZdeauTl5RmHDx82Xn31VePyyy93vtdmsxkLFy40DMMwduzYYRQWFhp1dXXGZ599ZvTv\n3994+eWXG8T58MMPndu7du0yAgICjBMnThhHjx41AgMDje3btzufv/76640333zTMAzDmDNnjjFs\n2DDjiy++MHbs2GHYbDZj/vz5TX6evLw848Ybbzzr542KijKKioqM2tpaZx633nqrsXPnTqO2ttZY\nu3at0atXL6OgoMCoqKgwfvrTnxpPPPFEg7xPfb1Ia6nYiM8LDQ01unTpYnTp0sUICAgwJk2a5Hxu\nzZo1xvDhwxu8vl+/fsbmzZsNw2hYbE738ssvGzfffLNz+2zFxjAM46677jKeeuopwzAMo6KiwggM\nDDSOHTtmGIZhDBo0yNiwYYPzve+8844xcuTIJuPm5eUZF1xwgdG1a1fn48c//nGDPP4T59Q8ioqK\nnPumTp1q/OpXv3JuFxQUGNdcc80ZXy/SWhpGE58XEBDAihUrqK6u5t1332Xx4sVs3boVgMLCQj7+\n+GO6devmfOzYsYOioqJGxzly5AhZWVnExMQQFBTEtGnT+Pzzz5udx/jx41m6dCkAr7/+OsnJyVx4\n4YUcPXqUjRs3kpSU5MzhnnvuYePGjWc8VmxsLIcOHXI+tm/f3uD5gQMHNnrPqfs2btzIgAEDnNsD\nBgxg27Zt1NTUnPUYIudLxUb8Rrt27Rg9ejRTp07l4YcfBmDYsGHYbLYGX9w1NTXO50+Vm5vL119/\nzVtvvcV3333H7NmzOXnypPP59u3bn/VcSkJCAvv37+ezzz7jjTfeYPz48QBcdNFFDBw4kA8++MCZ\nw3fffcehQ4fO+7NecMEFZ903aNAgPvnkE+f2J598QnR0NIGBgWc9hsj5UrERv/PII49QXFxMSUkJ\nw4cPZ9u2bfzpT3/i0KFD1NbWYrfb2bt3b6P3VVZW0q1bNy699FJKS0t58cUXGzw/YMAA/vrXv54x\nbocOHRg7diyPPPIIhw4danDiPTU1ld/+9rds2bKFkydPsnfvXtasWeO6D32aW265haVLl7J27Vp2\n7NjBs88+S3JystviiajYiN8JCQlhwoQJzJw5k3bt2mG32/n6668ZMGAAvXr14rnnnmuyhzJt2jSO\nHTtGaGgoDz/8MJmZmQ2uoZk8eTLvvfcewcHBPP/88wCNrrEZP348H374IWPHjqVdu//++U2cOJH0\n9HR++9vfEhwczPDhw6moqGgy/4CAADZt2tTgOpvAwMCzFrrT87DZbMyePZvf/e533Hrrrdxyyy1M\nnz79jK8Xaa0A42z9fi+wa9cuZsyYQXV1NcuWLbM6HRERaYLX92x69+7NggULrE5DRETOwiOLTXp6\nOt27dyc6OrrB/qKiIiIjIwkPDycnJ8ei7EREpKU8stikpaWRn5/faH9WVhbz5s2jsLCQ3Nxcqqqq\nLMhORERayiOLTXx8PN26dWuwr7q6GoDBgwcTGhpKYmIiJSUlHDx4kMmTJ7N161ZmzpxpRboiInIO\nXjORvrS0lIiICOd2VFQUxcXFJCUlMXfu3HO+PyQkhAMHDrgzRRERn3PxxRe7ZBTJI3s27nDgwAEM\nc3meNnkMGTJE8bwwluIpnuI1fLjqR7rXFJuYmBjKy8ud22VlZcTGxlqY0dmFhYUpnhfGUjzFUzz3\n8JpiExQUBJgz0hwOBwUFBR69dpOvN0AVG8VTPP+I5yoeWWxSUlKIi4ujoqKCnj17kpeXB0B2djYZ\nGRkkJCSQmZlJSEiIxZmemc1mUzwvjKV4iqd47uH1Kwg0V1veRVFExFe46rvTI3s2IiLiW1RsRETE\n7VRsRETE7VRsRETE7VRsRETE7VRsRETE7VRsRETE7bxmIU4RscaRI/D001BYaHUm4s1UbESkSYYB\ny5bBww/DsGHw0ktwgb4x/M7117vmOGo6ItJIeTk88ADs2wevvw7x8VZnJN5O52xExOnoUfjlL83i\nMmoUbNmiQiOuoWIjIhgGvP02REbC3r3w+eeQlaVhM3EdNSURP1dRAQ8+aBaZxYthyBCrMxJfpJ6N\niJ86ehT+938hLg5GjIBPP1WhEfdRsRHxM4YB77wDUVGwa5c5ZPaLX0CHDlZnJr5Mw2gifmTvXpg4\nERwOeOUVGDrU6ozEX6hnI+InjhyBkSPhuutg61YVGmlbPnGnzl27djFjxgyqq6tZtmxZk6/RnTrF\nn508CbffDsHBsGABBARYnZF4C92p8xS9e/dmwYIFVqch4rF+/Ws4cAD++EcVGrGGRxWb9PR0unfv\nTnR0dIP9RUVFREZGEh4eTk5OjkXZiXinJUvgjTfM62h+8AOrsxF/5VHFJi0tjfz8/Eb7s7KymDdv\nHoWFheTm5lJVVcXixYuZNm0alZWVFmQq4h1KSmDaNFi5Ei65xOpsxJ95VLGJj4+nW7duDfZVV1cD\nMHjwYEJDQ0lMTKSkpITU1FRmz55Njx49OHjwIJMnT2br1q3MnDnTitRFPM6338Jtt8HChdC3r9XZ\niL/z+KnPpaWlREREOLejoqIoLi4mKSnJuS84OJi5c+ee81g2m42wsDDCwsKw2WzYbDZ3pCxiuaNH\n4ZZbzCVnRo2yOhvxJna7HbvdjsPhwOFwuOy4Hl9sXMlut1udgojbnTwJEyZAdDRMn251NuJtTv8h\nHuCiGSUeNYzWlJiYGMrLy53bZWVlxMbGWpiRiGf7v/+DykqYN08zz8RzeHyxCQoKAswZaQ6Hg4KC\nAgYOHGhxViKe6c034dVXzeVoLrzQ6mxE/sujik1KSgpxcXFUVFTQs2dP8vLyAMjOziYjI4OEhAQy\nMzMJCQmxOFMRz1Naat7wbMUK6N7d6mxEGvKJFQSaQysIiC/buxcGDoScHEhOtjob8SVaQUBEAPjX\nv+DWWyEzU4VGPJd6NiJezDAgJQXat4fXXtOEAHE9V313+tXUZxFf8//+n3lPGrtdhUY8m4qNiJd6\n+214+WVzSZpOnazORuTsVGxEvND778PkyfDBB3DZZVZnI3JuKjYiXmT3bnjoIfjiC3jrLejf3+qM\nRJpHs9FEvMDx4/C738GAAeZj2zbdaVO8i3o2Ih5uzRrzYs3ISPPCzd69rc5IpOVUbEQ81J498Itf\nwKefwpw5cPPNVmckcv40jCbiYerq4JlnzPMx0dFQVqZCI95PPRsRD1JYaA6Z/fjHsHkzXHGF1RmJ\nuIaKjYgH+Nvf4OGHzQLzwgu64Zn4Hg2jiViovh6efRb69YOrrjKHzFRoxBepZyNiEcOA+++Hr76C\n4mJz6EzEV6nYiFgkJwc2bjQfgYFWZyPiXio2IhZYswZ+/3vYtEmFRvyDTxSb8vJy5syZQ11dHUlJ\nSdx2221WpyRyRuXlcNdd5kKaYWFWZyPSNnzqfjZ1dXVMmDCBpUuXNnpO97MRT3DwIMTGwi9/Cenp\nVmcjcm4+eafO9PR0unfvTnR0dIP9RUVFREZGEh4eTk5OTpPvXblyJUOHDmXcuHFtkapIi33/PYwb\nZ842U6ERf+NRPZuPP/6YLl26cPfdd7Nt2zbn/uuuu445c+YQGhrKiBEjWL9+PatXr2bLli1Mnz6d\nHj16OF87evRoVq5c2ejY6tmI1e6/37zR2apV5p01RbyBT96pMz4+HofD0WBfdXU1AIMHDwYgMTGR\nkpISUlNTSU1NBWDdunUsX74cwzAYO3Zsm+Ys0hwvvQQffWROCFChEX/kUcWmKaWlpURERDi3o6Ki\nKC4uJikpyblvyJAhDBkyxIr0RM7pww/hqadgwwYICrI6GxFreHyxcSWbzUZYWBhhYWHYbDZsNpvV\nKYmP274dxo+HN9+EPn2szkbk3Ox2O3a7HYfD0WikqTU8vtjExMQwffp053ZZWRk33XTTeR3Lbre7\nKCuRc/vuO3MywNNPg37XiLc4/Yd4QECAS47rUbPRmhL073GHoqIiHA4HBQUFDBw40OKsRM6uvh5+\n/nNITIRJk6zORsR6HlVsUlJSiIuLo6Kigp49e5KXlwdAdnY2GRkZJCQkkJmZSUhIiMWZipzdI4+Y\na589/7zVmYh4Bo+a+uxOmvosbeXll2HWLHNxzW7drM5GpHVc9d2pYiPiQuvWmRdufvwxXHml1dmI\ntJ5PriAg4s2++cY8T7NkiQqNyOlUbERc4PBhc+bZb34DCQlWZyPieTSMJtJKJ07A6NHQq5e5UoCL\nZoqKeAQNo4l4iMceg9paeOEFFRqRM/H4izpFPFleHqxYASUl0KGD1dmIeC4No4mcp/Xr4bbbzBlo\nkZFWZyPiHhpGE7GQwwFjx8Kf/qRCI9IcKjYiLVRTY04IeOwxOM9l+kT8jobRRFrg5ElIToZLL4X5\n8zUhQHyfT948TcTTPf64uZrzsmUqNCItoWIj0kyLF8Nbb8HmzfCDH1idjYh30TCaSDNs2mSep7Hb\n4eqrrc5GpO1oNppIG9mzB26/3bymRoVG5Pyo2IicxdGjcMst8ItfwM03W52NiPfSMJrIGZw8CWPG\nQFAQLFqkCQHinzQbTcTNnngC9u2DpUtVaERayyeG0ex2O/Hx8UyZMoV169ZZnY74gKVLzdlny5dD\nx45WZyPi/XyiZ9OuXTu6dOlCx44dueKKK6xOR7zc5s0wdSp8+KF58aaItJ5H9WzS09Pp3r070dHR\nDfYXFRURGRlJeHg4OTk5jd4XHx/P6tWreeihh5g1a1ZbpSs+qLTUXFxzwQK45hqrsxHxHR5VbNLS\n0sjPz2+0Pysri3nz5lFYWEhubi5VVVUsXryYadOmUVlZScC/B9SDg4M5evRoW6ctPuDAAcjIMGee\nPfec+V8RcR2PGkaLj4/H4XA02FddXQ3A4MGDAUhMTKSkpITU1FRSU1MBeOedd/jggw+or69nypQp\nbZqzeLeTJ2HhQvj1r+GOO+DLL6FrV6uzEvE9HlVsmlJaWkpERIRzOyoqiuLiYpKSkpz7kpOTSU5O\nPuexbDYbYWFhhIWFYbPZsNls7khZvMQnn8D998MFF8AHH0C/flZnJGI9u92O3W7H4XA0+vHfGh5f\nbFzJbrdbnYJ4gIMH4X//F955B555Bu6+G9p51ICyiHVO/yEe4KJ5/x7/JxYTE0N5eblzu6ysjNjY\nWAszEm/1nyGzqCho3x6++gruuUeFRqQteHzPJigoCDBnpPXq1YuCggKeeOIJi7MSb7NlizlkBvCX\nv0D//tbmI+JvPOo3XUpKCnFxcVRUVNCzZ0/y8vIAyM7OJiMjg4SEBDIzMwkJCbE4U/EWhw7BAw/A\nyJEwcSJs2KBCI2IFrY0mPskw4E9/gl/+Em69FWbMgOBgq7MS8T5aG03kLJ57zlw8c9UquP56q7MR\nEfVsxOe8/z5MmgTFxdCzp9XZiHg39WxEmlBWBmlpsHKlCo2IJ/GoCQIirVFVBaNGwfPPg2bHi3gW\nDaOJT6irg+HDIS4Ofv97q7MR8R2u+u5UsRGvZxjmOZp9+8xVAXSRpojr6JyNyL+98AKUlJjX0KjQ\niHgmFRvxavn55vpmxcUQGGh1NiJyJio24rW++spcRHP5cggNtTobETkbDTqIVzpwwJx59oc/wI03\nWp2NiJyLJgiI1/n+e7jpJrjuOtBdwEXcS7PRWkjFxndkZsLu3eaFm+3bW52NiG/TbDTxS7m5sG4d\nbNqkQiPiTdSzEa9RUACpqbBxI1xxhdXZiPgH9WzEr1RUwJ13wltvqdCIeKNmF5u6ujpWrVrFqlWr\nqK+vp3379tTU1BAcHExiYiJjxoyhna6oEzc4dMiceTZjBpxya3QR8SLNGkbbsGED77//PnfddRfh\n4eF06NDB+dyxY8f44osveO2115gwYQL9PfQ2iBpG80719eZdNqOiIDvb6mxE/E+bzUY7fvw427dv\np2/fvuc8WGlpKTExMa1OqqXWr1/PkiVLqK+v58svv2TDhg2NXqNi452mToWvvzbvUXOBBn1F2pxl\nU5//8pe/MHLkSHbv3k1oaCgVFRVceeWVrU7EFVasWMG+ffuYOHFio+dUbLzPvHkwe7a5FE3XrlZn\nI+KfXPXd2eKTLJs3bwbAbrcD8Le//a3VSfxHeno63bt3Jzo6usH+oqIiIiMjCQ8PJycn54zvf/31\n1xk/frzL8hHrfPQR/Pa35m2dVWhEvF+Li821117LsGHDWLNmDStXrmTbtm0uSyYtLY38/PxG+7Oy\nspg3bx6FhYXk5uZSVVXF4sWLmTZtGpWVlQDs2bOHoKAgLrroIpflI9bYsQNSUmDpUggPtzobEXGF\nFo+CJycnc/XVV5Ofn893333H5MmTXZZMfHw8Doejwb7q6moABg8eDEBiYiIlJSWkpqaSmprqfN2i\nRYtIT093WS5ijepqGD0anngChg2zOhsRcZVzFpvjx49TU1NDSEiIc9+VV17Z5HmanTt30qdPH5cm\nWFpaSkREhHM7KiqK4uJikpKSGrzuySefPOexbDYbYWFhhIWFYbPZsGkerUc5ccLs0QwdClOmWJ2N\niH+y2+3Y7XYcDkejH/+tcc5i07FjRwoKCjh8+DDJycl06tSp0Wv2799PdnY2NpvN5cXGlf5znkk8\n0/Tp5u2dNcVZxDqn/xAPCAhwyXGbNYx288038/e//53Zs2ezb98+amtrqa2t5fDhw1x44YX069eP\nRx99lKCgIJckdaqYmBimT5/u3C4rK+Omm25yeRyx1sKF8N575h03T7mMS0R8RLPP2XTs2JFrr72W\nTp060b9/f7q20RSh/xSwoqIievXqRUFBAU888USbxJa2UVQEv/qV+d9u3azORkTcoVmz0d5//30m\nTZpEUVEReXl5DBw4kKlTpzpP3rtKSkoKcXFxVFRU0LNnT/Ly8gDIzs4mIyODhIQEMjMzG5w/Eu+2\naxf8/Ofw2mtwyqk5EfExzbqoc9asWTzyyCPObcMwWLNmDXPnzmXJkiV07tzZrUm6gi7q9Dw1NRAX\nB5MmwYMPWp2NiDSlTS/qjIqKahR8xIgRzJ8/n7lz57Y6CfE/J07A+PFmsXngAauzERF3a1axOXLk\nCLW1tY32X3LJJVx66aUuT0p83+OPw5Ej8OKL4KLJLiLiwZpVbG688UYyMzPZunVro+c+//xzlycl\nvu3VV+Htt+HPf9bMMxF/0eyFOKurq7nzzjvZu3cvsbGxdO7cmdLSUjIyMrjzzjvdnWer6ZyNZ9i4\nEW69Fex287YBIuLZLFv1ubKyko8++ogjR45gs9m46qqrWp1EW1Cxsd7u3XDDDeY1NT/7mdXZiEhz\nWFZsvJWKjbW2b4fbb4e0NJg2zepsRKS5LLvFgEhL/Otf8JvfmD2atDR46CGrMxIRK+jeh+IWhgEr\nV5rFZeBA+Owz+NGPrM5KRKyiYiMut3OneTvnb76BBQvgpz+1OiMRsZqG0cRljh2DJ580ezJDhpi9\nGRUaEQH1bMRF3nvP7M1cfz18+in07Gl1RiLiSVRspFV27YKsLPj6a5g3D4YPtzojEfFEGkaT81Jb\nC089BTEx5kyzzz9XoRGRM1PPRlosP99cPPPaa2HLFujVy+qMRMTTqdhIi7zzjlloFi4E3TBVRJpL\nKwhIs23dag6VrV5tTgQQEd+nFQSkTf3zn3DLLZCbq0IjIi3nE8No3377LTNmzKBHjx7ExcWRkJBg\ndUo+pbYWkpPhnntg3DirsxERb+QTw2h5eXn06NGDESNGcM899/DKK680eo2G0c6PYcCECeYFm2++\nCe3UFxbxKz45jJaenk737t2Jjo5usL+oqIjIyEjCw8PJyclp9L7bb7+d9evX89hjj7Fjx462Stcv\nPPsslJWZNzxToRGR8+VRPZuPP/6YLl26cPfdd7Nt2zbn/uuuu445c+YQGhrKiBEjWL9+PatXr2bL\nli1Mnz6dHj16APD999+TkZHBokWLGh1bPZuWW7kSpkyBkhK4/HKrsxERK7jqu9OjztnEx8fjcDga\n7KuurgZg8ODBACQmJlJSUkJqaiqpqakA7N69mxkzZhAQEEBWVlab5uyrtm2De+81l6FRoRGR1vKo\nYtOU0tJSIiIinNtRUVEUFxeTlJTk3BcaGsr8+fPPeSybzUZYWBhhYWHYbDZsNps7UvZ6+/fD6NEw\nZ465qKaI+A+73Y7dbsfhcDT68d8aHl9sXMlut1udgserq4PbboPx482HiPiX03+IBwQEuOS4Hn/K\nNyYmhvLycud2WVkZsbGxFmbkuwzDPEcTEgJPP211NiLiSzy+2AQFBQHmjDSHw0FBQQEDNbbjFtnZ\n8Ne/wuLFmnkmIq7lUV8pKSkpxMXFUVFRQc+ePcnLywMgOzubjIwMEhISyMzMJCQkxOJMfc/q1eY0\n5xUroEsXq7MREV/jUVOf3UlTn8/syy/BZoN334W4OKuzERFP4pMXdUrbO3DAnHk2a5YKjYi4j3o2\nfqyuDkaMgJ/8BGbOtDobEfFErvruVLHxU/+ZebZ3rzl81r691RmJiCfyyRUEpO28+CKsXw+bNqnQ\niIj7qdj4oTVr4He/g40bITDQ6mxExB+o2PiZr7+G1FT485+hd2+rsxERf6HZaH7k4EEYNcrs1cTH\nW52NiPgTTRDwE99/Dz/7GVxzDTz/vNXZiIi30Gy0FvL3YvPAA7Bzp3nLAE0IEJHm0mw0abY//hHW\nrtXMMxGxjno2Pm7tWvNWARs2QJ8+VmcjIt5GPRs5p+3bISUF3nhDhUZErKXZaD7qu+/MNc+eegqG\nDrU6GxHxdxpG80H19XDzzXDllfDCC1ZnIyLeTKs+yxlNnw4nT2qKs4h4Dp2z8TELFsBf/gLFxXCB\n/t8VEQ+hYTQfsm4djBsHH39sDqGJiLSW3w6j7dq1i/vuu4+xY8c2ue2vvvkGfv5zWLJEhUZEPI/X\nFZvevXuzYMGCM277o8OHzTXPfvMbSEiwOhsRkcYsKzbp6el0796d6OjoBvuLioqIjIwkPDycnJwc\ni7LzHv/4B9x+OwweDJmZVmcjItI0y4pNWloa+fn5jfZnZWUxb948CgsLyc3NpaqqisWLFzNt2jQq\nKystyNQz1deb05qjo6F/f/N/BwRYnZWISNMsKzbx8fF069atwb7q6moABg8eTGhoKImJiZSUlJCa\nmsrs2bPp0aMHBw8eZPLkyWzdupWZM2c22vYHGzbAgAGwYgUUFcHMmdChg9VZiYicmUdNji0tLSUi\nIsK5HRUVRXFxMUlJSc59wcHBzJ07t8H7Tt8+E5vNRlhYGGFhYdhsNmw2m0vybiv79sGjj0JhITz3\nnDnzTL0ZEXElu92O3W7H4XDgcDhcdlyPKjbuZrfbrU7hvJw4AXPnwpNPwj33wFdf6XbOIuIep/8Q\nD3DRL1qPmo0WExNDeXm5c7usrIzY2FgLM7Lepk0QE2Pextluh2efVaEREe/jUcUmKCgIMGekORwO\nCgoKGDhwoMVZWWP/frj3Xhgzxlx+Zu1auPpqq7MSETk/lhWblJQU4uLiqKiooGfPnuTl5QGQnZ1N\nRkYGCQkJZGZmEhISYlWKljhxwrzZ2dVXQ1CQOWSWkqJzMyLi3bRcjQcpLYUpU6BzZ8jNNac1i4hY\nSTdP8zGffAIjR5orNd91l3oyIuJb1LPxAJWVMHCgeWFmcrLV2YiI/JffLsTpa44dg1tvhcmTVWhE\nxHepZ2Mhw4Dx46FdO3jtNQ2diYjn0TkbHzBjhnlrALtdhUZEfJuKjUXefhvmz4eSEujUyepsRETc\nS8NoFvj0U0hMhPx8c0FNERFPpQkCXuof/zAnBLz0kgqNiPgPFZs2VFtrFpr0dPDzu1iLiJ/RMFob\nMQy4+26oq4M33tCEABHxDpqN5mVmzjTXOSsqUqEREf+jYtMGVqyAF180Z5517mx1NiIibU/Fxs0+\n/xzuuw/efx9+9COrsxERsYYmCLjRvn0werS55tlPfmJ1NiIi1lGxcZPjx+G22yA11bwfjYiIP9Ns\nNDcwDHN6c3W1eTvndirpIuKlNBvNgz33HGzdCuvXq9CIiIAXDqPt2rWL++67j7H/viqyvLycKVOm\ncO+997J8+XKLszMnAjz/vDkD7aKLrM5GRMQzeO0w2tixY1m2bJlzu66ujgkTJrB06dImX98Ww2hl\nZTB0qFlobrjBraFERNqE16+Nlp6eTvfu3YmOjm6wv6ioiMjISMLDw8nJyWnWsVauXMnQoUMZN26c\nO1Jtlqoqc+bZc8+p0IiInM6yYpOWlkZ+fn6j/VlZWcybN4/CwkJyc3Opqqpi8eLFTJs2jcrKyiaP\nNXr0aDaaFw0fAAAKv0lEQVRs2EBeXp67025SXR2MGWOud5aaakkKIiIezbIJAvHx8Tgcjgb7qqur\nARg8eDAAiYmJlJSUkJqaSuq/v8UPHjzI448/ztatW3nmmWe44YYbWL58OYZhOM/jtCXDgAcegB/+\n0LwZmoiINOZRs9FKS0uJiIhwbkdFRVFcXExSUpJzX3BwMHPnzm3wviFDhjTr+DabjbCwMMLCwrDZ\nbNhstlbnnJMDmzbBxo3Qvn2rDyciYim73Y7dbsfhcDTqELSGRxUbd7Pb7S493po18Pvfm4UmMNCl\nhxYRscTpP8QDXLRysEdNfY6JiaG8vNy5XVZWRmxsrIUZndnXX8Ndd8Gbb0Lv3lZnIyLi2Tyq2AQF\nBQHmjDSHw0FBQQEDBw60OKvGDh6EUaPMXs2/Ty+JiMhZWFZsUlJSiIuLo6Kigp49ezpnkmVnZ5OR\nkUFCQgKZmZmEhIRYlWKTvv8exo2DpCS4916rsxER8Q5ee1FnS7nqwqQHH4QdO2DVKrjAr854iYg/\n0tpoFpg7FwoLobhYhUZEpCXUs2mmjz6CO+4wF9cMD3dhYiIiHszrl6vxJjt2mIVm6VIVGhGR86Fi\ncw7V1eaaZ08+CcOGWZ2NiIh30jDaWZw4YU5x7t0bcnPdlJiIiAfTMJqbHTkCkyaZi2xmZ1udjYiI\nd1OxOY1hmKsCREaahWbZMujQweqsRES8mybwnuKLL8zraA4ehNdfh/h4qzMSEfEN6tlgTgKYNs2c\nADBmDPz1ryo0IiKu5NfF5uRJePVViIgwz9GUlcH99+uCTRERV/Pbr9UtW8ybntXXw4oV8JOfWJ2R\niIjv8ruezYEDMGUKjBxpLqRZXKxCIyLibn5VbObNg6goc3bZV1+ZxaadX/0LiIhYw6+G0ZYsMe+u\nee21VmciIuJf/GoFgZMnDVx0h1MREb+gFQTOgwqNiIg1/KrYiIiINbyu2OzatYv77ruPsWPHAmC3\n24mPj2fKlCmsW7fO4uxERKQpXldsevfuzYIFC5zb7dq1o0uXLnTs2JErrrjCwswastvtiueFsRRP\n8RTPPSwrNunp6XTv3p3o6OgG+4uKioiMjCQ8PJycnJxzHic+Pp7Vq1fz0EMPMWvWLHel22K+3gBV\nbBRP8fwjnqtYVmzS0tLIz89vtD8rK4t58+ZRWFhIbm4uVVVVLF68mGnTplFZWdno9QH/PusfHBzM\n0aNH3Z53czkcDsXzwliKp3iK5x6WXWcTHx/f6B+turoagMGDBwOQmJhISUkJqamppKamAnDw4EEe\nf/xxtm7dyjPPPMNVV13FBx98QH19PVOmTGnTz3A2vt4AVWwUT/H8I56reNRFnaWlpURERDi3o6Ki\nKC4uJikpybkvODiYuXPnNnhfcnLyOY998cUXO3tBbUXxvDOW4ime4v3XxRdf7JLjeFSxcaeqqiqr\nUxAR8VseNRstJiaG8vJy53ZZWRmxsbEWZiQiIq7gUcUmKCgIMGekORwOCgoKGDhwoMVZiYhIa1lW\nbFJSUoiLi6OiooKePXuSl5cHQHZ2NhkZGSQkJJCZmUlISIhVKYqIiIv4zUKcIiJiHY8aRmtLpy97\n407l5eVMmTKFe++9l+XLl7s9Xlsv4bN+/XqmTJnCxIkTGTRokNvjffvtt0yePJmnnnqKwsJCt8Vp\nqo24s900dWx3tp2m4rmr7TQVy53tpql47mw3K1asYNKkSaSnp7N58+Yz5uDOeO5sK03Fa3FbMfzc\nmDFj2izW8ePHjTvuuMPtcdatW2fcdNNNRlZWlrFnzx63x/uPd99915g/f77b4yxatMjIz883DMMw\nJkyY4PZ4TbURd7abpo7tzrZzajx3t52mPps7282p8dqi3fzzn/80Jk+efMYc2iKeO9vKqfFa2lZ8\nqmfjqiVw3BFr5cqVDB06lHHjxrk9niuW8Dmff8vXX3+d8ePHuz3e7bffzvr163nsscfYsWOH2+K4\ngivitaTttDZeS9qOq/4tm9tuWhuvpe3mfOLNnDmTjIyMcx7bXfHc3VZOjdfi7xm3lD+LFBUVGVu2\nbDH69u3bYH+/fv2MdevWGQ6Hw7jqqquM/fv3O587318d5xPLMAxj1KhRbRavurrauPfee9sk3u7d\nu42JEyeeV6zziWcYhlFXV2ekpaW5PU5rejaubJPNaTuuitectuOKWC1pN676bM1tN82NV1VVZZw8\nedKYPn26UVhY2Og4rm4r54pnGK5tK+eK19zvGZ+6qLMlS+DccMMNzmVvZs6cyWOPPea2WF26dGH5\n8uUYhnHe47ctiVdXV9fqJXxaEi8pKYlFixaRnp5+XrFaGq9v377MmDGDgIAAsrKy3BanqTZy6nJJ\nzWk35xPv008/dR573bp1LWo7rY33zjvvNLvtuOLvrSXtprXxdu/e3aJ209x4xcXFfPPNN6xdu5aa\nmhp27NhBRkaG29rKmeK5q62cKV5L2gr4wQoCZ1sC5/Rlb9wV6+mnn2bIkCEujXWueM1ZwsdV8ZKS\nknjyySfbNN78+fPbJM7pbaSp5ZLcGW/IkCGtbjstiZecnNyqttPSv7fWtpuWxAsNDW11uznb39yD\nDz7Y4LXubCtNxXNnW2kqXkvbik+dsxEREc/k88WmLZfAaevldhTPu+L4Qzxf/myK1zo+X2zacgmc\ntl5uR/G8K44/xPPlz6Z4rXTOKQRe5I477jAuu+wy4wc/+IFx+eWXG4sWLTIMwzDsdrsRERFh9OnT\nx5gzZ47XxVI818Xz1c9lRTxf/myK5/p4Wq5GRETczueH0URExHoqNiIi4nYqNiIi4nYqNiIi4nYq\nNiIi4nYqNiIi4nYqNiIi4nYqNiIi4nYqNiIi4nYqNiIWO378+Hk9J+JNfP5+NiJWq6ysZP369XTp\n0oUvv/ySyMhIkpKSAHjvvfeIjY2lY8eOTb53+/btVFdXM2jQoLZMWcTl1LMRcaPq6mp+/etfM3r0\naEaOHMkjjzxCXl4eBw8e5O9//zuHDx8mJCTkjO/v27cv77//PnV1dW2YtYjrqdiIuNGdd97J1KlT\nufDCC537+vTpw6ZNm8jLy2vWnQ7vuusuVq1a5c40RdxOxUbETSorK9m7dy/9+vVrtL9Tp07s27eP\nTp06AXDixAmWLVvGjBkzWLhwIZMnT+abb74BIDw8XMVGvJ6KjYibfPTRR43ucmgYBps3b6Z///7U\n1tY693/22WeMHj2a0NBQ2rVrxx133MFll10GQIcOHaivr2/T3EVcTcVGxE2OHDlC586dG+xbs2YN\nI0aMoGvXrg2KTf/+/enYsSMlJSXYbDZsNpuz1wPQvn37NstbxB1UbETcxGazUVpa6tzev38/c+fO\n5emnnwbg8OHDzudKS0upqqriiy++oHfv3qxfv77BsWpqatomaRE30dRnETe56qqryMjI4NFHH+Wa\na65h3759LFmyxNnbOXXSQH5+Pl27dmXQoEG8++67/OhHP3I+d+zYMYKDg9s8fxFX0m2hRSzyhz/8\ngbS0NC655JKzvq60tJRdu3Yxbty4NspMxPU0jCZikYyMDLKzs8/5utdee40xY8a0QUYi7qOejYiF\nCgoKuOKKK+jTp0+Tz5eWltK+fXv69+/fxpmJuJaKjYiIuJ2G0URExO1UbERExO1UbERExO1UbERE\nxO1UbERExO1UbERExO1UbERExO1UbERExO3+P+W0wvBqPCzmAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 13 }, { "cell_type": "markdown", "metadata": {}, "source": [ "To sum up our observations, we can say that the error in our little subtraction increases with the magnitude of the number leading the decimal point, up to a point where it is constant, i.e., of the exact size as the subtracted number. In other words, our subtraction of $0.2$ has no effect for numbers larger than a certain magnitude. More on that later...\n", "\n", "A pressing question by now should be: How can we do reliable computations if any operation on floating introduces errors? Particularly since we know that so many safety critical aspects are based on computer calculations.\n", "The answer can be summarized by:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "1.2-1. == 1.2-1." ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 14, "text": [ "True" ] } ], "prompt_number": 14 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This, simply put, means that even though errors are part of computations, these errors are not random (like in experimental measurement apparati). They are systematic, reproducible and well understood. The answer is not that the errors are *small*, since this would be very much a matter of context.\n", "\n", "The two errors that significantly contribute to deviations from arithmetic in ${\\rm I\\!R}$ are due to mechanisms of:\n", "\n", "* rounding\n", "* truncation\n", "\n", "In brief, errors result from a *finite* representation of numbers on binary systems and rules to operate on them, which brings us to..." ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Binary number representation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In 1986, Link faced a problem: no matter his achievements, he could never accumulate more than 255 rupees in his wallet. This was not a consequence of reigning communism - he faced the same problem as Psy in the year 2014: the finite representation of numbers on binary systems.\n" ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Integer numbers" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we learn how computer systems represent, i.e., store, numbers, let's look at a systematic way of representing numeric values, here the natural number $123$, using our base 10 numeral system:\n", "\n", "$1\\times10^2 + 2\\times10^1 + 3\\times10^0 = 123$\n", "\n", "This can be ordered in a column manner, say for an arbitrary four columns:\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
$10^3$$10^2$$10^1$$10^0$
$0$$1$$2$$3$
\n", "\n", "so that, again, the expression yields:\n", "\n", "$0\\times10^3 +1\\times10^2 + 2\\times10^1 + 3\\times10^0 = 123$\n", "\n", "Since our factor ranges from 0 to 9 (10 values), we can represent all natural numbers $[0,9999]$ this way. If we allocated more cells, say five:\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
$10^4$$10^3$$10^2$$10^1$$10^0$
$0$$0$$1$$2$$3$
\n", "\n", "it would be $[0,99999]$, etc...\n", "\n", "Conventional computer systems, however, store data in terms of *low voltage* or *high voltage*, or *0* and *1*, *false* and *true*, *off* and *on*. So the obvious base for a binary number system is two:\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
$2^7$$2^6$$2^5$$2^4$$2^3$$2^2$$2^1$$2^0$
$0$$1$$1$$1$$1$$0$$1$$1$
\n", "\n", "Following the same rationale, we get:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "0*2**7 + 1*2**6 + 1*2**5 + 1*2**4 + 1*2**3 + 0*2**2 + 1*2**1 + 1*2**0" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 15, "text": [ "123" ] } ], "prompt_number": 15 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This concept is used by the operating system to allocate a fixed width of memory, here 8 bits, and interpret its values as multipliers for decreasing powers of two, bitwise, from the left to the right, to compute a number from $0 / I$ values. It immediately follows that the range of numbers is a function of the numbers of bits allocated to represent the integer value. Common bit chunk sizes and their limits are found here.\n", "\n", "This system has been extended to negative numbers, allowing for simple addition/subtraction operations - the two's complement, standard for integer representation and operations on modern computer systems. The system is simple, really:\n", "\n", "$123$ is, as derived above, given as:\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
$0$$1$$1$$1$$1$$0$$1$$1$
\n", "\n", "but now, the left most bit has a different meaning: it is the *significant* bit, or *sign* bit. A $0$ value indicates the following representation is for a positive number. To get $-123$ we first invert all bits\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
$1$$0$$0$$0$$0$$1$$0$$0$
\n", "\n", "and then add $1$ to get\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
$1$$0$$0$$0$$0$$1$$0$$1$
\n", "\n", "the two's complement representation of $-123$, with the signigicant bit of $1$ indicating a negative number. The operation of addition is very simply peformed in this system. For $123+123$, we first need to allocate more memory, to avoid overflow, since the result is not within $[-128,127]$. So for 16 bits:\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
$0$$0$$0$$0$$0$$0$$0$$0$$0$$1$$1$$1$$1$$0$$1$$1$
\n", "\n", "The addition of **int_1 + int_2 = resul**, where **int_1** = $123$ and **int_2** = $123$, can then be represented as\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
carry$0$$0$$0$$0$$0$$0$$0$$0$$0$$1$$1$$1$$1$$0$$1$$1$
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
int_1$0$$0$$0$$0$$0$$0$$0$$0$$0$$1$$1$$1$$1$$0$$1$$1$
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
int_2$0$$0$$0$$0$$0$$0$$0$$0$$0$$1$$1$$1$$1$$0$$1$$1$
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
resul$0$$0$$0$$0$$0$$0$$0$$0$$1$$1$$1$$1$$0$$1$$1$$0$
\n", "\n", "where **resul** = 246, check for yourself (we start at the first non-zero bit from the left):\n" ] }, { "cell_type": "code", "collapsed": false, "input": [ "1*2**7 + 1*2**6 + 1*2**5 + 1*2**4 + 0*2**3 + 1*2**2 + 1*2**1 + 0*2**0" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 16, "text": [ "246" ] } ], "prompt_number": 16 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Subtraction is simply implemented by adding the negative value, so $123-123$ becomes $123+(-123)$:\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
carry$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$$1$
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
int_1$0$$0$$0$$0$$0$$0$$0$$0$$0$$1$$1$$1$$1$$0$$1$$1$
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
int_2$1$$1$$1$$1$$1$$1$$1$$1$$1$$0$$0$$0$$0$$1$$0$$1$
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
resul$0$$0$$0$$0$$0$$0$$0$$0$$0$$0$$0$$0$$0$$0$$0$$0$
\n", "\n", "which yields $0$, as expected. This is all we'll look at here in terms of arithmetic operations on binary systems. The take-home message from this section is that integers are represented on computers in a two's complement manner (base 2 numeral system), and are limited in the *range* of representable numbers by the dedicated memory." ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Floating point number representation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Like integers, floating point numbers are limited in the *range* of numbers that can be represented, a function of the allocated memory. But they also face another limitation: one of *precision*.\n", "\n", "To illustrate the concept of a *floating* point, think about two applications:\n", "\n", "* a conventional hydrocarbon resevoir model on the km scale:\n", "\n", "\n", "\n", "* a pore model on the micron scale:\n", "\n", "\n", "\n", "Numerical values describing the model geometry and processes therein vary drastically in magnitude. Most likely, however, we're not concerned with mm scale errors on the field scale, and we're not expecting any values larger than mm on the pore scale. This motivates a numeral system that represents numbers with a *relative* precision, i.e., a *floating* point - as opposed to a *fixed* point. You are familiar with one system of this kind, namely the scientific notation, a decimal (base 10) floating point notation. For example\n", "\n", "$123000$ = 1.23**E3** = 1.23$\\times10^3$\n", "\n", "$0.00123$ = 1.23**E-3** = 1.23$\\times10^{-3}$\n", "\n", "A floating point representation has two components, a *fraction* or *significand* (here $1.23$) and an *exponent* (here $3$ and $-3$), around a base, here $10$ in the decimal system. It follwos\n", "\n", "$value = fraction \\times base^{exponent}$\n", "\n", "The standard bit representation of floating point numbers , here for single precison (32 bit), adheres to this concept:\n", "\n", "\n", "\n", "The base, again, is $2$. The numeric value of the bit pattern follows as\n", "\n", "\\begin{equation}\n", "value = (-1)^{sign} \\left( 1 + \\sum_{i=1}^{23} b_{23-i}2^{-i} \\right) \\times 2^{(exponent-127)}\n", "\\end{equation}\n", "\n", "similar to \n", "\n", "$value = sign \\times fraction \\times base^{exponent}$\n", "\n", "where the $fraction$ is also referred to as $mantissa$. The zero offset of $127$ in the exponent, also referred to as exponent bias, arises mostly to allow for the representation of special values and zero, but we won't go into details here. The special numbers are\n", "\n", "$NaN$ ... Not a Number\n", "\n", "and\n", "\n", "$-/+ \\infty$\n", "\n", "These values invoke non-standard arithmetic operations and significantly slow down the process compared to conventional binary numbers. Think about why it would make sense to initialize variables to NaN and not to zero! Look up the order of magnitude that results in $\\infty$.\n", "\n", "The floating point format allows for a **wide range** of numbers to be represented with **relative precision**. It is obvious that both are a function of allocated memory, i.e., 16bit (half-), 32bit (single-) or 64bit (double-precision). Specifically, the finiteness of the fraction limits the precision, and the finiteness of the exponent limits the range.\n", "\n", "Let's think about this, recalling our initial assessment of the relative error in the calculation\n", "\n", "$1.2 - 1., 10.2-10., 100.2-100...$\n" ] }, { "cell_type": "code", "collapsed": false, "input": [ "plt.plot(ls, relative_error)\n", "plt.xscale('log')\n", "plt.yscale('log')\n", "plt.xlabel('$O(x)$')\n", "plt.ylabel('$O(\\epsilon)$')\n", "plt.title('Relative Error')\n", "plt.ylim(None, 1.0E1)\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEkCAYAAAD5BJxYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9clfX9//EHmjNNhhLlJ5eCOQpIyjQmYejRITZJi1IX\nFhmUolSSK2vrs60+9XXLZYkRTU2lZmblstSaJGRH8gdIM8soQs2jTbYpaohOJPT6/nFtZyKoIOdw\nnR/P++12bu26zjnX63Xcm/M67/f1vt5XgGEYBiIiIm7UzuoERETE96nYiIiI26nYiIiI26nYiIiI\n26nYiIiI26nYiIiI26nYiJyFzWZj4cKF5/XePXv2EBgYiK4uEFGxET8QFhZG586dCQwM5Prrr+fx\nxx+ntra2We8NCAggICCg2XHWrl3r3O7Vqxc1NTXNfn9zvfLKK7Rv357AwEDn44c//CH/+Mc/XBpH\nxJVUbMTnBQQE8N5771FTU8Mrr7zC0qVL+fOf/+yWOG3Vixk0aBA1NTXOx+HDh/mf//mfRq+rr69v\n1r6zaenrRZqiYiN+pW/fvowYMYJVq1Y59+3cuZNHH32U0NBQJk6cyJdfftnke3fu3MmwYcMICQnh\nmmuuYebMmRw5cgSA1NRU9uzZw6hRowgMDGTWrFk4HA7atWvHyZMnefPNN4mJiWlwvNmzZ3PLLbcA\n5hf6W2+9xbBhw+jXrx8LFy6krq7ujJ/jbEUtLCyMl156ibi4OLp27crOnTtp164dy5Yto2/fvgwf\nPhyAlStXMnz4cKKjo5k7dy7/+te/AJx5n/56kVYxRHxcWFiYUVhYaBiGYXz66adGaGio8cILLxiG\nYRj19fXGpZdeauTl5RmHDx82Xn31VePyyy93vtdmsxkLFy40DMMwduzYYRQWFhp1dXXGZ599ZvTv\n3994+eWXG8T58MMPndu7du0yAgICjBMnThhHjx41AgMDje3btzufv/76640333zTMAzDmDNnjjFs\n2DDjiy++MHbs2GHYbDZj/vz5TX6evLw848Ybbzzr542KijKKioqM2tpaZx633nqrsXPnTqO2ttZY\nu3at0atXL6OgoMCoqKgwfvrTnxpPPPFEg7xPfb1Ia6nYiM8LDQ01unTpYnTp0sUICAgwJk2a5Hxu\nzZo1xvDhwxu8vl+/fsbmzZsNw2hYbE738ssvGzfffLNz+2zFxjAM46677jKeeuopwzAMo6KiwggM\nDDSOHTtmGIZhDBo0yNiwYYPzve+8844xcuTIJuPm5eUZF1xwgdG1a1fn48c//nGDPP4T59Q8ioqK\nnPumTp1q/OpXv3JuFxQUGNdcc80ZXy/SWhpGE58XEBDAihUrqK6u5t1332Xx4sVs3boVgMLCQj7+\n+GO6devmfOzYsYOioqJGxzly5AhZWVnExMQQFBTEtGnT+Pzzz5udx/jx41m6dCkAr7/+OsnJyVx4\n4YUcPXqUjRs3kpSU5MzhnnvuYePGjWc8VmxsLIcOHXI+tm/f3uD5gQMHNnrPqfs2btzIgAEDnNsD\nBgxg27Zt1NTUnPUYIudLxUb8Rrt27Rg9ejRTp07l4YcfBmDYsGHYbLYGX9w1NTXO50+Vm5vL119/\nzVtvvcV3333H7NmzOXnypPP59u3bn/VcSkJCAvv37+ezzz7jjTfeYPz48QBcdNFFDBw4kA8++MCZ\nw3fffcehQ4fO+7NecMEFZ903aNAgPvnkE+f2J598QnR0NIGBgWc9hsj5UrERv/PII49QXFxMSUkJ\nw4cPZ9u2bfzpT3/i0KFD1NbWYrfb2bt3b6P3VVZW0q1bNy699FJKS0t58cUXGzw/YMAA/vrXv54x\nbocOHRg7diyPPPIIhw4danDiPTU1ld/+9rds2bKFkydPsnfvXtasWeO6D32aW265haVLl7J27Vp2\n7NjBs88+S3JystviiajYiN8JCQlhwoQJzJw5k3bt2mG32/n6668ZMGAAvXr14rnnnmuyhzJt2jSO\nHTtGaGgoDz/8MJmZmQ2uoZk8eTLvvfcewcHBPP/88wCNrrEZP348H374IWPHjqVdu//++U2cOJH0\n9HR++9vfEhwczPDhw6moqGgy/4CAADZt2tTgOpvAwMCzFrrT87DZbMyePZvf/e533Hrrrdxyyy1M\nnz79jK8Xaa0A42z9fi+wa9cuZsyYQXV1NcuWLbM6HRERaYLX92x69+7NggULrE5DRETOwiOLTXp6\nOt27dyc6OrrB/qKiIiIjIwkPDycnJ8ei7EREpKU8stikpaWRn5/faH9WVhbz5s2jsLCQ3Nxcqqqq\nLMhORERayiOLTXx8PN26dWuwr7q6GoDBgwcTGhpKYmIiJSUlHDx4kMmTJ7N161ZmzpxpRboiInIO\nXjORvrS0lIiICOd2VFQUxcXFJCUlMXfu3HO+PyQkhAMHDrgzRRERn3PxxRe7ZBTJI3s27nDgwAEM\nc3meNnkMGTJE8bwwluIpnuI1fLjqR7rXFJuYmBjKy8ud22VlZcTGxlqY0dmFhYUpnhfGUjzFUzz3\n8JpiExQUBJgz0hwOBwUFBR69dpOvN0AVG8VTPP+I5yoeWWxSUlKIi4ujoqKCnj17kpeXB0B2djYZ\nGRkkJCSQmZlJSEiIxZmemc1mUzwvjKV4iqd47uH1Kwg0V1veRVFExFe46rvTI3s2IiLiW1RsRETE\n7VRsRETE7VRsRETE7VRsRETE7VRsRETE7VRsRETE7bxmIU4RscaRI/D001BYaHUm4s1UbESkSYYB\ny5bBww/DsGHw0ktwgb4x/M7117vmOGo6ItJIeTk88ADs2wevvw7x8VZnJN5O52xExOnoUfjlL83i\nMmoUbNmiQiOuoWIjIhgGvP02REbC3r3w+eeQlaVhM3EdNSURP1dRAQ8+aBaZxYthyBCrMxJfpJ6N\niJ86ehT+938hLg5GjIBPP1WhEfdRsRHxM4YB77wDUVGwa5c5ZPaLX0CHDlZnJr5Mw2gifmTvXpg4\nERwOeOUVGDrU6ozEX6hnI+InjhyBkSPhuutg61YVGmlbPnGnzl27djFjxgyqq6tZtmxZk6/RnTrF\nn508CbffDsHBsGABBARYnZF4C92p8xS9e/dmwYIFVqch4rF+/Ws4cAD++EcVGrGGRxWb9PR0unfv\nTnR0dIP9RUVFREZGEh4eTk5OjkXZiXinJUvgjTfM62h+8AOrsxF/5VHFJi0tjfz8/Eb7s7KymDdv\nHoWFheTm5lJVVcXixYuZNm0alZWVFmQq4h1KSmDaNFi5Ei65xOpsxJ95VLGJj4+nW7duDfZVV1cD\nMHjwYEJDQ0lMTKSkpITU1FRmz55Njx49OHjwIJMnT2br1q3MnDnTitRFPM6338Jtt8HChdC3r9XZ\niL/z+KnPpaWlREREOLejoqIoLi4mKSnJuS84OJi5c+ee81g2m42wsDDCwsKw2WzYbDZ3pCxiuaNH\n4ZZbzCVnRo2yOhvxJna7HbvdjsPhwOFwuOy4Hl9sXMlut1udgojbnTwJEyZAdDRMn251NuJtTv8h\nHuCiGSUeNYzWlJiYGMrLy53bZWVlxMbGWpiRiGf7v/+DykqYN08zz8RzeHyxCQoKAswZaQ6Hg4KC\nAgYOHGhxViKe6c034dVXzeVoLrzQ6mxE/sujik1KSgpxcXFUVFTQs2dP8vLyAMjOziYjI4OEhAQy\nMzMJCQmxOFMRz1Naat7wbMUK6N7d6mxEGvKJFQSaQysIiC/buxcGDoScHEhOtjob8SVaQUBEAPjX\nv+DWWyEzU4VGPJd6NiJezDAgJQXat4fXXtOEAHE9V313+tXUZxFf8//+n3lPGrtdhUY8m4qNiJd6\n+214+WVzSZpOnazORuTsVGxEvND778PkyfDBB3DZZVZnI3JuKjYiXmT3bnjoIfjiC3jrLejf3+qM\nRJpHs9FEvMDx4/C738GAAeZj2zbdaVO8i3o2Ih5uzRrzYs3ISPPCzd69rc5IpOVUbEQ81J498Itf\nwKefwpw5cPPNVmckcv40jCbiYerq4JlnzPMx0dFQVqZCI95PPRsRD1JYaA6Z/fjHsHkzXHGF1RmJ\nuIaKjYgH+Nvf4OGHzQLzwgu64Zn4Hg2jiViovh6efRb69YOrrjKHzFRoxBepZyNiEcOA+++Hr76C\n4mJz6EzEV6nYiFgkJwc2bjQfgYFWZyPiXio2IhZYswZ+/3vYtEmFRvyDTxSb8vJy5syZQ11dHUlJ\nSdx2221WpyRyRuXlcNdd5kKaYWFWZyPSNnzqfjZ1dXVMmDCBpUuXNnpO97MRT3DwIMTGwi9/Cenp\nVmcjcm4+eafO9PR0unfvTnR0dIP9RUVFREZGEh4eTk5OTpPvXblyJUOHDmXcuHFtkapIi33/PYwb\nZ842U6ERf+NRPZuPP/6YLl26cPfdd7Nt2zbn/uuuu445c+YQGhrKiBEjWL9+PatXr2bLli1Mnz6d\nHj16OF87evRoVq5c2ejY6tmI1e6/37zR2apV5p01RbyBT96pMz4+HofD0WBfdXU1AIMHDwYgMTGR\nkpISUlNTSU1NBWDdunUsX74cwzAYO3Zsm+Ys0hwvvQQffWROCFChEX/kUcWmKaWlpURERDi3o6Ki\nKC4uJikpyblvyJAhDBkyxIr0RM7pww/hqadgwwYICrI6GxFreHyxcSWbzUZYWBhhYWHYbDZsNpvV\nKYmP274dxo+HN9+EPn2szkbk3Ox2O3a7HYfD0WikqTU8vtjExMQwffp053ZZWRk33XTTeR3Lbre7\nKCuRc/vuO3MywNNPg37XiLc4/Yd4QECAS47rUbPRmhL073GHoqIiHA4HBQUFDBw40OKsRM6uvh5+\n/nNITIRJk6zORsR6HlVsUlJSiIuLo6Kigp49e5KXlwdAdnY2GRkZJCQkkJmZSUhIiMWZipzdI4+Y\na589/7zVmYh4Bo+a+uxOmvosbeXll2HWLHNxzW7drM5GpHVc9d2pYiPiQuvWmRdufvwxXHml1dmI\ntJ5PriAg4s2++cY8T7NkiQqNyOlUbERc4PBhc+bZb34DCQlWZyPieTSMJtJKJ07A6NHQq5e5UoCL\nZoqKeAQNo4l4iMceg9paeOEFFRqRM/H4izpFPFleHqxYASUl0KGD1dmIeC4No4mcp/Xr4bbbzBlo\nkZFWZyPiHhpGE7GQwwFjx8Kf/qRCI9IcKjYiLVRTY04IeOwxOM9l+kT8jobRRFrg5ElIToZLL4X5\n8zUhQHyfT948TcTTPf64uZrzsmUqNCItoWIj0kyLF8Nbb8HmzfCDH1idjYh30TCaSDNs2mSep7Hb\n4eqrrc5GpO1oNppIG9mzB26/3bymRoVG5Pyo2IicxdGjcMst8ItfwM03W52NiPfSMJrIGZw8CWPG\nQFAQLFqkCQHinzQbTcTNnngC9u2DpUtVaERayyeG0ex2O/Hx8UyZMoV169ZZnY74gKVLzdlny5dD\nx45WZyPi/XyiZ9OuXTu6dOlCx44dueKKK6xOR7zc5s0wdSp8+KF58aaItJ5H9WzS09Pp3r070dHR\nDfYXFRURGRlJeHg4OTk5jd4XHx/P6tWreeihh5g1a1ZbpSs+qLTUXFxzwQK45hqrsxHxHR5VbNLS\n0sjPz2+0Pysri3nz5lFYWEhubi5VVVUsXryYadOmUVlZScC/B9SDg4M5evRoW6ctPuDAAcjIMGee\nPfec+V8RcR2PGkaLj4/H4XA02FddXQ3A4MGDAUhMTKSkpITU1FRSU1MBeOedd/jggw+or69nypQp\nbZqzeLeTJ2HhQvj1r+GOO+DLL6FrV6uzEvE9HlVsmlJaWkpERIRzOyoqiuLiYpKSkpz7kpOTSU5O\nPuexbDYbYWFhhIWFYbPZsNls7khZvMQnn8D998MFF8AHH0C/flZnJGI9u92O3W7H4XA0+vHfGh5f\nbFzJbrdbnYJ4gIMH4X//F955B555Bu6+G9p51ICyiHVO/yEe4KJ5/x7/JxYTE0N5eblzu6ysjNjY\nWAszEm/1nyGzqCho3x6++gruuUeFRqQteHzPJigoCDBnpPXq1YuCggKeeOIJi7MSb7NlizlkBvCX\nv0D//tbmI+JvPOo3XUpKCnFxcVRUVNCzZ0/y8vIAyM7OJiMjg4SEBDIzMwkJCbE4U/EWhw7BAw/A\nyJEwcSJs2KBCI2IFrY0mPskw4E9/gl/+Em69FWbMgOBgq7MS8T5aG03kLJ57zlw8c9UquP56q7MR\nEfVsxOe8/z5MmgTFxdCzp9XZiHg39WxEmlBWBmlpsHKlCo2IJ/GoCQIirVFVBaNGwfPPg2bHi3gW\nDaOJT6irg+HDIS4Ofv97q7MR8R2u+u5UsRGvZxjmOZp9+8xVAXSRpojr6JyNyL+98AKUlJjX0KjQ\niHgmFRvxavn55vpmxcUQGGh1NiJyJio24rW++spcRHP5cggNtTobETkbDTqIVzpwwJx59oc/wI03\nWp2NiJyLJgiI1/n+e7jpJrjuOtBdwEXcS7PRWkjFxndkZsLu3eaFm+3bW52NiG/TbDTxS7m5sG4d\nbNqkQiPiTdSzEa9RUACpqbBxI1xxhdXZiPgH9WzEr1RUwJ13wltvqdCIeKNmF5u6ujpWrVrFqlWr\nqK+vp3379tTU1BAcHExiYiJjxoyhna6oEzc4dMiceTZjBpxya3QR8SLNGkbbsGED77//PnfddRfh\n4eF06NDB+dyxY8f44osveO2115gwYQL9PfQ2iBpG80719eZdNqOiIDvb6mxE/E+bzUY7fvw427dv\np2/fvuc8WGlpKTExMa1OqqXWr1/PkiVLqK+v58svv2TDhg2NXqNi452mToWvvzbvUXOBBn1F2pxl\nU5//8pe/MHLkSHbv3k1oaCgVFRVceeWVrU7EFVasWMG+ffuYOHFio+dUbLzPvHkwe7a5FE3XrlZn\nI+KfXPXd2eKTLJs3bwbAbrcD8Le//a3VSfxHeno63bt3Jzo6usH+oqIiIiMjCQ8PJycn54zvf/31\n1xk/frzL8hHrfPQR/Pa35m2dVWhEvF+Li821117LsGHDWLNmDStXrmTbtm0uSyYtLY38/PxG+7Oy\nspg3bx6FhYXk5uZSVVXF4sWLmTZtGpWVlQDs2bOHoKAgLrroIpflI9bYsQNSUmDpUggPtzobEXGF\nFo+CJycnc/XVV5Ofn893333H5MmTXZZMfHw8Doejwb7q6moABg8eDEBiYiIlJSWkpqaSmprqfN2i\nRYtIT093WS5ijepqGD0anngChg2zOhsRcZVzFpvjx49TU1NDSEiIc9+VV17Z5HmanTt30qdPH5cm\nWFpaSkREhHM7KiqK4uJikpKSGrzuySefPOexbDYbYWFhhIWFYbPZsGkerUc5ccLs0QwdClOmWJ2N\niH+y2+3Y7XYcDkejH/+tcc5i07FjRwoKCjh8+DDJycl06tSp0Wv2799PdnY2NpvN5cXGlf5znkk8\n0/Tp5u2dNcVZxDqn/xAPCAhwyXGbNYx288038/e//53Zs2ezb98+amtrqa2t5fDhw1x44YX069eP\nRx99lKCgIJckdaqYmBimT5/u3C4rK+Omm25yeRyx1sKF8N575h03T7mMS0R8RLPP2XTs2JFrr72W\nTp060b9/f7q20RSh/xSwoqIievXqRUFBAU888USbxJa2UVQEv/qV+d9u3azORkTcoVmz0d5//30m\nTZpEUVEReXl5DBw4kKlTpzpP3rtKSkoKcXFxVFRU0LNnT/Ly8gDIzs4mIyODhIQEMjMzG5w/Eu+2\naxf8/Ofw2mtwyqk5EfExzbqoc9asWTzyyCPObcMwWLNmDXPnzmXJkiV07tzZrUm6gi7q9Dw1NRAX\nB5MmwYMPWp2NiDSlTS/qjIqKahR8xIgRzJ8/n7lz57Y6CfE/J07A+PFmsXngAauzERF3a1axOXLk\nCLW1tY32X3LJJVx66aUuT0p83+OPw5Ej8OKL4KLJLiLiwZpVbG688UYyMzPZunVro+c+//xzlycl\nvu3VV+Htt+HPf9bMMxF/0eyFOKurq7nzzjvZu3cvsbGxdO7cmdLSUjIyMrjzzjvdnWer6ZyNZ9i4\nEW69Fex287YBIuLZLFv1ubKyko8++ogjR45gs9m46qqrWp1EW1Cxsd7u3XDDDeY1NT/7mdXZiEhz\nWFZsvJWKjbW2b4fbb4e0NJg2zepsRKS5LLvFgEhL/Otf8JvfmD2atDR46CGrMxIRK+jeh+IWhgEr\nV5rFZeBA+Owz+NGPrM5KRKyiYiMut3OneTvnb76BBQvgpz+1OiMRsZqG0cRljh2DJ580ezJDhpi9\nGRUaEQH1bMRF3nvP7M1cfz18+in07Gl1RiLiSVRspFV27YKsLPj6a5g3D4YPtzojEfFEGkaT81Jb\nC089BTEx5kyzzz9XoRGRM1PPRlosP99cPPPaa2HLFujVy+qMRMTTqdhIi7zzjlloFi4E3TBVRJpL\nKwhIs23dag6VrV5tTgQQEd+nFQSkTf3zn3DLLZCbq0IjIi3nE8No3377LTNmzKBHjx7ExcWRkJBg\ndUo+pbYWkpPhnntg3DirsxERb+QTw2h5eXn06NGDESNGcM899/DKK680eo2G0c6PYcCECeYFm2++\nCe3UFxbxKz45jJaenk737t2Jjo5usL+oqIjIyEjCw8PJyclp9L7bb7+d9evX89hjj7Fjx462Stcv\nPPsslJWZNzxToRGR8+VRPZuPP/6YLl26cPfdd7Nt2zbn/uuuu445c+YQGhrKiBEjWL9+PatXr2bL\nli1Mnz6dHj16APD999+TkZHBokWLGh1bPZuWW7kSpkyBkhK4/HKrsxERK7jqu9OjztnEx8fjcDga\n7KuurgZg8ODBACQmJlJSUkJqaiqpqakA7N69mxkzZhAQEEBWVlab5uyrtm2De+81l6FRoRGR1vKo\nYtOU0tJSIiIinNtRUVEUFxeTlJTk3BcaGsr8+fPPeSybzUZYWBhhYWHYbDZsNps7UvZ6+/fD6NEw\nZ465qKaI+A+73Y7dbsfhcDT68d8aHl9sXMlut1udgserq4PbboPx482HiPiX03+IBwQEuOS4Hn/K\nNyYmhvLycud2WVkZsbGxFmbkuwzDPEcTEgJPP211NiLiSzy+2AQFBQHmjDSHw0FBQQEDNbbjFtnZ\n8Ne/wuLFmnkmIq7lUV8pKSkpxMXFUVFRQc+ePcnLywMgOzubjIwMEhISyMzMJCQkxOJMfc/q1eY0\n5xUroEsXq7MREV/jUVOf3UlTn8/syy/BZoN334W4OKuzERFP4pMXdUrbO3DAnHk2a5YKjYi4j3o2\nfqyuDkaMgJ/8BGbOtDobEfFErvruVLHxU/+ZebZ3rzl81r691RmJiCfyyRUEpO28+CKsXw+bNqnQ\niIj7qdj4oTVr4He/g40bITDQ6mxExB+o2PiZr7+G1FT485+hd2+rsxERf6HZaH7k4EEYNcrs1cTH\nW52NiPgTTRDwE99/Dz/7GVxzDTz/vNXZiIi30Gy0FvL3YvPAA7Bzp3nLAE0IEJHm0mw0abY//hHW\nrtXMMxGxjno2Pm7tWvNWARs2QJ8+VmcjIt5GPRs5p+3bISUF3nhDhUZErKXZaD7qu+/MNc+eegqG\nDrU6GxHxdxpG80H19XDzzXDllfDCC1ZnIyLeTKs+yxlNnw4nT2qKs4h4Dp2z8TELFsBf/gLFxXCB\n/t8VEQ+hYTQfsm4djBsHH39sDqGJiLSW3w6j7dq1i/vuu4+xY8c2ue2vvvkGfv5zWLJEhUZEPI/X\nFZvevXuzYMGCM277o8OHzTXPfvMbSEiwOhsRkcYsKzbp6el0796d6OjoBvuLioqIjIwkPDycnJwc\ni7LzHv/4B9x+OwweDJmZVmcjItI0y4pNWloa+fn5jfZnZWUxb948CgsLyc3NpaqqisWLFzNt2jQq\nKystyNQz1deb05qjo6F/f/N/BwRYnZWISNMsKzbx8fF069atwb7q6moABg8eTGhoKImJiZSUlJCa\nmsrs2bPp0aMHBw8eZPLkyWzdupWZM2c22vYHGzbAgAGwYgUUFcHMmdChg9VZiYicmUdNji0tLSUi\nIsK5HRUVRXFxMUlJSc59wcHBzJ07t8H7Tt8+E5vNRlhYGGFhYdhsNmw2m0vybiv79sGjj0JhITz3\nnDnzTL0ZEXElu92O3W7H4XDgcDhcdlyPKjbuZrfbrU7hvJw4AXPnwpNPwj33wFdf6XbOIuIep/8Q\nD3DRL1qPmo0WExNDeXm5c7usrIzY2FgLM7Lepk0QE2Pextluh2efVaEREe/jUcUmKCgIMGekORwO\nCgoKGDhwoMVZWWP/frj3Xhgzxlx+Zu1auPpqq7MSETk/lhWblJQU4uLiqKiooGfPnuTl5QGQnZ1N\nRkYGCQkJZGZmEhISYlWKljhxwrzZ2dVXQ1CQOWSWkqJzMyLi3bRcjQcpLYUpU6BzZ8jNNac1i4hY\nSTdP8zGffAIjR5orNd91l3oyIuJb1LPxAJWVMHCgeWFmcrLV2YiI/JffLsTpa44dg1tvhcmTVWhE\nxHepZ2Mhw4Dx46FdO3jtNQ2diYjn0TkbHzBjhnlrALtdhUZEfJuKjUXefhvmz4eSEujUyepsRETc\nS8NoFvj0U0hMhPx8c0FNERFPpQkCXuof/zAnBLz0kgqNiPgPFZs2VFtrFpr0dPDzu1iLiJ/RMFob\nMQy4+26oq4M33tCEABHxDpqN5mVmzjTXOSsqUqEREf+jYtMGVqyAF180Z5517mx1NiIibU/Fxs0+\n/xzuuw/efx9+9COrsxERsYYmCLjRvn0werS55tlPfmJ1NiIi1lGxcZPjx+G22yA11bwfjYiIP9Ns\nNDcwDHN6c3W1eTvndirpIuKlNBvNgz33HGzdCuvXq9CIiIAXDqPt2rWL++67j7H/viqyvLycKVOm\ncO+997J8+XKLszMnAjz/vDkD7aKLrM5GRMQzeO0w2tixY1m2bJlzu66ujgkTJrB06dImX98Ww2hl\nZTB0qFlobrjBraFERNqE16+Nlp6eTvfu3YmOjm6wv6ioiMjISMLDw8nJyWnWsVauXMnQoUMZN26c\nO1Jtlqoqc+bZc8+p0IiInM6yYpOWlkZ+fn6j/VlZWcybN4/CwkJyc3Opqqpi8eLFTJs2jcrKyiaP\nNXr0aDaaFw0fAAAKv0lEQVRs2EBeXp67025SXR2MGWOud5aaakkKIiIezbIJAvHx8Tgcjgb7qqur\nARg8eDAAiYmJlJSUkJqaSuq/v8UPHjzI448/ztatW3nmmWe44YYbWL58OYZhOM/jtCXDgAcegB/+\n0LwZmoiINOZRs9FKS0uJiIhwbkdFRVFcXExSUpJzX3BwMHPnzm3wviFDhjTr+DabjbCwMMLCwrDZ\nbNhstlbnnJMDmzbBxo3Qvn2rDyciYim73Y7dbsfhcDTqELSGRxUbd7Pb7S493po18Pvfm4UmMNCl\nhxYRscTpP8QDXLRysEdNfY6JiaG8vNy5XVZWRmxsrIUZndnXX8Ndd8Gbb0Lv3lZnIyLi2Tyq2AQF\nBQHmjDSHw0FBQQEDBw60OKvGDh6EUaPMXs2/Ty+JiMhZWFZsUlJSiIuLo6Kigp49ezpnkmVnZ5OR\nkUFCQgKZmZmEhIRYlWKTvv8exo2DpCS4916rsxER8Q5ee1FnS7nqwqQHH4QdO2DVKrjAr854iYg/\n0tpoFpg7FwoLobhYhUZEpCXUs2mmjz6CO+4wF9cMD3dhYiIiHszrl6vxJjt2mIVm6VIVGhGR86Fi\ncw7V1eaaZ08+CcOGWZ2NiIh30jDaWZw4YU5x7t0bcnPdlJiIiAfTMJqbHTkCkyaZi2xmZ1udjYiI\nd1OxOY1hmKsCREaahWbZMujQweqsRES8mybwnuKLL8zraA4ehNdfh/h4qzMSEfEN6tlgTgKYNs2c\nADBmDPz1ryo0IiKu5NfF5uRJePVViIgwz9GUlcH99+uCTRERV/Pbr9UtW8ybntXXw4oV8JOfWJ2R\niIjv8ruezYEDMGUKjBxpLqRZXKxCIyLibn5VbObNg6goc3bZV1+ZxaadX/0LiIhYw6+G0ZYsMe+u\nee21VmciIuJf/GoFgZMnDVx0h1MREb+gFQTOgwqNiIg1/KrYiIiINbyu2OzatYv77ruPsWPHAmC3\n24mPj2fKlCmsW7fO4uxERKQpXldsevfuzYIFC5zb7dq1o0uXLnTs2JErrrjCwswastvtiueFsRRP\n8RTPPSwrNunp6XTv3p3o6OgG+4uKioiMjCQ8PJycnJxzHic+Pp7Vq1fz0EMPMWvWLHel22K+3gBV\nbBRP8fwjnqtYVmzS0tLIz89vtD8rK4t58+ZRWFhIbm4uVVVVLF68mGnTplFZWdno9QH/PusfHBzM\n0aNH3Z53czkcDsXzwliKp3iK5x6WXWcTHx/f6B+turoagMGDBwOQmJhISUkJqamppKamAnDw4EEe\nf/xxtm7dyjPPPMNVV13FBx98QH19PVOmTGnTz3A2vt4AVWwUT/H8I56reNRFnaWlpURERDi3o6Ki\nKC4uJikpybkvODiYuXPnNnhfcnLyOY998cUXO3tBbUXxvDOW4ime4v3XxRdf7JLjeFSxcaeqqiqr\nUxAR8VseNRstJiaG8vJy53ZZWRmxsbEWZiQiIq7gUcUmKCgIMGekORwOCgoKGDhwoMVZiYhIa1lW\nbFJSUoiLi6OiooKePXuSl5cHQHZ2NhkZGSQkJJCZmUlISIhVKYqIiIv4zUKcIiJiHY8aRmtLpy97\n407l5eVMmTKFe++9l+XLl7s9Xlsv4bN+/XqmTJnCxIkTGTRokNvjffvtt0yePJmnnnqKwsJCt8Vp\nqo24s900dWx3tp2m4rmr7TQVy53tpql47mw3K1asYNKkSaSnp7N58+Yz5uDOeO5sK03Fa3FbMfzc\nmDFj2izW8ePHjTvuuMPtcdatW2fcdNNNRlZWlrFnzx63x/uPd99915g/f77b4yxatMjIz883DMMw\nJkyY4PZ4TbURd7abpo7tzrZzajx3t52mPps7282p8dqi3fzzn/80Jk+efMYc2iKeO9vKqfFa2lZ8\nqmfjqiVw3BFr5cqVDB06lHHjxrk9niuW8Dmff8vXX3+d8ePHuz3e7bffzvr163nsscfYsWOH2+K4\ngivitaTttDZeS9qOq/4tm9tuWhuvpe3mfOLNnDmTjIyMcx7bXfHc3VZOjdfi7xm3lD+LFBUVGVu2\nbDH69u3bYH+/fv2MdevWGQ6Hw7jqqquM/fv3O587318d5xPLMAxj1KhRbRavurrauPfee9sk3u7d\nu42JEyeeV6zziWcYhlFXV2ekpaW5PU5rejaubJPNaTuuitectuOKWC1pN676bM1tN82NV1VVZZw8\nedKYPn26UVhY2Og4rm4r54pnGK5tK+eK19zvGZ+6qLMlS+DccMMNzmVvZs6cyWOPPea2WF26dGH5\n8uUYhnHe47ctiVdXV9fqJXxaEi8pKYlFixaRnp5+XrFaGq9v377MmDGDgIAAsrKy3BanqTZy6nJJ\nzWk35xPv008/dR573bp1LWo7rY33zjvvNLvtuOLvrSXtprXxdu/e3aJ209x4xcXFfPPNN6xdu5aa\nmhp27NhBRkaG29rKmeK5q62cKV5L2gr4wQoCZ1sC5/Rlb9wV6+mnn2bIkCEujXWueM1ZwsdV8ZKS\nknjyySfbNN78+fPbJM7pbaSp5ZLcGW/IkCGtbjstiZecnNyqttPSv7fWtpuWxAsNDW11uznb39yD\nDz7Y4LXubCtNxXNnW2kqXkvbik+dsxEREc/k88WmLZfAaevldhTPu+L4Qzxf/myK1zo+X2zacgmc\ntl5uR/G8K44/xPPlz6Z4rXTOKQRe5I477jAuu+wy4wc/+IFx+eWXG4sWLTIMwzDsdrsRERFh9OnT\nx5gzZ47XxVI818Xz1c9lRTxf/myK5/p4Wq5GRETczueH0URExHoqNiIi4nYqNiIi4nYqNiIi4nYq\nNiIi4nYqNiIi4nYqNiIi4nYqNiIi4nYqNiIi4nYqNiIWO378+Hk9J+JNfP5+NiJWq6ysZP369XTp\n0oUvv/ySyMhIkpKSAHjvvfeIjY2lY8eOTb53+/btVFdXM2jQoLZMWcTl1LMRcaPq6mp+/etfM3r0\naEaOHMkjjzxCXl4eBw8e5O9//zuHDx8mJCTkjO/v27cv77//PnV1dW2YtYjrqdiIuNGdd97J1KlT\nufDCC537+vTpw6ZNm8jLy2vWnQ7vuusuVq1a5c40RdxOxUbETSorK9m7dy/9+vVrtL9Tp07s27eP\nTp06AXDixAmWLVvGjBkzWLhwIZMnT+abb74BIDw8XMVGvJ6KjYibfPTRR43ucmgYBps3b6Z///7U\n1tY693/22WeMHj2a0NBQ2rVrxx133MFll10GQIcOHaivr2/T3EVcTcVGxE2OHDlC586dG+xbs2YN\nI0aMoGvXrg2KTf/+/enYsSMlJSXYbDZsNpuz1wPQvn37NstbxB1UbETcxGazUVpa6tzev38/c+fO\n5emnnwbg8OHDzudKS0upqqriiy++oHfv3qxfv77BsWpqatomaRE30dRnETe56qqryMjI4NFHH+Wa\na65h3759LFmyxNnbOXXSQH5+Pl27dmXQoEG8++67/OhHP3I+d+zYMYKDg9s8fxFX0m2hRSzyhz/8\ngbS0NC655JKzvq60tJRdu3Yxbty4NspMxPU0jCZikYyMDLKzs8/5utdee40xY8a0QUYi7qOejYiF\nCgoKuOKKK+jTp0+Tz5eWltK+fXv69+/fxpmJuJaKjYiIuJ2G0URExO1UbERExO1UbERExO1UbERE\nxO1UbERExO1UbERExO1UbERExO1UbERExO3+P+W0wvBqPCzmAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 17 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We recognize that our error is expressed relative to the result $0.2$. We see that for a number of $10^{16}$ and greater the relative error is $1$, meaning the subtraction has no effect, i.e.,\n", "\n", "$ 10^{17} - 0.2 = 10^{17}$\n", "\n", "Let's express the error relative to the magnitude of the involved numbers." ] }, { "cell_type": "code", "collapsed": false, "input": [ "relative_error_2 = np.abs((x_tilde-x)/ls)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 18 }, { "cell_type": "code", "collapsed": false, "input": [ "plt.plot(ls, relative_error_2)\n", "plt.xscale('log')\n", "plt.yscale('log')\n", "plt.xlabel('$O(x)$')\n", "plt.ylabel('$O(\\epsilon)$')\n", "plt.title('Absolute Error')\n", "plt.ylim(None, 1.0E1)\n", "plt.annotate('$\\epsilon_{machine}$', xy=(1.1E1, 1.1E-16), xycoords='data', xytext=(1E6, 1E-10), textcoords='data', \n", " arrowprops=dict(arrowstyle=\"->\", connectionstyle=\"arc3\"), fontsize=20 )\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEkCAYAAAD5BJxYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XdYVNfWBvB3RIOoiUgEY0UlFrChXhTFghXUGEs0yrUg\n3ihqAqgJRnNjj7liil2QGyXW3E+NLRJR0AgSIxYgIsYgKilYESkiiMD+/jgPIwjiDMyZ+v6ehyeZ\ndtYaPJw1e885eymEEAJEREQyqqbrBIiIyPix2BARkexYbIiISHYsNkREJDsWGyIikh2LDRERyY7F\nhkzWlClTsHDhQo1u89tvv0Xv3r01uk0iY8BiQ0bP1dUVVlZWyM/PL3W/QqGAQqHQUVZVL0zffvst\nzMzM8Oqrryp/XnvtNdy5c0eDWRJpBosNGbWUlBScO3cONjY2OHz4sK7T0TgXFxdkZ2crf7KysvDG\nG2+UeV5BQYFK91VE3ecTlcRiQ0Zt+/btGDhwICZNmoRt27aVeTw7OxujR4/GG2+8gY8//hgPHjxQ\nPrZy5Uo4Ojqibt266NixIxITEwEAjx8/RmBgIDp06IDBgwfjhx9+KDd2SkoKqlWrhqKiIuV9rq6u\n2LJlC65evYoZM2bgl19+wauvvgorKysA0gF9z5496N+/PxwdHbFly5YyI7KSKloApHnz5ti0aRN6\n9uwJS0tLXL9+HdWqVcPevXvRvn17DBo0CABw+PBhDBo0CB06dEBQUBAeP35cKv/nn09UKYLIiNnZ\n2YmdO3eKpKQkUaNGDXH37l3lY56enqJ27dpi27ZtIjU1VYwfP16MHz9eCCHE5cuXRZs2bURqaqoQ\nQoirV6+K27dvCyGEWLRokejXr5/4/fffxYkTJ0Tz5s3FTz/9JIQQIiQkRPTq1UsIIcTNmzeFQqEQ\nhYWFypiurq5iy5YtQgghvv32W+Vzi61du1b0799fXL58WSQnJwtXV1cRHBxc7nsrGas8zZs3Fw4O\nDiIqKkrk5eUp8xk5cqS4fv26yMvLEydPnhTNmjUT4eHhIikpSQwYMEAsXry4VP4ln09UWRzZkNGK\njo5Gamoq3n77bbRq1QoODg7YvXt3qed06dIFkydPRqNGjbB06VIcO3YMRUVFKCwsRF5eHq5du4ai\noiK0adNGOT116NAhzJ8/H61bt0b//v0xYcIEHDhwQO38RDmjkj179mD58uVo164d7Ozs4Ofnh4MH\nD75wG2fPnkW9evWUP61atSr1+Pjx49G7d2+Ym5sr75s7dy5atmwJc3NzHDx4EBMmTMDAgQPRqlUr\nzJ8/v8x7Kfl8ospisSGjtW3bNgwePBivvvoqAGDs2LGlptIUCgU6deqkvN26dWs8ffoUv/32Gzp2\n7IgVK1Zg/vz5aNy4MRYtWoTHjx8jOzsbly5dQteuXZWv69q1K06fPl3lfHNycnDmzBkMGzZMWTym\nTJmCM2fOvPA1zs7OePjwofLn2rVrpR7v3r17mdeUvO/MmTNl3ktCQgKys7Mr3AaRulhsyCjl5uZi\nz549OHnyJBo2bIiGDRviq6++wq+//opLly4BkEYW8fHxytf8/vvvqFGjBuzt7QEAEyZMwC+//IKz\nZ8/i+PHjCAkJwauvvoqOHTviwoULytdduHABffr0KZODtbU1atSooTw7rKCgAAkJCcrHzczMSo1u\nateuje7du+PYsWPK4pGRkYGHDx9W+vdQvXr1Cu9zcXEp8146dOigLNAv2gaRulhsyCgdPHgQ1atX\nx2+//YZff/0Vv/76K3777Tf07t0b27dvVz4vLi4Ou3btwq1bt7Bs2TK4u7ujWrVquHDhAmJiYvD0\n6VNYWFigevXqygPwiBEj8MUXXyApKQmnTp3Cd999h5EjR5bJoXbt2nB2dkZwcDDS09Pxn//8p9SI\noWvXrrh27RoePXqkvG/SpElYtGgRYmNjUVRUhNTUVBw/fly239OIESPw3Xff4eTJk0hOTsYXX3yB\nUaNGyRaPTBeLDRml7du3Y+rUqWjSpAlsbGxgY2ODBg0a4IMPPsDu3btRWFgIhUKB6dOn4/vvv0eX\nLl3QuHFjrFu3DgCQlZWF6dOnw8rKCv369UO3bt0wceJEAMC8efMwcuRIjB49GitWrMDXX3+Nvn37\nAih77c7KlSvxyy+/oEOHDigqKoKLi4vyMQcHB4wcORLt2rWDjY0NAGDatGmYOnUqFi1aBCsrKwwa\nNAhJSUnlvkeFQqE8m63kz8WLF1/4e3n+uiJXV1esXr0an3/+OUaOHIkRI0bA39//hc8nqiyFKO9b\nSiIiIg3iyIaIiGTHYkNERLJjsSEiItmx2BARkexM5gT6N998E9evX9d1GkREBsXOzg7JyclV3o7J\njGyuX78OIYTWfvr27ct4BhrPmN8b4zGeuj+a+pBuMsVG25o3b854BhrPmN8b4zGerhj8NNqTJ0+w\nYMEC5ObmYsSIEXB3d9d1SgCMfwc05njG/N4Yj/F0xeBHNj///DOcnJwQGBiI/fv36zodJVdXV8Yz\n0HjG/N4Yj/F0RS9XEJg6dSpCQ0NhY2NTauHCqKgoeHt7o6CgAL6+vvDx8cHatWvRo0cPdOvWDZMm\nTcKOHTvK3aZCoYAevlUiIr2mqWOnXo5svLy8EBYWVuZ+Pz8/bN68GREREdi4cSPS0tLQsWNH3Lhx\nAwBQq1YtbadKREQq0Mti07t3b9SrV6/UfZmZmQCAPn36wNbWFoMHD0ZMTAx69uyJCxcuwMfHB6NH\nj9ZFukRE9BIGc4LA+fPn0bZtW+VtBwcHnD17FsOGDcOXX36p0jZcXV3RvHlzNG/eHK6urgY790lE\nJJdTp07h1KlTSElJQUpKisa2azDFRhNOnTql6xSIiPTa8x/ENdVmQi+n0crj5OSEq1evKm8nJibC\n2dlZhxkREZGqDKbY1K1bF4B0RlpKSgrCw8PZG52IyEDoZbHx8PBAz549kZSUhKZNmyIkJAQAsGbN\nGnh7e2PgwIGYNWsW6tevr+NMiYhIFXp5nY0ceJ0NEZH6jPo6GyIiMi4sNkREJDsWGyIikh2LDRER\nyY7FhoiIZGcUKwgcOnQIoaGhKCgowIwZM9CtWzddp0RERCUY1anP9+7dw+LFixEYGFjmMZ76TESk\nPqM89Xnq1Klo0KABOnToUOr+qKgo2Nvbo1WrVli/fv0LXx8QEABvb2+50yQiIjXpVbFRp4/Njh07\nMGfOHNy6dQtCCMybNw9Dhw6Fo6OjDjInIqKK6NV3Nr179y6zpHXJPjYAlH1sJk2ahEmTJgEA1q1b\nh5MnTyI7OxvJyckc3RAR6Rm9KjblqaiPTTFfX1/4+vrqIj0iIlKB3hcbTWLzNCKiipls8zQnJyf4\n+/srbycmJsLd3b1S22LzNCKiipls8zT2sSEiMnx6VWzYx4aIyDgZ1UWdFeFFnURE6jPKizqJiMg4\nsdgQEZHsWGyIiEh2LDZERCQ7FhsiIpIdiw0REcnOaIpNTk4OnJycEBoaqutUiIjoOUZTbFatWoVx\n48bpOg0iIiqHXhWbyjZPCw8Ph4ODA6ytrbWVKhERqUGvVhA4ffo06tSpg8mTJyMhIUF5f+fOnbF2\n7VrY2trCzc0N0dHROHr0KGJjY+Hv749NmzYhJycHV65cgYWFBQ4cOFBm8TiuIEBEpD5NHTv1atXn\nyjZP++yzzwAA27Ztg7W1tcZWKSUiIs3Qq2JTHlWapxXz9PSscFvsZ0NEVDGT7WejSexnQ0RUMZPt\nZ+Pk5ISrV68qbycmJsLZ2VmHGRERkbr0vtiweRoRkeHTq2LD5mlERMZJr059lhNPfSYiUh+bpxER\nkcFgsSEiItmx2BARkexYbIiISHYsNkREJDujWUFgw4YNuHHjBhwdHTF58mRdp0NERCUYxcgmLi4O\nx44dg5mZGezt7XWdDhERPUevik1l+9lER0ejX79+WLVqFTZt2qStdImISEV6VWy8vLwQFhZW5n4/\nPz9s3rwZERER2LhxI9LS0rBjxw7MmTMHt27dQseOHWFlZQWFQoHCwkIdZE5ERBXRq+9sKtvPpn79\n+ggPD8fcuXPLbT1ARES6pVfFpjyq9LN55ZVXlA3UiIhI/+h9sdEkNk8jIqqYyTZPc3Jygr+/v/J2\nYmIi3N3dK7UtNk8jIqqYyTZPYz8bIiLDp1fFhv1siIiME/vZEBHRC7GfDRERGQwWGyIikh2LDRER\nyY7FhoiIZMdiQ0REsmOxISIi2en9CgKqyMzMxOzZs2FpaQkHBwdMmzZN1ykREVEJRjGyiYmJQY8e\nPbB69WpEREToOh0iInqOXhWbyjZPc3Z2xq5duzBgwAAMGTJEW+kSEZGK9GoFgdOnT6NOnTqYPHky\nEhISlPd37twZa9euha2tLdzc3BAdHY2jR48iNjYW/v7+CA0NRZMmTTBkyBCMGTMG+/btK7NtriBA\nRKQ+TR079eo7m8o2Txs8eDCWLl2KkydPolu3blrNmYiIXk6vik15VGmeZmtri61bt750W+xnQ0RU\nMZPtZ6NJ7GdDRFQxk+1n4+TkhKtXrypvJyYmwtnZWYcZERGRuvS+2LB5GhGR4dOrYsPmaURExkmv\nTn2WE099JiJSH5unERGRwWCxoVIuXrwIDw8PDBs2DEOHDsWxY8d0nRIRGQGTOvWZKrZy5Urs3LkT\n+/fvR+vWrXWdDhEZEX5nQwCAiIgIDBkyBHFxcWjfvr2u0yEiPaGpYyeLDQGQFjO9ceMGXFxclPeN\nGDECU6ZM0V1SRKRzLDZqYrF5sfz8fNStWxezZs3CV199pet0iEiPmOzZaDdv3sR7772HsWPHAgDy\n8vIwd+5czJw5E2FhYTrOzjDl5+ejoKAAjRo10mked+/exciRI2FnZ/fC51y+fBmOjo5azIqINMHg\nik2LFi3wzTffKG+fOXMGTk5OCAwMxP79+3WYmeGqU6cOOnTogNTUVJ3m0aBBAwwfPhzu7u4vfE7r\n1q1x+PBhLWZFRJqgs2JT2UZpz0tISFB+Es7NzZUlV1OwdOlS7N+/H+np6cr7Tp48qfXf6bFjx+Dm\n5vbCx1955RU0a9ZMixkRkSbo7NRnLy8v+Pj4YPLkyaXu9/Pzw+bNm5WN0jw8PEo1Snt+qqdjx464\nceMGunXrhlq1amnzLRiV4cOHo6ioCP7+/qhfvz6aNm0Ke3t7WFhYID09HZs3b8aFCxewcOFCXL58\nGQ8fPsSjR4/Qr18/XL58Gbdv30azZs3g6ekJABBCYPXq1Xj99ddx584dJCcnY8GCBWjZsiUA4L//\n/S8ePHiAZs2awcLCAqNGjUJRURGioqLg5+eH48ePY8eOHfD29kavXr1QVFSErVu34vz58/Dz84OD\ngwPS09MRFBSE2NhYLFy4EImJifj7779x7949fPnll6XeX1BQECwsLGBubo4nT54o8yQiLRE6dPPm\nTdG+fXvl7YyMDOHo6Ki87ePjI44cOVLqNQ8ePBDe3t7Czs5OrFy5UuTl5YkPP/xQfPDBByIsLOyF\nsXT8Vg1acHCwePLkiWjbtq3YsWOHEEKIx48fCwsLC3H8+HEhhBBxcXGia9euytf8+9//FmvWrBFC\nCJGWlibq1aunfGzRokVi/vz5QgghDh06JFxdXYUQQpw7d040atRIpKamCiGEWLNmjfDz8xNCCLF/\n/35x79494enpKfbt21cmr++++04IIURmZqaoW7duqfwXLFggAgIChBBCZGVlic8//1yDvx0i46ap\nY6deXdSpSqM0KysrBAUFlXrd859iX4TN0ypn3LhxSE9PR3Z2NiZOnAgAiI+PR5cuXTBo0CAAQGxs\nrPLf7ubNm1izZg3u3r0LALh06RL69u0LALh+/To2btyIP//8EwAwYMAA5b/D8ePHMWvWLOXo9eLF\ni8oOrQMHDoQQAidOnMDmzZtL5ZWTk4Px48crX1PyOqEHDx5g3bp1CAwMxO7du2FmZgZfX1/ZfldE\nho7N0zSAzdMq57XXXsPhw4cxYMAA5X0nTpzAwIEDlbf37NmDmTNnIiMjAydOnICzszNq166tfG7/\n/v2RkZGB8PDwUlOexc8BpAtL//Of/wCQvn8LCwvD119/jYyMDFhaWiIwMBCjRo1CYWEhCgoKys1r\n7969GDNmDLKysvDaa6/h559/hr29vbKFOBFVzCSap7FRmv6KiIgodVAvefvhw4eIjY3FsGHDsG/f\nPjRu3BhNmjQBADx69Aj79+9Hz549sXfvXjRq1Ah16tRRbqegoAB79uxBTk4OEhMT0a1bNwDA0aNH\n0b17d1hYWGDv3r0AgF27dmHy5MnYs2eP8g/g+bz27t2L8ePH43//+x8A6ey1oqIi5eNPnz5VjoyI\nSHv0amRTslFas2bNEB4ejsWLF+s4KwKA5ORkBAQEAJC+/E9LS0OPHj0ASP9uffr0wa5duzBgwAA0\na9YMZ8+exd69e9G4cWN4e3tj3759GDJkCPr06YPz589j06ZNsLS0RG5uLsaOHYuEhAS4u7ujWjXp\n80/79u1Rr149bNiwAXPmzAEA9O3bFzExMejatSvMzMzK5AUAb731FkJDQ5VFq23btvD09MQXX3yB\nxo0bIy8vTzkVSETao7MVBDw8PBAZGYkHDx7AxsYGy5Ytg5eXFyIjIzFjxgw8ffoUvr6+Gptf5woC\nRETq43I1amKxISJSn8kuV0NERIaHxYaIiGTHYkNERLJjsaEy/v77b3zwwQc4cOCArlMhIiPBYkNK\nd+7cwezZs9GpUyfUrl2bKywQkcYYZLF5vqfNoUOHMH36dEydOhXnzp3TcXaGJy0tDfPmzUO7du2g\nUChw5coVBAQEoF69erpOjYiMhEGf+jx27Fjl1eUAcO/ePSxevBiBgYFlnstTn8t6+PAhvvrqKwQG\nBmLcuHH45JNPlFf+ExEBmjt2qryCQH5+Pn744Qf88MMPKCgogJmZGbKzs2FlZYXBgwdjzJgxyqu/\nVTV16lSEhobCxsYGCQkJyvujoqLg7e2NgoIC+Pr6wsfHR6XtBQQEwNvbW60cTFFWVhbWrl2LdevW\n4e2338bFixfRvHlzXadFRMZMlaWho6OjxYIFC0RiYqLIz88v9djjx4/FuXPnhK+vr7h48aJaS05H\nRUWJ2NjYUm0GhBDC0dFRREZGipSUFNGmTRtx//59sX37djF79mzl8vNCCDFmzBghhBBFRUXC399f\nREREvDCWim/VqD169EisXLlSWFtbi4kTJ4pr167pOiUi0nOaOna+dGTz5MkT1K1bF59//nm5j1tY\nWMDJyQlOTk44f/68WoWud+/eZZawzszMBADl0vKDBw9GTEwMJk2apFy5Nz09HZ988gni4+OxcuVK\n1K5dGydPnkR2djaSk5M5unlOXl4egoKCEBAQgD59+iAyMhL29va6TouITMhLi425uXmp/iA//vgj\nhg4dij/++AO2trZISkpC69atAUirNldVZXvaqDLVZmr9bPLz87FlyxasWLEC//jHPxAWFoZOnTrp\nOi0i0mN608/m3LlzGDp0KE6dOgVPT0/8/fffymKj70yln83Tp0+xfft2LF++HA4ODjh48CD+8Y9/\n6DotIjIAcvWzUbvYdOrUCf3790fDhg1Rr1493Lx5E/3799dIMoA0OvL391feTkxMhLu7u8a2b8wK\nCwuxe/duLF26FLa2tti9ezd69uyp67SIiNQvNqNGjUK7du0QFhaGjIwMzJgxQ6MJsaeN+oqKirBv\n3z4sWbIEr7/+Or755hujnyIkIsPy0utsnjx5guzsbNSvX/+lG7t+/Trs7OxUDq7NnjbGeJ2NEAKH\nDh3C4sWLYW5ujuXLl2Pw4MEaG/YSEWm1n82RI0eQlZWFUaNGwcLCoszj9+/fx5o1a+Dq6opBgwZV\nOSk5GFOxEUIgLCwMCxcuRGFhIZYtW4a33nqLRYaINE7rzdNu376NkJAQ3Lt3D3l5ecjLy0NWVhZq\n1qwJR0dHeHt7K6fA9JGxFJuTJ0/i008/RWZmJpYtW4ZRo0apfTEtEZGqtL6CgLm5OTp16gQLCwt0\n6dIFlpaWVQ5OqouOjsbChQuRmpqKJUuWYNy4cTAzM9N1WkREKlFpZBMaGoqQkBDY2dnh1q1bOHfu\nHNzc3LB8+XK9Hs2UZKgjm3PnzmHhwoVISkrC4sWLMXHiRFSvrvZ5HURElaLVkc1vv/2Gffv2KW8L\nIXD8+HFMmTIFu3btQq1ataqcCJUWFxeHRYsWIT4+Hp9++im8vLzwyiuv6DotIqJKUWmy38HBodRt\nhUIBNzc3BAcHl7mSn6omMTERY8aMwbBhwzB48GBcu3YN3t7eLDREZNBUKjaPHj1CXl5emfutra1h\nY2Oj8aRMUVJSEv75z3+if//+cHZ2RnJyMnx8fFCzZk1dp0ZEVGUqFZtevXph1qxZiI+PL/PYpUuX\nNJ5URZ5vnAYAOTk5cHJyQmhoqFZz0YSbN2/Cy8sLLi4uaN++PZKTk/HRRx9xapKIjIpK39k0atQI\nq1evxoQJE5CamgpnZ2fUqlUL58+f1/oKyy1atMA333xTqtisWrUK48aN02oeVfXXX3/hs88+w/ff\nf4/3338f165d4xl+RGS0VL5Ao27dujhy5AhCQ0PRq1cvtG7dGv/9738xYcKESgWeOnUqGjRogA4d\nOpS6PyoqCvb29mjVqhXWr1//0u2Eh4fDwcEB1tbWlcpD227fvg1fX184OjrCysoKv//+O5YuXcpC\nQ0RGTe1zaBs1alTpAlOSl5cXfHx8MHny5FL3+/n5YfPmzbC1tYWbmxs8PDxw9OhRxMbGwt/fH40a\nNSr1/MjISOTk5ODKlSuwsLDA0KFD9fJK+vv37yMgIAAhISHw9PTElStX0KBBA12nRUSkFTq7YEMT\njdMCAgLw2WefAQC2bdsGa2trvSs06enp+OqrrxAUFAQPDw8kJCSUKZhERMZOr64OrGzjNADw9PR8\n6fa12TwtKysLq1evxvr16zFq1CjExsbC1tZWtnhERJqgN83TDJk2mqc9evQIGzZswNdff40hQ4Yg\nJiZGrZWwiYh0Sa7maXq1gqOTkxOuXr2qvJ2YmAhnZ2cdZqS63NxcfP3113jzzTcRHx+PqKgobNu2\njYWGiAh6VmxKNk5LSUlBeHg4unfvruOsKvbkyRNs3LgRb775JqKjoxEeHo7//e9/paYDiYhMnc6K\njYeHB3r27ImkpCQ0bdoUISEhAIA1a9bA29sbAwcOxKxZs1Rq2qYLT58+xTfffIPWrVvjxx9/xOHD\nh7F///4yp3ITEZEa/WwMnaZWLi0sLMSuXbuwbNkytGjRAsuWLUOPHj00kCERkf7Rej8bY/fw4UPU\nq1fvhY8XFRVhz549WLJkCWxsbLBlyxb07dtXixkSERkuFhsAZ8+exejRo5GamlrmzAshBA4ePIjF\nixejVq1aWL9+PQYOHKh31/MQEekzky82GRkZ8PDwQGBgYKkCIoTAjz/+iEWLFkEIgc8//xzDhg1j\nkSEiqgST/s5GCIGxY8eiYcOGynXYhBA4ceIEFi5ciOzsbCxbtgwjR45EtWp6deIeEZFW8DsbDQgO\nDkZycjJ27twJQDrleuHChbhz5w6WLFmCd999F2ZmZjrOkojI8BncyObmzZtYsWIFMjMzsXfvXgDA\nhg0bcOPGDTg6OpZZ2LPY89U5ISEB/fv3R3R0NDIyMrBw4UIkJydj8eLFmDBhAqpXN+k6TEQEQHMj\nG4ObGyruZ1MsLi4Ox44dg5mZGezt7VXaRk5ODsaNG4cPPvgAH374IcaOHYsxY8bg999/h6enJwsN\nEZGG6azYaKqfzenTp9GvXz+sWrUKmzZtUim2p6cnHj16hKCgIAwcOBAJCQmYNGkSnjx5gqysLDx8\n+BBPnz6t1PsiIqKydFZsvLy8EBYWVub+4n42ERER2LhxI9LS0rBjxw7MmTMHt27dKvP8Tp06wcrK\nCgqFAoWFhS+NGxcXh++//x5//fUX7ty5g48++gjW1tawsrJCgwYN0KRJE7Ro0QJLlizRxNskIiIY\ncD+buLg4BAQEYM6cOQgPD8fcuXNLtSJ4kc6dO6OgoAAKhYJnmBERaYlefTlR2X42xQ3UXkab/Wwq\n4/FjwMIC0OdLeXJzgTt3gKZNAX39ais/H/jrLyA9HWjXDqhVS9cZlVVQAKSmAn/8If27OzkBr7+u\n66yI2M9GI6ysTqFTJ6BPH6BTJ+3GFgJIS5MOLi/6yckBXnsN6NVL+undG+jcGahRQ3t5ZmRUnGNm\nJlC/vvTf7t2f5dm9O1CnjnZyfPSo4hzv3wcaNQIsLYFr14AOHZ7l6eIi5S+33Fzgzz9fnOPt24C1\nNWBrC9SsCVy8CDRpIuVY/O9va6vfHzzIOMnVz0anpz6npKRg+PDhSEhIACBNo7m6uiIuLg4A4OPj\nA3d3d5Wmx15GoVBg926BqCggKgr4+2+gRw+p8PTuLX2yrFmz8tsvLARu3XrxweXPPwFzc+kA8qKf\n+vWlvKKjpZ/Tp4GbN4Fu3Z4dgHr0qPxBvagIuHu34hyFKJtXs2bP/v+NN4Bq1aRRw5kzUo7R0UB8\nvDSKKM6zVy/Axkb9HIUAHjyouJjk5pbO6fmfRo2ejboePwbOn3/2+/zlF6Bx49IFvXlz9Q/qLyvK\nGRlS8XhRjk2bAq+88mx7BQXApUvP8jx9WvqQUTLPdu0AdS77yskpP7eCAmk/6tULcHTU7ocZMjya\nOvVZr4oNIH2nsnbtWjRr1gzu7u6Ijo7WSJuB539haWnSH3ZUlPSH/dtvQNeu0h91nz7SH+Orrz57\nffHUzB9/ACkpZf+AU1OlaZCKiknJ7anq4UPpAFl8UI+NBeztnx2AevUCGjSQnvv0qVSsXnQA/Osv\naeRUUY6WlpX7NJ2b++ygHh0tFaIGDUrnaWcnFbyqFOVmzaQRQWU/bBUWPjuoFx/YFYrSebZv//JR\naFFRxUW5YUOpKFeWEMD166XzvHcP6NnzWQFq3frZVNyLRsrl/Q4VCuDnn6XtpqQ8+zCj7REqGQaD\nLzYeHh6IjIzEgwcPYGNjg2XLlsHLywuRkZGYMWMGnj59Cl9fX/j6+mok3st+YdnZ0kG9eOQTGwu0\nbSt9+iw5I2d6AAAQ1klEQVQ5NdO8edk/3ubNpU+q5uYaSbVCeXnAhQvPDkBnzgBWVlIxvHtXGnlU\ndKDW1vcXhYXA5cvPiuTp09Io4/FjeYpyZQkhjR6Lc4yOlqberKzKFpCSP/XqaX+K6949qUgU53n9\nesWjJxubl+f48GHpEWpcHODg8Kzwurg8+zBDpsngi422qfsLKz6oF08rlZya0SdFRcDvv0snFjRu\nrL9TIkJIB0tLS+0U5aooKNDPf2ttyMsrPe1YcoQ6ciTw1lv8HsnUsNioSVO/MCJTUlgIJCZKhSco\nSBo5r1snTeWSaTDZ5WqISHvMzICOHYH335em2IYPl77T/OgjICtL19mRIWGxISKVVK8O+PpK38Wl\np0ujm507pSlSopfhNBoRVcrZs8AHH0iXDGzYIJ1GTcaH02hEpFPOzkBMDODpCbi5SVNt6em6zor0\nFYsNEVWamRkwbZp0nZpCIU2tBQdLJxYQlWRw02jPN0/LzMzE7NmzYWlpCQcHB0ybNq3c13EajUh+\n8fHS1FpenjS15uys64yoqkx2Gu355mkxMTHo0aMHVq9ejYiICB1mRkSOjtJp0rNnA++8A3h5SRcb\nExl88zRnZ2fs2rULAwYMwJAhQ+RKl4hUpFAAEydKU2v160vL/6xdKy2nRKbL4Jun/d///R/mz5+P\nEydO4MiRI9pInYhU8NprwBdfSMs/HTkirWD+00+6zop0xWCbp8XHxyMgIADjx4/H0qVLcfLkSXTr\n1k2r74GIXs7eHjh+HDhwQJpWc3YGvvxSWteNTIderQBV2eZpW7duVWn7+t48jchYKRTA6NGAuzsQ\nECB9t/Phh8Dcufq/Vp6pYfM0DTh16pSuUyAyabVqAUuXStfmzJ377PucoUN1nRkVk6t5ml6djebk\n5ISrV68qbycmJsKZ504SGZ2WLYGDB6VFPWfPBt5+W2qZQMZLr4pN3bp1AUhnpKWkpCA8PBzdu3fX\ncVZEJJchQ4CEBKlvTvfuwMKFUs8jMj46KzYeHh7o2bMnkpKS0LRpU4SEhAAA1qxZA29vbwwcOBCz\nZs3SSJdOItJf5ubAxx9LF4QmJ0snFOzbxwU+jY3BrSBQWVxBgMgwnDoF+PhITdvWrZM6h5LumOwK\nAkRk3Fxdpd45I0YAffuyd46xYLEhIr1Tvbo0uklMBB4+BNq2BXbs4NSaIeM0GhHpvZgYaYFPc3Ng\n/XppNQLSDk6jEZHJ6N5dKjhTpkgXhs6axd45hobFhogMQrVqwHvvSQt8VqsmnbW2eTN75xgKg5tG\nO3ToEEJDQ1FQUIAZM2agU6dOWLBgAXJzczFixAi4u7uX+zpOoxEZl19/labWHj+Weuf06KHrjIyT\npo6dBldsit27dw+LFy/Gu+++izt37sDDwwPTp09HcHBwuc9nsSEyPkIAu3ZJ1+kMHgysXCmdMk2a\nY/Df2VS1n01AQAC8vb1x6dIl2NnZAQByc3NlzZmI9MvzvXPatQPWrGHvHH1kcP1shBCYN28ehg4d\nCkdHR3Ts2BE3btwAANSqVUvbb4OI9EBx75zTp4HQUPbO0UcG189m3bp1OHnyJLKzs5GcnIwpU6bg\n3//+N37++WeMHj1aq++BiPTL871zuneXeuc0barrzEivWgyo0s/G19cXvr6+pV735ZdfqrR99rMh\nMn4v6p3z4YfsnaMK9rPRAPazITIdJXvnzJnD3jmqYj8bIqJKaNkSOHSIvXN0Ta+KDfvZEJFcinvn\n9OwJdOvG3jnaxn42RGQyzM2B+fOlC0KLe+d8/z0X+NQGg72oU128qJOInleyd8769VLxodIM/qJO\nIiJdK+6dM3w40KcPe+fIicWGiExa9eqAnx9w+bK0krS9PXvnyIHTaEREJZw9Ky3wWbMme+cAnEYj\nIpKFs7PUO8fTk71zNInFhojoOWZmwLRp7J2jSZxGIyJ6ifh4aWotN9f0eueYdD+b5xuo3b59u9Tt\nbt26lXkNiw0RVYWp9s4x6WJTrLiBWmBgYLm3S2KxISJNyMoCli8HQkKATz8F3n8fqFFD11nJxyhO\nENBUA7UX3SYi0rSSvXN+/JG9c1Sl05HN6dOnUadOHUyePBkJCQnK+zt37oy1a9fC1tYWbm5uiI6O\nxtGjRxEbGwt/f380bNgQH3/8Mdzc3DBgwAAIIUrdLg9HNkSkaUJIvXPmzjXe3jmaOnbqtMWAphqo\n5efnl7rN0Q0RacPzvXM6d5YKD3vnlKV3/Wwq20DNx8fnpdtm8zQikkPJ3jlz5xp27xw2T9MANk8j\nIjm1bAkcPAgcPSotgRMUBKxeDdjZ6Toz1ZlE8zSADdSIyPAV985xcZG+y2HvHD0sNmygRkTGwNxc\nuiYnPv5Z75x9+0x3gU+dFhs2UCMiY9ekCfDdd8C2bdL3OoMGAVeu6Dor7TPoizrVwVOfiUjXCgqA\nTZuki0InTwYWL5au29FnRnFRJxGRKaleHfD1BRITgYwMoG1b0+mdw5ENEZGOxMRIC3yam+tv7xyO\nbIiIDFz37lLBmTLF+HvnsNgQEelQtWrAe+8Zf+8cTqMREekRfeudY7LTaIcOHcL06dMxdepUnDt3\nDgCQk5MDJycnhIaG6jg7IqKqcXSUVpSeMwcYMwbw8gLu3tV1VlVncMVmxIgRCA4OxsqVK5XX5axa\ntQrjxo3TcWZERJqhUAATJ0pTa/XrA+3aAWvWAE+f6jqzytNZsdFUL5vw8HA4ODjA2tpa7pSJiLSq\nZO+c0FDD7p2js+9sNNXL5tNPP0VOTg6uXLkCCwsLHDhwoNyF4/idDREZMl31zjH4fjaa6mXz2Wef\nAQC2bdsGa2trja1QSkSkT57vnePoKPXNMZTeOXrVYqCyvWwAwNPT86XbZz8bIjJ0JXvnzJmj+d45\n7GejAexnQ0TGomVL4NAhzffOMYl+NuxlQ0SknuLeOT17St/lLFqkn71z9KrYsJcNEZH6zM2B+fOl\nC0KvXZNWIfj+e/1a4FNnxYa9bIiINKtk75wlS4DBg6VrdfQBl6shIjJCJXvneHpK02uV6Z1jssvV\nEBHRyxX3zrl8WVpJ2t4e2LlTd1NrHNkQEZmAs2elBT5r1pQW+HR0VO11HNkQEZHKnJ2l3jmenoCb\nG/D++9rtncNiQ0RkIszMgGnTpJMGFAppai04WDu9cziNRkRkoop75+TlSVNr5V3WyGk0IiKqkuLe\nObNnA++8I2/vHIMrNuU1T9uwYQPmzp2L7du36zg7IiLD8nzvnOK11jTdO8fgis3zzdPi4+Nx7Ngx\nmJmZwd7eXtfpKWl7HTbGM8xYjMd4+hKvuHdOVBRw5Ijme+cYfPO006dPo1+/fli1ahU2bdokd9oq\nM5Yd0BTjGfN7YzzGexl7e+D4cWDZMmDKFM1tV2fFxsvLC2FhYWXu9/Pzw+bNmxEREYGNGzciLS0N\nO3bswJw5c3Dr1i0IITBv3jwMHToUjo6O6NixI6ysrKBQKFCojVMqVKTJpbkZT7vxjPm9MR7jqaK4\nd05Skua2afDN07y8vBAeHo65c+eW6nuja8a4A5pKPGN+b4zHeOrQZFM2vepnU9nmacXdOivy+uuv\na72LJ+MZbjxjfm+Mx3jqeP311zWyHb0qNnJKS0vTdQpERCZLr85GY/M0IiLjpFfFhs3TiIiME5un\nERGR7ExmbTQiItIdvZpG06abN2/ivffew9ixY2WPdfXqVcycORP/+te/sH//ftnjnTp1Cr1798bM\nmTMRGRkpe7zo6GjMnDkT06ZNg4uLi+zx/vrrL8yYMQPLli1DRESEbHHK20fk3G/K27ac+0558eTc\nd8qLJ9e+U14sOfeb8pbRknNfKS+enPtKefHU3leEiRszZozWYj158kSMHz9e9jiRkZHC3d1d+Pn5\niT///FP2eMUOHjwogoODZY+zdetWERYWJoQQwtPTU/Z45e0jcu435W1bzn2nZDxt7DvlvT+59p2S\nsbSx39y9e1fMmDHjhTloI56c+0rJeOruK0Y1sqnqEjhyxjp8+DD69euHd999V/Z4vXv3xtGjRzF7\n9mx8+eWXsscrtnv3bvzzn/+UPd4777yD6OhofPzxx0hOTpYtjiZoIp46+05V46m772jq96nKvlPV\nWOruN5WJV7yMVmVoIp7c+0rJeGofZ2QpfzoSFRUlYmNjRfv27Uvd7+joKCIjI0VKSopo06aNuH//\nvvKxyn7qqEwsIYQYPny41uJlZmaKf/3rX1qJ98cff4hp06ZVKlZl4gkhRH5+vvDy8pI9TlVGNprc\nJ1XZdzQVT9V9RxPxVN13NPXeVN1vVI2XlpYmioqKhL+/v4iIiCizHU3vKy+LJ4Rm95WXxVN1XzGq\nizrVWQKnR48e+OSTTxAfH4+AgAB8/PHHssWqU6cO9u/fDyFEpedv1YmXn5+PY8eOoaCgADNnzpQ9\n3rBhw7B161ZMnTq1UrHUjde+fXusWLECCoUCfn5+ssUpbx9JT09Xa7+pTLy4uDjltiMjI9Xad6oa\n78CBA2rtO5r4m1N136lqrD/++EOt/UbVeGfPnsWNGzdKLaPl7e0t277yonhy7SsviqfuvmJUxaY8\nFS2BExQUpJVYy5cvR9++fTUa62XxRo0apbV4w4YNw5IlS7QaLzg4WCtxnt9HrKysqrzfqBOvb9++\nVd531Ik3atSoKu876v7NVWXfUSeWra1tlfebiv7mfHx8Sj1Xzn2lvHhy7ivlxVN3XzGq72yIiEg/\nGX2x0eYSONpebofxDCsO4xl2LMarGqMvNtpcAkfby+0wnmHFYTzDjsV4VfTSUwgMyPjx40XDhg3F\nK6+8Ipo0aSK2bt0qhBDi1KlTom3btsLOzk6sXbvW4GIxnubiGev7MoV4xvzeTCEel6shIiLZGf00\nGhER6R6LDRERyY7FhoiIZMdiQ0REsmOxISIi2bHYEBGR7FhsiIhIdiw2REQkOxYbIiKSHYsNkY49\nefKkUo8RGRKj72dDpGu3bt1CdHQ06tSpgytXrsDe3h7Dhg0DABw5cgTOzs4wNzcv97XXrl1DZmYm\nXFxctJkykcZxZEMko8zMTHz66ad4++23MXToUHz00UcICQlBeno6bt++jaysLNSvX/+Fr2/fvj1C\nQ0ORn5+vxayJNI/FhkhGEyZMgK+vL2rWrKm8z87ODr/88gtCQkJU6nQ4ceJE/PDDD3KmSSQ7Fhsi\nmdy6dQupqalwdHQsc7+FhQXu3bsHCwsLAEBhYSH27t2LFStWYMuWLZgxYwZu3LgBAGjVqhWLDRk8\nFhsimfz0009luhwKIXDu3Dl06dIFeXl5yvt//fVXvP3227C1tUW1atUwfvx4NGzYEABQo0YNFBQU\naDV3Ik1jsSGSyaNHj1CrVq1S9x0/fhxubm6wtLQsVWy6dOkCc3NzxMTEwNXVFa6urspRDwCYmZlp\nLW8iObDYEMnE1dUV58+fV96+f/8+goKCsHz5cgBAVlaW8rHz588jLS0Nly9fRosWLRAdHV1qW9nZ\n2dpJmkgmPPWZSCZt2rSBt7c35s2bh44dO+LevXvYtWuXcrRT8qSBsLAwWFpawsXFBQcPHkTjxo2V\nj+Xm5sLKykrr+RNpEttCE+nIqlWr4OXlBWtr6wqfd/78edy8eRPvvvuuljIj0jxOoxHpiLe3N9as\nWfPS5+3cuRNjxozRQkZE8uHIhkiHwsPD0bJlS9jZ2ZX7+Pnz52FmZoYuXbpoOTMizWKxISIi2XEa\njYiIZMdiQ0REsmOxISIi2bHYEBGR7FhsiIhIdiw2REQkOxYbIiKSHYsNERHJ7v8BR+mYcGbz54EA\nAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 19 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Expressed this way, the relative error is constant up to $10^{16}$, in the order of magnitude $10^{16}$.\n", "\n", "The explanation for these observations is the following: floating point numbers are represented to a *relative* precision, denoted by \n", "\n", "$\\epsilon_{machine} = 2.220446049250313 \\times 10^{\u221216}$ for 64bit floating point numbers\n", "\n", "It represents the smallest *relative* gap between representable numbers, so that, for example, all numbers between\n", "\n", "$[1,1+\\epsilon_{machine}]$\n", "\n", "are rounded to either one.\n", "\n", "We can now see why the absolute error of $(10^{17}+0.2) - 10^{17}$ is *exactly* $0.2$: we cannot represent the absolute value of $0.2$ at this magnitude, since the absolute gap between representable number is \n" ] }, { "cell_type": "code", "collapsed": false, "input": [ "1E17*2.220446049250313E-16" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 20, "text": [ "22.20446049250313" ] } ], "prompt_number": 20 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This gap, in absolute terms, is much larger than $0.2$, the value we're trying to subtract." ] }, { "cell_type": "code", "collapsed": false, "input": [ "1e17 - 0.2" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 21, "text": [ "1e+17" ] } ], "prompt_number": 21 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is also referred to as catastrophic cancellation, when there is devastating loss of precision when small numbers are computed from large numbers, which themselves are subject to round-off error. So is" ] }, { "cell_type": "code", "collapsed": false, "input": [ "1e17 - 2" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 22, "text": [ "1e+17" ] } ], "prompt_number": 22 }, { "cell_type": "markdown", "metadata": {}, "source": [ "but not" ] }, { "cell_type": "code", "collapsed": false, "input": [ "1e17 - 20" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 23, "text": [ "9.999999999999998e+16" ] } ], "prompt_number": 23 }, { "cell_type": "markdown", "metadata": {}, "source": [ "since here the subtraction is *significant* enough to cause a round-off to $10^{17}-10^{17}\\times\\epsilon_{machine}$\n", "\n", "To summarize, computers use a floating point number representation, where **range** *and* **precision** are limited by the allocated memory (bits). The precision is relative and given by $\\epsilon_{machine}$. A nice illustration of the absolute spacing between representable floating point numbers is shown below\n", "\n", "\n", "\n", "for a reduced system with only three bits for the fraction and exponent. The numbers on the x-axis are powers of two, and the red lines indicate the scaling of gaps between representable numbers with increasing magnitude. The spacing is the absolute value of $\\epsilon_{machine}$ in that region." ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Hexadecimal (base 16) numeral system" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another numeral system encountered in computer engineering is the hexadecimal, short hex, system. In very simple terms, it extends the deximal system\n", "\n", "$0, 1, 2, 3, 4, 5, 6, 7, 8, 9$ by $A, B, C, D, E, F$ (case insensitive), to yield\n", "\n", "$0, 1, 2, 3, 4, 5, 6, 7, 8, 9, A, B, C, D, E, F$ (16 values)\n", "\n", "so that the 8-bit hex number $0001E078$ can be represented in exponent format as (starting with the left most non-zero bit)\n", "\n", "$(1\u2009\u00d7\u200916^4) + (E\u2009\u00d7\u200916^3) + (0\u2009\u00d7\u200916^2) + (7\u2009\u00d7\u200916^1) + (8\u2009\u00d7\u200916^0)$\n", "\n", "equal to (using decimal counterparts)\n", "\n", "$(1\u2009\u00d7\u200916^4) + (14\u2009\u00d7\u200916^3) + (0\u2009\u00d7\u200916^2) + (7\u2009\u00d7\u200916^1) + (8\u2009\u00d7\u200916^0)$\n", "\n", "to yield" ] }, { "cell_type": "code", "collapsed": false, "input": [ "1*16**4 + 14*16**3 + 0*16**2 + 7*16**1 + 8*16**0" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 24, "text": [ "123000" ] } ], "prompt_number": 24 }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Floating point (FP) arithemtic \u2260 arithmetic in ${\\rm I\\!R}$, but..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We revisit the starting point of this module:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "1.2-1. == 10.2-10." ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 25, "text": [ "False" ] } ], "prompt_number": 25 }, { "cell_type": "markdown", "metadata": {}, "source": [ "But we know now that it should be within some \\epsilon_{machine}. So the should not indicate whether numbers are bitwise equal (what the '==' operatore looks at), but whether they are the close enough to account for floating point errors. Thus, the equivalence of numbers in FP is checked for within a given *tolerance*.\n", "Fortunately, we're not the first ones to come across that (and chances are for most problems you face in the near future you aren't either). So there is a function out there for a more appropriate comparison:\n" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.allclose(1.2-1., 10.2-10., rtol=1e-016, atol=1e-15)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 26, "text": [ "True" ] } ], "prompt_number": 26 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This method returns $rtol \\times |b| +atol < |a-b|$ with $a=1.2-1.$ and $b=10.2-10.$." ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Problems" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Download the notebook by using the button in the nbviewer bar on the top left.\n", "\n", "1. \n", "The FP calculation $v^2$ results in $\\infty$ for what $O(v)$ in double-precision?\n", "\n", "2. After how many iterations do you expect \n", "```\n", "i = 0.\n", "while i != 1.:\n", " i = i+0.1\n", "```\n", "to stop?\n", "\n", "3. Does the order of floating point operations matter? Check for example \n", "```\n", "1e-16 + 1 - 1e-16 == 1e-16 - 1e-16 + 1 \n", "```\n", "\n", "4. Write a function that compares two numbers `a, b` to some relative tolerance `rtol`!\n", "```\n", "def is_equal(a, b, rtol):\n", " pass # here your code\n", "```\n", "\n", "5. What is the FP result of $\\sqrt{1e-16 + 1} - 1$?\n", "\n", "6. With respect to FP arithmetic, what would be an advantage of using dimensionless numbers?\n", "\n", "7. You store a surface map in the form of a two-dimensional array with `1E6 x 1E6` height measurement points. How much memory do you need for half-, single- and double-precision? If the measurement accuracy is in the order of 20 cm, and minimum/maximum height measurements are in the order of 100 meters, what precision do you recommend?\n", "\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**The FP calculation $v^2$ results in $\\infty$ for what $O(v)$ in double-precision?**" ] }, { "cell_type": "code", "collapsed": false, "input": [ "v = 1.0E154\n", "v**2 # largest representable number" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 34, "text": [ "1e+308" ] } ], "prompt_number": 34 }, { "cell_type": "code", "collapsed": false, "input": [ "v = 1.0E155\n", "v**2 # overflow" ], "language": "python", "metadata": {}, "outputs": [ { "ename": "OverflowError", "evalue": "(34, 'Result too large')", "output_type": "pyerr", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mOverflowError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[0mv\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m1.0E155\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mv\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[1;31mOverflowError\u001b[0m: (34, 'Result too large')" ] } ], "prompt_number": 35 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The largest representable number in double-presicion os of `O(1E308)`, so that `v = O(1E154)`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**After how many iterations do you expect**\n", "```\n", "i = 0.\n", "while i != 1.:\n", " i = i+0.1\n", "```\n", "**to stop?**\n", "\n", "Never - as has been shown, floating point operations are only accurate within some precision. The statement `while i != 1.` will never yield `False`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Does the order of floating point operations matter? Check for example **" ] }, { "cell_type": "code", "collapsed": false, "input": [ "1e-16 + 1 - 1e-16 == 1e-16 - 1e-16 + 1" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 36, "text": [ "False" ] } ], "prompt_number": 36 }, { "cell_type": "markdown", "metadata": {}, "source": [ "yes, order matters!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Write a function that compares two numbers a, b to some relative tolerance rtol!**" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def is_equal(a, b, rtol):\n", " return abs(a-b) < max(a,b)*rtol" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 47 }, { "cell_type": "code", "collapsed": false, "input": [ "is_equal(1.2-1., 10.2-10., 1.0E-14)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 48, "text": [ "True" ] } ], "prompt_number": 48 }, { "cell_type": "markdown", "metadata": {}, "source": [ "**What is the FP result of $\\sqrt{1e-16 + 1} - 1$?**" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.sqrt(1.0E-16 + 1.) - 1." ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 49, "text": [ "0.0" ] } ], "prompt_number": 49 }, { "cell_type": "markdown", "metadata": {}, "source": [ "**With respect to FP arithmetic, what would be an advantage of using dimensionless numbers?**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dimensionless numbers are mostly $O(1)$, so floating point operations involving numbers of large differences in order of mangnitude can be avoided. Also the limited range and precision of single-precision is mitigated." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**You store a surface map in the form of a two-dimensional array with `1E6 x 1E6` height measurement points. How much memory do you need for half-, single- and double-precision? If the measurement accuracy is in the order of 20 cm, and minimum/maximum height measurements are in the order of 100 meters, what precision do you recommend?**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Memory:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# numbers bits bits/byte\n", "bytes = 1E6**2 * 64 / 8\n", "kbytes = bytes/1024\n", "mbytes = kbytes/1024\n", "gbytes = mbytes/1024\n", "tbytes = gbytes/1024" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 62 }, { "cell_type": "code", "collapsed": false, "input": [ "print bytes\n", "print kbytes\n", "print mbytes\n", "print gbytes\n", "print tbytes" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "8e+12\n", "7812500000.0\n", "7629394.53125\n", "7450.58059692\n", "7.27595761418\n" ] } ], "prompt_number": 63 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Necessary precision:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "epsilon_half_16 = 4.88e-04\n", "epsilon_single_32 = 5.96e-08\n", "epsilon_double_64 = 1.11e-16\n", "order = 100.\n", "print order*epsilon_half_16\n", "print order*epsilon_single_32\n", "print order*epsilon_double_64" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "0.0488\n", "5.96e-06\n", "1.11e-14\n" ] } ], "prompt_number": 65 }, { "cell_type": "markdown", "metadata": {}, "source": [ "So that half precision suffices to represent numbers with an absolute precision of down to about 4 [cm], which is less than our absolte measurement precision of 20 [cm]. This way, we use " ] }, { "cell_type": "code", "collapsed": false, "input": [ "# numbers bits bits/byte\n", "bytes = 1E6**2 * 16 / 8\n", "kbytes = bytes/1024\n", "mbytes = kbytes/1024\n", "gbytes = mbytes/1024\n", "tbytes = gbytes/1024" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 66 }, { "cell_type": "code", "collapsed": false, "input": [ "print bytes\n", "print kbytes\n", "print mbytes\n", "print gbytes\n", "print tbytes" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "2e+12\n", "1953125000.0\n", "1907348.63281\n", "1862.64514923\n", "1.81898940355\n" ] } ], "prompt_number": 67 }, { "cell_type": "markdown", "metadata": {}, "source": [ "and save a fourth the memory." ] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }