{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Accuracy of operations involving floating point numbers" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 1: Equality test for two floats\n", "\n", "Equality tests involving two floating point numbers can be tricky\n", "\n", "Consider $x = 1.1 + 2.2$.\n", "\n", "The answer should be $x = 3.3$ but due to round-off error one can only assume $x = 3.3 + \\varepsilon_M$ where e.g. $\\varepsilon_M \\sim 10^{-15}$ is the machine precision for 64-bit floating point numbers.\n", "\n", "For this reason an equality test $x == 3.3$ might give some unexpected results..." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x = 3.3000000000000003\n", "3.3\n", "x == 3.3 is False\n" ] } ], "source": [ "x = 1.1 + 2.2\n", "\n", "print(\"x = \",x)\n", "print(3.3)\n", "\n", "if (x == 3.3):\n", " print(\"x == 3.3 is True\")\n", "else:\n", " print(\"x == 3.3 is False\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A safer way to compare two floats is to check the equality only wihtin a certain precision $\\varepsilon$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x = 3.3000000000000003\n", "x == 3.3 to a precision of 1e-15 is True\n" ] } ], "source": [ "print(\"x = \",x)\n", "\n", "# The desired precision\n", "eps = 1.e-15\n", "\n", "# The comparison\n", "if (abs(x-3.3) < eps):\n", " print(\"x == 3.3 to a precision of\",eps,\"is True\")\n", "else:\n", " print(\"x == 3.3 to a precision of\",eps,\"is False\")" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "5.0\n", "x == 3.3 is False\n" ] } ], "source": [ "import numpy as np\n", "\n", "print(np.sqrt(25))\n", "\n", "if (np.sqrt(36.+1.e-13) == np.round(np.sqrt(36))):\n", " print(\"x == 3.3 is True\")\n", "else:\n", " print(\"x == 3.3 is False\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 2: Subtracting two large numbers with a small difference\n", "\n", "Let us have $x = 1$ and $y = 1 + \\delta \\sqrt{2}$\n", "\n", "It follows that $$\\delta^{-1} (y - x) = \\sqrt{2} = 1.41421356237\\dots$$\n", "\n", "Let us test this relation on a computer for a very small value of $\\delta = 10^{-14}$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x = 1.0\n", "y = 1.0000000000000142\n", "(y-x) = 1.4210854715202004e-14\n", "(1/delta) * (y-x) = 1.4210854715202004\n", "The accurate value is sqrt(2) = 1.4142135623730951\n", "The difference is 0.006871909147105226\n" ] } ], "source": [ "from math import sqrt\n", "\n", "delta = 1.e-14\n", "x = 1.\n", "y = 1. + delta * sqrt(2)\n", "print(\"x = \", x)\n", "print(\"y = \", y)\n", "res = (1./delta)*(y-x)\n", "print(\"(y-x) = \",y-x)\n", "print(\"(1/delta) * (y-x) = \",res)\n", "print(\"The accurate value is sqrt(2) = \", sqrt(2))\n", "print(\"The difference is \", res - sqrt(2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Try smaller/bigger values of $\\delta$ and observe the behavior, e.g. $\\delta = 10^{-5}$ or $\\delta = 10^{-16}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 3: Roots of the quadratic equation\n", "\n", "The quadratic equation $$ax^2+bx+c = 0$$ has the following two roots\n", "$$\n", "x_{1,2} = \\frac{-b \\pm \\sqrt{b^2 - 4ac}}{2a}\n", "$$\n", "\n", "Let us calculate the roots for $a = 10^{-4}$, $b = 10^4$, and $c = 10^{-4}$" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x1 = -9.094947017729282e-09\n", "x2 = -100000000.0\n" ] } ], "source": [ "a = 1.e-4\n", "b = 1.e4\n", "c = 1.e-4\n", "\n", "x1 = (-b + sqrt(b*b - 4.*a*c)) / (2.*a)\n", "x2 = (-b - sqrt(b*b - 4.*a*c)) / (2.*a)\n", "\n", "print(\"x1 = \", x1)\n", "print(\"x2 = \", x2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do the results look accurate to you?\n", "\n", "The value of $x_1$ is not accurate due to subtracting two large numbers with small difference $b$ and $\\sqrt{b^2-4ac}$.\n", "\n", "Consider another form of the solution.\n", "By multiplying the numerator and denominator of the above expression for $x_{1,2}$ by $(-b\\mp\\sqrt{b^2-4ac})$ one obtains\n", "\n", "$$\n", "x_{1,2} = \\frac{2c}{-b \\mp \\sqrt{b^2-4ac}}\n", "$$\n", "\n", "Let us see what we get now" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x1 = -1e-08\n", "x2 = -109951162.7776\n" ] } ], "source": [ "x1 = 2*c / (-b - sqrt(b*b-4.*a*c))\n", "x2 = 2*c / (-b + sqrt(b*b-4.*a*c))\n", "\n", "print(\"x1 = \", x1)\n", "print(\"x2 = \", x2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This time $x_1$ is fine, but not $x_2$.\n", "\n", "One, therefore, has to combine the two forms to get accurate results for both $x_1$ and $x_2$.\n", "\n", "Consider writing a function which avoids large round-off errors for both $x_1$ and $x_2$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise:\n", "\n", "Consider now the quadratic equation with parameters $a = 10^{-4}$, $b = -10^4$, and $c = 10^{-4}$, where the sign of $b$ is now reversed. \n", "\n", "- Can you calculate both roots $x_{1,2}$ accurately in this case?\n", "- Which form should you use for $x_1$ and which for $x_2$ and why?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 4: Numerical derivative\n", "\n", "Consider a function \n", "$$f(x) = x(x-1)~.$$\n", "Its derivative is \n", "$$f'(x) = \\lim_{h \\to 0} \\frac{f(x+h)-f(x)}{h} = 2x - 1$$\n", "\n", "Let us calculate the derivative numerically by using small but finite values of $h$ ranging from $1$ down to $10^{-16}$ at $x = 1$\n", "and compare it to the correct result, $f'(1) = 1$." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "h f'(1) Relative error \n", "1 2.0 1.0 \n", "0.1 1.100000000000001 0.10000000000000098 \n", "0.01 1.010000000000001 0.010000000000000897\n", "0.001 1.0009999999998895 0.0009999999998895337\n", "0.0001 1.0000999999998899 9.999999988985486e-05\n", "1e-05 1.0000100000065513 1.0000006551269536e-05\n", "1e-06 1.0000009999177333 9.99917733279787e-07\n", "1e-07 1.0000001005838672 1.0058386723521551e-07\n", "1e-08 1.0000000039225287 3.922528746258536e-09\n", "1e-09 1.000000083740371 8.374037108183074e-08\n", "1e-10 1.000000082840371 8.284037100736441e-08\n", "1e-11 1.000000082750371 8.275037099991778e-08\n", "1e-12 1.0000889005833413 8.890058334132256e-05\n", "1e-13 0.9992007221627407 0.0007992778372593046\n", "1e-14 0.9992007221626509 0.0007992778373491216\n", "1e-15 1.1102230246251577 0.11022302462515765 \n", "1e-16 0.0 1.0 \n", "1e-17 0.0 1.0 \n", "1e-18 0.0 1.0 \n", "1e-19 0.0 1.0 \n" ] } ], "source": [ "def f(x):\n", " return x*(x-1.)\n", "\n", "def df_exact(x):\n", " return 2.*x - 1.\n", "\n", "def df_numeric(x,h):\n", " return (f(x+h) - f(x)) / h\n", "\n", "print(\"{:<10} {:<20} {:<20}\".format('h',\"f'(1)\",\"Relative error\"))\n", "x0 = 1.\n", "\n", "arr_h = []\n", "arr_df = []\n", "arr_err = []\n", "\n", "for i in range(0,-20,-1):\n", " h = 10**i\n", " df_val = df_numeric(x0,h)\n", " df_err = abs(df_numeric(x0,h) - df_exact(x0)) / df_exact(x0)\n", " print(\"{:<10} {:<20} {:<20}\".format(h,df_val,df_err))\n", " arr_h.append(h)\n", " arr_df.append(df_val)\n", " arr_err.append(df_err)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The accurace of our numerical derivative first increases as $h$ becomes smaller, as expected, but then increases again. This is due to large round-off error when $h$ become very small compared to $f$.\n", "\n", "Let us plot the dependence of the relative accuracy vs $h$" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAHTCAYAAADVvEfIAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVLJJREFUeJzt3Xtc1FX+P/DXh0EGb6CAIMg1Sw1JjRHvqGhhaKghWVles2SjTWL9lmZluiW/VdfLrkNlmaab5g3NdWuNCoQNSyVpNcrLBnIXcRUEFXA8vz9mmRgHmHGYYS68no/HPHDOHM55z8dh5j3nnM/5SEIIASIiIiJqloOlAyAiIiKydkyYiIiIiPRgwkRERESkBxMmIiIiIj2YMBERERHpwYSJiIiISA8mTERERER6MGEiIiIi0oMJExEREZEeTJjI7vzlL3+BJEkICQmxdCh2Kz8/H5MmTYKbmxskSUJCQkKzdVeuXIkDBw7olG/duhWSJOHEiRPmC9TOBQYGYs6cOWZrf+zYsRg7dqzRv2+O+PLz8yFJErZu3WrSdhvLzc3FW2+9hfz8fJ3H5syZg8DAQLP1TdbL0dIBEJnaRx99BAD46aef8P3332Po0KEWjsj+vPzyy/j+++/x0UcfoWfPnvD29m627sqVKxEbG4upU6e2XYDtxP79++Hi4mLpMNqUt7c3jh49it69e5utj9zcXCxfvhxjx47VSY7eeOMNLFy40Gx9k/ViwkR25cSJE/jxxx8xadIk/OMf/8DmzZutNmG6fv06OnXqZOkwjHL69GkMGTKESZCF3LhxAx07dsSDDz5o6VDajEqlwq1btyCXyzFs2DCLxWHORI2sG6fkyK5s3rwZAPD//t//w4gRI/Dpp5/i+vXrOvWKi4vx/PPPw8/PD05OTvDx8UFsbCwuXryoqXP16lX84Q9/wD333AO5XA5PT09MnDgRv/zyCwAgPT0dkiQhPT1dq+2mpgzmzJmDLl264NSpU4iMjETXrl0xfvx4AEBqaiqmTJkCX19fODs7495778WCBQtQUVGhE/cvv/yCp556Cl5eXpDL5fD398esWbNQW1uL/Px8ODo6IikpSef3MjIyIEkS9uzZ0+LxKygowDPPPANPT0/I5XLcf//9+POf/4zbt29rPefz58/jiy++gCRJkCSpyakLAJAkCTU1Nfj44481de+c4rl27Rp+97vfwcPDA+7u7oiJiUFJSYlOW7t27cLw4cPRuXNndOnSBRMmTMDJkydbfD7Ab1N/aWlpevuRJAlvvfWWTht3Ti01tPnNN9/gueeeg7u7O1xcXDBr1izU1NSgrKwM06dPR7du3eDt7Y1Fixahvr5eq826ujq8/fbb6NevH+RyOXr06IG5c+fi0qVLOn0/+uijSElJwYMPPghnZ2csX768ybgA/a9bAFi+fDmGDh0KNzc3uLi4IDQ0FJs3b4ax12Kvr6/HK6+8gp49e6JTp04YNWoUjh071mTdsrIyLFiwAL6+vnByckJQUBCWL1+OW7duaeo0/A2tWrUKb7/9NoKCgiCXy5GWlqbz93XgwAFIkoSvv/5ap693330XkiTh3//+NwD1F6onn3wSgYGB6NixIwIDA/HUU0/hwoULmt/ZunUrHn/8cQBARESE5nXb0N+dU3IPPvggwsPDdfpWqVTo1asXYmJiNGWG/p+TdeIIE9mNGzduYOfOnQgLC0NISAjmzZuH+fPnY8+ePZg9e7amXnFxMcLCwlBfX4/XXnsNAwYMwOXLl3H48GFcuXIFXl5euHbtGkaNGoX8/Hy8+uqrGDp0KKqrq5GRkYHS0lL069fvruOrq6vD5MmTsWDBAixevFjzAfGf//wHw4cPx/z58+Hq6or8/HysXbsWo0aNwqlTp9ChQwcAwI8//ohRo0bBw8MDK1aswH333YfS0lIcPHgQdXV1CAwMxOTJk/Hee+/hlVdegUwm0/S9ceNG+Pj44LHHHms2vkuXLmHEiBGoq6vDH//4RwQGBuLQoUNYtGgR/vOf/yA5ORmhoaE4evQoHnvsMfTu3Rtr1qwBgGan5I4ePYpx48YhIiICb7zxBgDoTCHNnz8fkyZNwo4dO1BYWIj/+7//wzPPPINvvvlGU2flypV4/fXXMXfuXLz++uuoq6vD6tWrER4ejmPHjiE4OFjv8Tekn7s1f/58xMTE4NNPP8XJkyfx2muv4datWzhz5gxiYmLw/PPP46uvvsKf/vQn+Pj4IDExEQBw+/ZtTJkyBZmZmXjllVcwYsQIXLhwAcuWLcPYsWNx4sQJdOzYUdPPDz/8gJ9//hmvv/46goKC0Llz5ybjMfR1m5+fjwULFsDf3x8A8N133+H3v/89iouL8eabb971cXjuueewbds2LFq0CA8//DBOnz6NmJgYXLt2TateWVkZhgwZAgcHB7z55pvo3bs3jh49irfffhv5+fnYsmWLVv2//OUv6NOnD9asWQMXFxfcd999On0/+uij8PT0xJYtWzRfQhps3boVoaGhGDBggOZ59+3bF08++STc3NxQWlqKd999F2FhYcjNzYWHhwcmTZqElStX4rXXXoNSqURoaCiA5keW5s6di4ULF+LcuXNa8X355ZcoKSnB3LlzAdz9/zlZIUFkJ7Zt2yYAiPfee08IIcS1a9dEly5dRHh4uFa9efPmiQ4dOojc3Nxm21qxYoUAIFJTU5utk5aWJgCItLQ0rfK8vDwBQGzZskVTNnv2bAFAfPTRRy0+h9u3b4v6+npx4cIFAUB89tlnmsfGjRsnunXrJsrLy/XGtH//fk1ZcXGxcHR0FMuXL2+x78WLFwsA4vvvv9cq/93vfickSRJnzpzRlAUEBIhJkya12F6Dzp07i9mzZ+uUb9myRQAQL7zwglb5qlWrBABRWloqhBCioKBAODo6it///vda9a5duyZ69uwppk+f3mL/hvYjhBAAxLJly3TaCAgI0HoODW3eGdPUqVMFALF27Vqt8kGDBonQ0FDN/Z07dwoAYt++fVr1jh8/LgCI5ORkrb5lMpnW8W8uLkNet3dSqVSivr5erFixQri7u4vbt29rHhszZowYM2ZMi7//888/CwDi5Zdf1ir/5JNPBACt+BYsWCC6dOkiLly4oFV3zZo1AoD46aefhBC//Q317t1b1NXVadVt6u8rMTFRdOzYUVy9elVTlpubKwCIv/71r83GfuvWLVFdXS06d+4sNmzYoCnfs2dPk3/bQqj/lgMCAjT3KyoqhJOTk3jttde06k2fPl14eXmJ+vp6IcTd/Z+TdeKUHNmNzZs3o2PHjnjyyScBAF26dMHjjz+OzMxMnDt3TlPviy++QEREBO6///5m2/riiy/Qp08fPPTQQyaNcdq0aTpl5eXliIuLg5+fHxwdHdGhQwcEBAQAAH7++WcA6vVOR44cwfTp09GjR49m2x87diwGDhwIpVKpKXvvvfcgSRKef/75FmP75ptvEBwcjCFDhmiVz5kzB0KIVo3EtGTy5Mla9xtGAxqmSQ4fPoxbt25h1qxZuHXrlubm7OyMMWPG6EyJGtuPMR599FGt+w2vqUmTJumUN+7n0KFD6NatG6Kjo7We06BBg9CzZ0+d5zRgwAD06dNHbzyGvm6/+eYbPPTQQ3B1dYVMJkOHDh3w5ptv4vLlyygvL9fbT2NpaWkAgKefflqrfPr06XB01J7EOHToECIiIuDj46P1vKOiogAAR44c0ao/efJkzQhrS+bNm4cbN25g165dmrItW7ZALpdjxowZmrLq6mq8+uqruPfee+Ho6AhHR0d06dIFNTU1mr+1u+Xu7o7o6Gh8/PHHmqnrK1eu4LPPPsOsWbM0x+Bu/8/J+jBhIrtw/vx5ZGRkYNKkSRBC4OrVq7h69SpiY2MB/HbmHKCeevL19W2xPUPq3K1OnTrpTEfdvn0bkZGRSElJwSuvvIKvv/4ax44dw3fffQdAPc0IqN+AVSqVQTG99NJL+Prrr3HmzBnU19fjgw8+QGxsLHr27Nni712+fLnJqTUfHx/N4+bg7u6udV8ulwP47bk3rCsLCwtDhw4dtG67du1qcq2XMf0Yw83NTeu+k5NTs+U3b97U3L948SKuXr0KJycnnedUVlam85xaOguxMUNet8eOHUNkZCQA4IMPPsC3336L48ePY+nSpQDu/ng0vC7ufH05OjrqHPOLFy/i73//u85z7t+/PwAY/bz79++PsLAwzZSeSqXC3/72N0yZMkXr/2LGjBnYuHEj5s+fj8OHD+PYsWM4fvw4evTo0arXwbx581BcXIzU1FQAwM6dO1FbW6u1vuxu/8/J+nANE9mFjz76CEII7N27F3v37tV5/OOPP8bbb78NmUyGHj16oKioqMX2DKnj7OwMAKitrdUqb+6NT5IknbLTp0/jxx9/xNatW7XWWZ0/f16rnpubG2Qymd6YAPWHwquvvgqlUolhw4ahrKwM8fHxen/P3d0dpaWlOuUNC6M9PDz0tmEODf3u3btXM/JmLnK5XOf/EzB9stiw8Pyf//xnk4937dpV635Tr52mGPK6/fTTT9GhQwccOnRI8xoG0OReWYZoSIrKysrQq1cvTfmtW7d0jpuHhwcGDBiAd955p8m2GpLzBoY+b0C9luiFF17Azz//jF9//RWlpaWa9UMAUFlZiUOHDmHZsmVYvHixpry2thb//e9/De6nKRMmTICPjw+2bNmCCRMmYMuWLRg6dKjW2rq7/T8n68OEiWyeSqXCxx9/jN69e+PDDz/UefzQoUP485//jC+++AKPPvoooqKisH37dpw5cwZ9+/Ztss2oqCi8+eab+OabbzBu3Lgm6zScKfPvf/8bEyZM0JQfPHjQ4NgbPhAaRjsavP/++1r3O3bsiDFjxmDPnj145513WkxenJ2d8fzzz2Pjxo3IysrCoEGDMHLkSL2xjB8/HklJSfjhhx80C10BYNu2bZAkCREREQY/r8bkcnmrvr1PmDABjo6O+M9//tPklKYpBQYGas6oavDNN9+gurrapP08+uij+PTTT6FSqUy67YUhr1tJkuDo6Kh1UsCNGzewfft2o/psOOvxk08+gUKh0JTv3r1b68w3QP28P//8c/Tu3Rvdu3c3qr/mPPXUU0hMTMTWrVvx66+/olevXpqRNED9vIUQOn9rH374IVQqlVbZ3Y4+ymQyzJw5E+vXr0dmZiZOnDih8zdsrv9zajtMmMjmffHFFygpKcGf/vSnJnclDgkJwcaNG7F582Y8+uijWLFiBb744guMHj0ar732Gh544AFcvXoV//znP5GYmIh+/fohISEBu3btwpQpU7B48WIMGTIEN27cwJEjR/Doo48iIiICPXv2xEMPPYSkpCR0794dAQEB+Prrr5GSkmJw7P369UPv3r2xePFiCCHg5uaGv//975qh/cYazpwbOnQoFi9ejHvvvRcXL17EwYMH8f7772t9Q33hhRewatUqZGdnN5lENuXll1/Gtm3bMGnSJKxYsQIBAQH4xz/+geTkZPzud78zaA1NUx544AGkp6fj73//O7y9vdG1a9dmE9WmBAYGYsWKFVi6dCl+/fVXPPLII+jevTsuXryIY8eOoXPnzprT7Ftr5syZeOONN/Dmm29izJgxyM3NxcaNG+Hq6mqS9hs8+eST+OSTTzBx4kQsXLgQQ4YMQYcOHVBUVIS0tDRMmTKlxTMam2PI63bSpElYu3YtZsyYgeeffx6XL1/GmjVrdBIJQ91///145plnsH79enTo0AEPPfQQTp8+rTmzrbEVK1YgNTUVI0aMwEsvvYS+ffvi5s2byM/Px+eff4733nvP6Knwbt264bHHHsPWrVtx9epVLFq0CA4Ov606cXFxwejRo7F69Wp4eHggMDAQR44cwebNm9GtWzetthquErBp0yZ07doVzs7OCAoK0plibGzevHn405/+hBkzZqBjx4544okntB431/85tSGLLjknMoGpU6cKJyenFs8ee/LJJ4Wjo6MoKysTQghRWFgo5s2bJ3r27Ck6dOggfHx8xPTp08XFixc1v3PlyhWxcOFC4e/vLzp06CA8PT3FpEmTxC+//KKpU1paKmJjY4Wbm5twdXUVzzzzjDhx4kSTZ8l17ty5ydhyc3PFww8/LLp27Sq6d+8uHn/8cVFQUNDkGVu5ubni8ccfF+7u7sLJyUn4+/uLOXPmiJs3b+q0O3bsWOHm5iauX79uyGEUQghx4cIFMWPGDOHu7i46dOgg+vbtK1avXi1UKpVWvbs5Sy4nJ0eMHDlSdOrUSQDQnHXVcKbZ8ePHteo3d/bhgQMHREREhHBxcRFyuVwEBASI2NhY8dVXX7XY/930U1tbK1555RXh5+cnOnbsKMaMGSNycnKaPUvuzjaXLVsmAIhLly5plTf1/19fXy/WrFkjBg4cKJydnUWXLl1Ev379xIIFC8S5c+c09Vo61nfGJYRhr9uPPvpI9O3bV8jlcnHPPfeIpKQksXnzZgFA5OXlaeoZcpZcw3H7wx/+IDw9PYWzs7MYNmyYOHr0aJPxXbp0Sbz00ksiKChIdOjQQbi5uQmFQiGWLl0qqqurhRC/nQm3evVqnb6aOkuuwZdffikACADi7NmzOo8XFRWJadOmie7du4uuXbuKRx55RJw+fbrJONevXy+CgoKETCbT6u/Os+QaGzFihAAgnn766SYfN/T/nKyTJISRO5URkdUqLy9HQEAAfv/732PVqlWWDoeIyOZxSo7IjhQVFeHXX3/F6tWr4eDgwGteERGZCLcVILIjH374IcaOHYuffvoJn3zyidZZS0REZDxOyRERERHpwREmIiIiIj2YMBERERHpwYSJiIiISA+eJWcit2/fRklJCbp27XpX2/kTERGR5QghcO3aNfj4+GhtdnonJkwmUlJSAj8/P0uHQUREREYoLCxscad5Jkwm0nBZisLCQp3LARAREZF1qqqqgp+fn94LIDNhMpGGaTgXFxcmTERERDZG33IaLvomIiIi0oMJUysplUoEBwcjLCzM0qEQERGRmXCnbxOpqqqCq6srKisrOSVHRERkIwz9/OYIExEREZEeTJiIiIiI9GDCRERERKQHEyYiIiIiPZgwEREREenBhKmRQ4cOoW/fvrjvvvvw4YcfWjocQKUC0tOBnTvVP1Uqtt+W7RMRkeVZy3u9ICGEEPX19eK+++4TRUVFoqqqStx7773i8uXLBv9+ZWWlACAqKytNE9C+fUL4+goB/Hbz9VWXs33zt09ERJbXBu/1hn5+c4Tpf44dO4b+/fujV69e6Nq1KyZOnIjDhw9bJpiUFCA2Figq0i4vLlaXp6SwfXO2T0RElmdl7/V2kzBlZGQgOjoaPj4+kCQJBw4c0KmTnJyMoKAgODs7Q6FQIDMzU/NYSUkJevXqpbnv6+uL4uLitghdm0oFLFyozqPv1FCWkGD8kCTbJyIia2eF7/V2kzDV1NRg4MCB2LhxY5OP79q1CwkJCVi6dClOnjyJ8PBwREVFoaCgAAAgmvhPaelCfLW1taiqqtK6mURmpm423ZgQQGGhuh7bN337RERkeVb4Xm83CVNUVBTefvttxMTENPn42rVr8eyzz2L+/Pm4//77sX79evj5+eHdd98FAPTq1UtrRKmoqAje3t7N9peUlARXV1fNzc/PzzRPpLTUtPXYPhER2RorfK+3m4SpJXV1dcjOzkZkZKRWeWRkJLKysgAAQ4YMwenTp1FcXIxr167h888/x4QJE5ptc8mSJaisrNTcCgsLTRNsC0maUfXYPhER2RorfK9vFwlTRUUFVCoVvLy8tMq9vLxQVlYGAHB0dMSf//xnRERE4MEHH8T//d//wd3dvdk25XI5XFxcsH37dgwbNgzjx483TbDh4YCvL9DcdKAkAX5+6nps3/TtExGR5Vnhe327SJga3LkmSQihVTZ58mScPXsW58+fx/PPP29Qm/Hx8cjNzcXx48dNE6RMBmzY0BCw9mMN99evV9dj+6Zvn4iILM8K3+vbRcLk4eEBmUymGU1qUF5erjPqdLeUSiWCg4MRFhbWqna0xMQAe/cCjc7aA6DOtvfuVT/O9s3XPhERWZ6VvddLoqnTw2ycJEnYv38/pk6dqikbOnQoFAoFkpOTNWXBwcGYMmUKkpKSWt1nVVUVXF1dUVlZCRcXl1a3B0B9umRmpnpRm7e3eujRlNk02yciImtn5vd6Qz+/HU3Wo4VVV1fj/Pnzmvt5eXnIycmBm5sb/P39kZiYiJkzZ2Lw4MEYPnw4Nm3ahIKCAsTFxVkwaj1kMmDsWLZvqfaJiMjyrOS93m4SphMnTiAiIkJzPzExEQAwe/ZsbN26FU888QQuX76MFStWoLS0FCEhIfj8888REBDQqn6VSiWUSiVU3CiRiIjIbtnllJwlmGVKjoiIiMzK0M/vdrHo25zMsuibiIiIrApHmEyEI0xERES2hyNMRERERCbChKmVOCVHRERk/zglZyKckiMiIrI9nJIjIiIiMhEmTERERER6MGFqJa5hIiIisn9cw2QiXMNERERke7iGiYiIiMhEmDARERER6cGEiYiIiEgPJkytxEXfRERE9o+Lvk2Ei76JiIhsDxd9ExEREZkIEyYiIiIiPZgwEREREenBhImIiIhIDyZMRERERHowYWolbitARERk/7itgIlwWwEiIiLbw20FiIiIiEyECRMRERGRHkyYiIiIiPRgwkRERESkh6OlAyCySioVkJkJlJYC3t5AeDggk1k6KiIishCOMDXy2GOPoXv37oiNjbV0KGRJKSlAYCAQEQHMmKH+GRioLiciIl0qFZCeDuzcqf6pUlk6IpNjwtTISy+9hG3btlk6DLKklBQgNhYoKtIuLy5WlzNpIiLS1k6+ZDJhaiQiIgJdu3a1dBhkKSoVsHAh0NTWZA1lCQl2+c2JiMgo7ehLps0kTBkZGYiOjoaPjw8kScKBAwd06iQnJyMoKAjOzs5QKBTIzMxs+0DJdmVm6v7RNyYEUFiorkdE1N61sy+ZNpMw1dTUYODAgdi4cWOTj+/atQsJCQlYunQpTp48ifDwcERFRaGgoEBTR6FQICQkROdWUlJy1/HU1taiqqpK60Y2rrTUtPWIiOxZO/uSaTNnyUVFRSEqKqrZx9euXYtnn30W8+fPBwCsX78ehw8fxrvvvoukpCQAQHZ2tsniSUpKwvLly03WHlkBb2/T1iMismft7EumzYwwtaSurg7Z2dmIjIzUKo+MjERWVpZZ+lyyZAkqKys1t8LCQrP0Q20oPBzw9QUkqenHJQnw81PXIyJq79rZl0y7SJgqKiqgUqng5eWlVe7l5YWysjKD25kwYQIef/xxfP755/D19cXx48ebrSuXy+Hi4oLt27dj2LBhGD9+vNHxk5WQyYANG9T/vjNpari/fj33YyIiAtrdl0y7SJgaSHf8pwkhdMpacvjwYVy6dAnXr19HUVERwsLC9P5OfHw8cnNzW0yuyIbExAB79wK9emmX+/qqy2NiLBMXEZG1aWdfMu0iYfLw8IBMJtMZTSovL9cZdTI1pVKJ4OBgg5IrshExMUB+PpCWBuzYof6Zl8dkiYjoTu3oS6YkRFPnA1o3SZKwf/9+TJ06VVM2dOhQKBQKJCcna8qCg4MxZcoUzaJvc6qqqoKrqysqKyvh4uJi9v6IiIishg1fTsrQz2+bOUuuuroa58+f19zPy8tDTk4O3Nzc4O/vj8TERMycORODBw/G8OHDsWnTJhQUFCAuLs6scSmVSiiVSqjsZJ8JIiKiuyaTAWPHWjoKs7KZEab09HRERETolM+ePRtbt24FoN64ctWqVSgtLUVISAjWrVuH0aNHt0l8HGEiIiKyPYZ+fttMwmTtmDARERHZHkM/v+1i0bclcdE3ERGR/eMIk4lwhImIiMj2cISJiIiIyESYMLUSp+SIiIjsH6fkTIRTckRERLaHU3JEREREJsKEiYiIiEgPJkytxDVMRERE9o9rmEyEa5iIiIhsD9cwEREREZkIEyYiIiIiPZgwEREREenBhKmVuOibiIjI/nHRt4lw0TcREZHt4aJvIiIiIhNhwkRERESkBxMmIiIiIj2YMBERERHpwYSJiIiISA8mTK3EbQWIiIjsH7cVMBFuK0BERGR7uK0AERERkYkwYSIiIiLSgwkTERERkR5MmIiIiIj0YMJEREREpAcTpv8pLCzE2LFjERwcjAEDBmDPnj2WDomIiMg0VCogPR3YuVP9U6WydEQ2x9HSAVgLR0dHrF+/HoMGDUJ5eTlCQ0MxceJEdO7c2dKhERERGS8lBVi4ECgq+q3M1xfYsAGIibFcXDaGI0z/4+3tjUGDBgEAPD094ebmhv/+97+WDYqIiKg1UlKA2FjtZAkAiovV5SkplonLBtlMwpSRkYHo6Gj4+PhAkiQcOHBAp05ycjKCgoLg7OwMhUKBzMxMo/o6ceIEbt++DT8/v1ZGTUREZCEqlXpkqan9qRvKEhI4PWcgm0mYampqMHDgQGzcuLHJx3ft2oWEhAQsXboUJ0+eRHh4OKKiolBQUKCpo1AoEBISonMrKSnR1Ll8+TJmzZqFTZs2tRhPbW0tqqqqtG5ERERWIzNTd2SpMSGAwkJ1PdLLJi+NIkkS9u/fj6lTp2rKhg4ditDQULz77ruasvvvvx9Tp05FUlKSQe3W1tbi4YcfxnPPPYeZM2e2WPett97C8uXLdcp5aRQiIrIKO3cCM2bor7djB/DUU+aPx0q1q0uj1NXVITs7G5GRkVrlkZGRyMrKMqgNIQTmzJmDcePG6U2WAGDJkiWorKzU3AoLC42KnYiIyCy8vU1br52zi4SpoqICKpUKXl5eWuVeXl4oKyszqI1vv/0Wu3btwoEDBzBo0CAMGjQIp06dara+XC6Hi4sLtm/fjmHDhmH8+PGteg5EREQmFR6uPhtOkpp+XJIAPz91PdLLrrYVkO54UQghdMqaM2rUKNy+ffuu+4yPj0d8fLxmSI+IiMgqyGTqrQNiY9XJUeMVOA2fjevXq+uRXnYxwuTh4QGZTKYzmlReXq4z6mRqSqUSwcHBCAsLM2s/REREdy0mBti7F+jVS7vc11ddzn2YDGYXCZOTkxMUCgVSU1O1ylNTUzFixAiz9h0fH4/c3FwcP37crP0QEREZJSYGyM8H0tLUC7zT0oC8PCZLd8lmpuSqq6tx/vx5zf28vDzk5OTAzc0N/v7+SExMxMyZMzF48GAMHz4cmzZtQkFBAeLi4swal1KphFKphIr7WBARkbWSyYCxYy0dhU2zmW0F0tPTERERoVM+e/ZsbN26FYB648pVq1ahtLQUISEhWLduHUaPHt0m8Rl6WiIRERFZD0M/v41OmLZt2wYAmDBhgtnXCdkCJkxERES2x9DPb6On5ObMmQNHR0dcvXrV2CbsAqfkiIiI7J/RI0weHh4A1HsgEUeYiIiIbJHZd/ru168fKisrUV1dbWwTRERERDbB6IRpzpw5UKlU+PDDD00Zj83hPkxERET2r1VnyU2fPh2fffYZ1q1bh+effx6OjjazS4HJcUqOiIjI9pj9LLl58+ZBCIF9+/ahpqYG3bt3R1hYGDw9PSFrZpt1SZKwefNmY7qzekyYiIiIbI/ZEyYHBwdIkgRDfr2hniRJdns2GRMmIiIi22P2bQVmzZpl8IVt7Rm3FSAiIrJ/NrPTt7XjCBMREZHtMfu2AkRERETtBRMmIiIiIj1Msg/AxYsXsXfvXpw4cQLl5eWQJAk9evRAWFgYpk2bxmvNERERkU1r1RomlUqFN954A2vXrkV9fT0AaM6aa1gQ3qFDB/zhD3/AihUrmt1uwJY1XvR99uxZrmEiIiKyIWbfVgAAnn76aXz66acQQkAul2Pw4MHw9fUFABQVFeHEiROora2FJEmYMWMGtm/fbmxXVo+LvomIiGyP2Rd9HzhwADt37oQQAomJiSgtLUVmZiZ27tyJnTt3IjMzE2VlZVi0aBGEENixYwcOHjxobHdEREREFmN0wrR582ZIkoSlS5dizZo16Natm04dV1dXrFq1CkuXLoUQAh988EFrYiUiIiKyCKOn5Hr27InLly+joqICrq6uLdatrKyEu7s7PDw8UFZWZlSg1o5TckRERLbH7FNyV65cgaurq95kCYCm3pUrV4ztjoiIiMhijE6YunfvjsrKSlRVVemtW1lZicrKSnTv3t3Y7oiIiIgsxuiEKSwsDLdv38a6dev01l23bh1u376NwYMHG9ud1VIqlQgODkZYWJilQyEiIiIzMTphmjt3LoQQ+OMf/4g33ngD1dXVOnWuXbuG119/HX/84x8hSRLmz5/fqmCtUXx8PHJzc3H8+HFLh0JERERm0qp9mJ588kns3r0bkiTB2dkZYWFh6NWrFyRJQmFhIU6cOIGbN29CCIEnnngCO3fuNGXsVoWLvomIiGyPoZ/frbo0yvbt2+Hr64u//OUvuHHjBjIyMjQ7fDfkYY6Ojli4cCFWrlzZmq6IiIiILKZVI0wNSkpKsG/fPs215ADA09MTgwcPxrRp0+Dj49PqQK0dR5iIiIhsj9kvjZKRkQEAGDBgQJObVrY3TJiIiIhsj9mn5MaOHQsHBwdcunTJ2CaIiIiIbILRCZOrqytkMpnd7K107do1jBs3DvX19VCpVHjppZfw3HPPWTosIiJqD1QqIDMTKC0FvL2B8HBAJrN0VNSI0QnTvffei3//+9+ora2FXC43ZUwW0alTJxw5cgSdOnXC9evXERISgpiYGLi7u1s6NCIismcpKcDChUBR0W9lvr7Ahg1ATIzl4iItRu/D9OSTT6K+vh67d+82ZTwWI5PJ0KlTJwDAzZs3oVKpYIL18ERERM1LSQFiY7WTJQAoLlaXp6RYJi7SYXTCtHDhQgwfPhwvvvgiPv/8c1PG1KSMjAxER0fDx8cHkiThwIEDOnWSk5MRFBQEZ2dnKBQKZGZm3lUfV69excCBA+Hr64tXXnkFHh4eJoqeiIjoDiqVemSpqS/nDWUJCep6ZHFGT8mtXLkSo0ePxqlTpxAdHY3+/ftj5MiR8PT0hKyFedc333zTqP5qamowcOBAzJ07F9OmTdN5fNeuXUhISEBycjJGjhyJ999/H1FRUcjNzYW/vz8AQKFQoLa2Vud3v/zyS/j4+KBbt2748ccfcfHiRcTExCA2NhZeXl5NxlNbW6vVliHX1CPSMPd6Ba6HILJ+mZm6I0uNCQEUFqrrjR3bZmFRM4SRJEkSDg4OQpIkzc3BwUHvzRQAiP3792uVDRkyRMTFxWmV9evXTyxevNioPuLi4sTu3bubfXzZsmUCgM6tsrLSqP6oHdm3TwhfXyHUb4fqm6+vutwW2ici09ixQ/vvtLnbjh2WjtSuVVZWGvT5bfQI0+jRozW7eltaXV0dsrOzsXjxYq3yyMhIZGVlGdTGxYsX0bFjR7i4uKCqqgoZGRn43e9+12z9JUuWIDExUXO/qqoKfn5+xj0Baj8a1ivcOQTfsF5h797WLfI0d/tEZDre3qatR2ZldMKUnp5uwjBap6KiAiqVSmf6zMvLC2VlZQa1UVRUhGeffRZCCAgh8OKLL2LAgAHN1pfL5ZDL5VAqlVAqlVBxjpn00bdeQZLU6xWmTDFu+szc7RORaYWHq8+GKy5u+u9WktSPh4e3fWykw+iEadu2bQCACRMmNLvOp63dOeIlhDB4FEyhUCAnJ+eu+4yPj0d8fLxmp1CiZpl7vQLXQxDZFplMvXVAbKw6OWqcNDV8dq1fzy84VsLos+TmzJmDZ599Fl27djVlPEbx8PCATCbTGU0qLy83ezKnVCoRHByMsLAws/ZDdqC01LT12rp9IjK9mBj1VHmvXtrlvr6cQrcyRidMbm5ucHV11exdZElOTk5QKBRITU3VKk9NTcWIESPM2nd8fDxyc3Nx/Phxs/ZDdsDc6xW4HoLINsXEAPn5QFoasGOH+mdeHpMlK2P0lFy/fv3w/fffo7q6Gl26dDFlTE2qrq7G+fPnNffz8vKQk5MDNzc3+Pv7IzExETNnzsTgwYMxfPhwbNq0CQUFBYiLizNrXFzDRAYz93oFrocgsl0yGafKrZ2xp+F98MEHQpIksW7dOmObuCtpaWlNnsY/e/ZsTR2lUikCAgKEk5OTCA0NFUeOHGmT2IQw/LREauf27RNCktS3xqcNN5S19tR/c7dPRGRnDP38loQw/vof06dPx2effYZ169bh+eefh6Oj0QNWNq9h0XdlZSVcXFwsHQ5Zs6auG+Xnp17caYoheHO3T0RkRwz9/DY6YZo3bx6EENi3bx9qamrQvXt3hIWFtbjTtyRJ2Lx5szHdWa3GU3Jnz55lwkSG4U7fRERWwewJk4ODAyRJMugCtQ31JEmy27U+HGEiIiKyPYZ+fhs9hzZr1iyr2embiIiIyJyMTpi2bt1qwjBsF8+SIyIisn+tWvRNv+GUHBERke0x9PPb6I0riYiIiNoLk+wDcPDgQRw+fBgXLlzAjRs38PXXX2seq6mpwY8//ghJkjB8+HBTdEdERETUplqVMBUWFiImJgY//PADgKYvdiuXy/HUU0+hqKgIOTk5eOCBB1rTpdXhGiYiIiL7Z/SU3PXr1xEZGYns7Gz06tUL8fHx6Ny5s049R0dHzJ8/H0IIfPbZZ60K1hrxWnJERET2z+iESalU4syZMwgNDcXPP/+Mv/zlL81eU27KlCkAgC+//NLY7oiIiIgsxuiEae/evZAkCWvXrm1yZKmxkJAQODo64uzZs8Z2R0RERGQxRidMZ86cgUwmw8iRI/V34uAAV1dXXLlyxdjuiIiIiCzG6ISptrYWHTt2bPa6cXeqqamBXC43tjurpVQqERwcjLCwMEuHQkRERGZidMLk6emJ6upqXL16VW/dH3/8ETdv3oSvr6+x3VktLvomIiKyf0YnTCNGjAAA7N69W2/dd955B5IkYcyYMcZ2R0RERGQxRidMcXFxEELgrbfeQm5ubpN1rl+/jvj4eOzdu1fzO0RERES2xuiNK8eMGYNnn30WmzdvxtChQzFp0iTU1NQAAFavXo1Tp07hH//4h2bKLiEhAQMHDjRJ0ERERERtqVUX31WpVPjDH/6Av/71r2hopvFO3w07f7/88stYvXq1zi7g9oQX3yUiIrI9hn5+typhavDTTz/hww8/xLfffouSkhKoVCr07NkTI0eOxHPPPdcuRpaYMBEREdmeNk2Y2rPG15I7e/YsEyYiIiIbwoSpjXGEiYiIyPYY+vlt9FlyRERERO0FEyYiIiIiPZgwEREREenBhImIiIhIDyZMRERERHowYbrD9evXERAQgEWLFlk6FCIishYqFZCeDuzcqf6pUlk6ImpjRl8axV698847GDp0qKXDICIia5GSAixcCBQV/Vbm6wts2ADExFguLmpTHGFq5Ny5c/jll18wceJES4dCZL34TZvak5QUIDZWO1kCgOJidXlKimXiojZnsoTp0qVLOHHiBDIyMkzVpJaMjAxER0fDx8cHkiThwIEDOnWSk5MRFBQEZ2dnKBQKZGZm3lUfixYtQlJSkokiJrJDKSlAYCAQEQHMmKH+GRjIDw2yTyqVemSpqf2dG8oSEviloZ1odcJ08OBBhIaGomfPnhg6dCjGjRun9fiVK1fwyCOP4JFHHkFNTY3R/dTU1GDgwIHYuHFjk4/v2rULCQkJWLp0KU6ePInw8HBERUWhoKBAU0ehUCAkJETnVlJSgs8++wx9+vRBnz59DIqntrYWVVVVWjciu8Zv2tTeZGbqvt4bEwIoLFTXI/snWiEpKUk4ODgISZI0NwcHB516jz32mHBwcBC7d+9uTXcaAMT+/fu1yoYMGSLi4uK0yvr16ycWL15sUJuLFy8Wvr6+IiAgQLi7uwsXFxexfPnyZusvW7ZMANC5VVZW3vXzIbJ6t24J4esrhPojQvcmSUL4+anrEdmLHTuaf803vu3YYelIqRUqKysN+vw2eoTp+++/x9KlS+Ho6Ih169ahoqICXl5eTdZ95plnIITAwYMHje2uRXV1dcjOzkZkZKRWeWRkJLKysgxqIykpCYWFhcjPz8eaNWvw3HPP4c0332y2/pIlS1BZWam5FRYWtuo5EFk1ftOm9sjb27T1yKYZfZbchg0bAKgTh4ULF7ZYd8yYMQCA48ePG9tdiyoqKqBSqXQSNi8vL5SVlZmlT7lcDrlcDqVSCaVSCRXnsMmelZaath6RLQgPV58NV1zc9DomSVI/Hh7e9rFRmzM6YfrXv/4FAHjxxRf11nV3d0eXLl1QXFxsbHcGkSRJ674QQqfMEHPmzDG4bnx8POLj4zVXOyayS/ymTe2RTKbeOiA2Vp0cNU6aGj5b1q9X1yO7Z/SUXHl5Obp27QoPDw+D6nfo0AF1dXXGdtciDw8PyGQyndGk8vLyZqcJTUWpVCI4OBhhYWFm7YfIohq+aTf3BUSSAD8/ftMm+xMTA+zdC/TqpV3u66su5z5M7YbRCVOnTp1w/fp13L59W2/dqqoqXL16Fd27dze2uxY5OTlBoVAgNTVVqzw1NRUjRowwS58N4uPjkZuba7bpRiKr0PBNG9BNmvhNm+xdTAyQnw+kpQE7dqh/5uUxWWpnjE6Y+vTpA5VKhX//+9966+7btw9CCAwcONDY7lBdXY2cnBzk5OQAAPLy8pCTk6PZNiAxMREffvghPvroI/z88894+eWXUVBQgLi4OKP7NARHmKjd4Ddtas9kMmDsWOCpp9Q/+eWg/TH2NLy3335bSJIknnjiCU1Zz549dbYVOHfunPD09BQODg4iOTnZ2O5EWlpak6fxz549W1NHqVSKgIAA4eTkJEJDQ8WRI0eM7u9uGXpaIpHNu3VLiLQ09anUaWncSoCIbJqhn9+SEE0t/devqqoKffr0waVLlzBr1iy8+uqriIiIQHl5OVQqFX799Vd8+umnWL16NSorKxEUFITc3FzI5XKTJXvWpGHRd2VlJVxcXCwdDhERERnA0M9voxMmQL0X0yOPPKKzy3XHjh1x48YNAOoz1dzd3fHVV1+1akrOWjXeVuDs2bNMmIiIiGxImyRMAHDhwgUsWrQI+/fv11kALkkSpkyZgj//+c8ICgpqTTdWjyNMREREtsfQz2+j92FqEBAQgD179uDKlSs4evQoSkpKoFKp0LNnT4wYMQI9evRobRdEREREFtXqhKlB9+7dMXHiRFM1ZzO40zcREZH9M3pKLisry+x7HNkSTskRERHZHkM/v43eh2nUqFG47777sGLFCvz666/GNkNERERk9YweYXJwUOdaDddqGz58OGbPno3p06e3y2uqcYSJiIjI9ph9hOn8+fNYtmwZevfuDSEEsrKyEBcXB29vbzzxxBP4+9//3i7W9XCnbyIiIvvX6m0FAOC7777Dtm3bsHv3bvz3v/9VNyxJ8PDwwFNPPYVnnnkGgwcPbnWw1owjTERERLanzfZhaqy+vh6ff/45tm3bhs8//xy1tbWaKbu+ffti9uzZePXVV03VnVVhwkRERGR7LJIwNXb16lXs2rULf/vb3/Dtt9+qO5Mku52mY8JERERke9ps48rmdOvWDfPmzUOPHj1w8+ZNZGdnm6srIiIiIrMyS8J09OhRbN++Hbt378aVK1c05fa46zc3riQiIrJ/JpuS+/XXX/G3v/0N27dv1+zLJISAXC5HdHQ0Zs2ahaioKMhkMlN0Z3U4JUdERGR72mRK7urVq9i9eze2bduGo0ePAlAnSQAwcuRIzJo1q93uy0RERET2w+iEKTY2Fv/4xz9QV1enSZJ69+6NmTNnYubMmQgKCjJZkERERESWZHTClJKSAkB90d3p06dj5syZvLYcERER2SWjE6aGdUnR0dFwcnIyZUxEREREVsXohOmzzz4zZRxEREREVsvoa8mRGq8lR0REZP/MttN3e8NtBYiIiGyPoZ/fBo0wyWQyyGQy9O/fX6fsbm6OjmbbWJyIiIjIbAzKYBoGoRoPRnFgioiIiNoLgxKmtLQ0AECnTp10yoiIiIjsHdcwmQjXMBEREdkek65hIiIiImrPjE6Yxo0bh8cff9zg+k899RTGjx9vbHdtwtHREYMGDcKgQYMwf/58S4dDRER3Q6UC0tOBnTvVP1UqS0dEdsTo09bS09PRs2dPg+t/9913KCgoMLa7NtGtWzfk5ORYOgwiIrpbKSnAwoVAUdFvZb6+wIYNQEyM5eIiu9FmU3K3b9+GJElt1R0REbUXKSlAbKx2sgQAxcXq8v9d+5SoNdokYVKpVCgvL0fnzp2NbiMjIwPR0dHw8fGBJEk4cOCATp3k5GQEBQXB2dkZCoUCmZmZd9VHVVUVFAoFRo0ahSNHjhgdKxERtRGVSj2y1NT5Sw1lCQmcnqNWM3hKrqqqClevXtUqU6lUKCwsbHZPJiEErl69ii1btqC2thYDBgwwOtCamhoMHDgQc+fOxbRp03Qe37VrFxISEpCcnIyRI0fi/fffR1RUFHJzc+Hv7w8AUCgUqK2t1fndL7/8Ej4+PsjPz4ePjw9Onz6NSZMm4dSpU82umK+trdVqq6qqyujnRkRERsrM1B1ZakwIoLBQXW/s2DYLi+yPwQnTunXrsGLFCq2yiooKBAYGGvT7kiRh5syZdxVcY1FRUYiKimr28bVr1+LZZ5/VLNZev349Dh8+jHfffRdJSUkAgOzs7Bb78PHxAQCEhIQgODgYZ8+exeDBg5usm5SUhOXLlxvzVIiIyFRKS01bj6gZdzUlJ4TQ3CRJ0rrf0s3HxwcrVqzAiy++aJYnUVdXh+zsbERGRmqVR0ZGIisry6A2rly5ohkxKioqQm5uLu65555m6y9ZsgSVlZWaW2FhofFPgIiIjOPtbdp6RM0weIQpISEBc+bMAaBOnO655x706NEDx44da/Z3HBwc4OLiAldX11YH2pKKigqoVCp4eXlplXt5eaGsrMygNn7++WcsWLAADg4OkCQJGzZsgJubW7P15XI55HI5lEollEolVJwfJyJqe+Hh6rPhioubXsckSerHw8PbPjayKwYnTK6urlqJz+jRo+Hh4YGAgACzBGaMO8/CaxgJM8SIESNw6tSpu+4zPj4e8fHxmp1CiYioDclk6q0DYmPVyVHjpKnh/X/9enU9olYw+iy59PR07N2715SxGM3DwwMymUxnNKm8vFxn1MnUlEolgoODERYWZtZ+iIioGTExwN69QK9e2uW+vupy7sNEJmAXl0ZxcnKCQqFAamqqVnlqaipGjBhh1r7j4+ORm5uL48ePm7UfIiJqQUwMkJ8PpKUBO3aof+blMVkikzF6p+/G6urqkJOTg6KiItTU1DS7zQAAzJo1y6g+qqurcf78ec39vLw85OTkwM3NDf7+/khMTMTMmTMxePBgDB8+HJs2bUJBQQHi4uKM6s9QXMNERGQlZDJuHUBmI4mWshs9amtrsXTpUmzatAk1NTX6O5Mk3Lp1y6i+0tPTERERoVM+e/ZsbN26FYB648pVq1ahtLQUISEhWLduHUaPHm1Uf3fL0KsdExERkfUw9PPb6ITp1q1beOihh5CZmQkhBDw9PVFeXg4HBwf4+PigoqICN2/eBAB06dIF7u7uANQjQ/aICRMREZHtMfTz2+g1TJs3b0ZGRgZ8fHxw4sQJzYJrT09PFBQUoLq6GmlpaRgxYgRu3bqFt99+2y6TJS76JiIisn9GjzCNHTsWmZmZ2LJli2ZdkoODA3r27ImSkhJNPZVKhaioKBw5cgRHjx5FaGioaSK3MhxhIiIisj1mH2E6ffo0ACA2Nlar/M7FzzKZDGvXrkV9fT3WrFljbHdEREREFmN0wnTt2jW4uLigU6dOmjK5XI7q6mqduiEhIejatSsyMzON7c5qcUqOiIjI/hk9Jefn54fy8nLN9dcAoFevXigrK0NpaSk8PT015UIIdOrUCUIIzUJwe8MpOSIiIttj9ik5X19f3Lp1S2t37ZCQEADAP//5T6266enpqK2t5aVDiIiIyCYZnTCNGTMGAPDtt99qyh577DEIIbBo0SLs2bMH586dw969ezF79mxIkoRx48a1PmIiIiKiNmb0lFxWVhZGjRqFqVOnIiUlBQBQX18PhUKB06dPa130VgiBLl264NixY+jXr59pIrcSjXf6Pnv2LKfkiIiIbIjZN64UQqCgoACOjo7o1eiCh5cuXUJCQgL279+PmzdvQpIkjBw5EuvXr7fbLQUArmEiIiKyRWZPmPS5desWLl26BBcXF3Tu3NkcXVgVJkxERES2x9DPb5NcfLfJhh0d4e3tba7miYiIiNqM0Yu+iYiIiNoLJkytxI0riYiI7J9Ba5hkMplpOpMk3Lp1yyRtWRuuYSIiIrI9Jl3DZKZ14UREREQ2waCEKS0tzdxxEBEREVktgxKmhl29iYiIiNojLvomIiIi0oMJExEREZEerU6YioqKkJiYiP79+6NLly5wdNSe5bty5QpWrlyJpKQk3L59u7XdWR1uK0BERGT/WnVplK+++gqPP/44qqqqNGfSSZIElUqlVW/o0KE4ceIEvvrqK0RERLQuYivFbQWIiIhsj6Gf30aPMBUWFiI2NhaVlZWIjo7G3r170b179ybrzps3D0IIHDhwwNjuiIiIiCzG6IRp7dq1qKqqwvTp03HgwAHExMTAycmpyboTJkwAAPzrX/8ytjsiIiIiizE6YTp8+DAkScIf//hHvXUDAwMhl8uRl5dnbHdEREREFmN0wnThwgV07NgR9913n0H1u3TpgurqamO7IyIiIrIYoxMmBwcHg896q6+vR2VlJbp27Wpsd0REREQWY3TC5Ofnh9raWhQWFuqtm5aWhlu3buHee+81trs2kZeXh4iICAQHB+OBBx5ATU2NpUMiIrIfKhWQng7s3Kn+eccZ1UTWzOiEafz48QCA999/v8V6N27cwOLFiyFJEh555BFju2sTc+bMwYoVK5Cbm4sjR45ALpdbOiQiIvuQkgIEBgIREcCMGeqfgYHqciIbYHTC9PLLL8PR0RFr1qzBxx9/3GSdrKwsjB49Gjk5OejUqRNeeOEFowM1t59++gkdOnRAeHg4AMDNzU1nE04iIjJCSgoQGwsUFWmXFxery5k0kQ0wOmG65557sHHjRtTV1WHevHnw8fHBlStXAACTJk2Cv78/wsPDkZ2dDUmSsGnTJnh5eRkdaEZGBqKjo+Hj4wNJkprc0yk5ORlBQUFwdnaGQqFAZmamwe2fO3cOXbp0weTJkxEaGoqVK1caHSsREf2PSgUsXAg0tUdyQ1lCAqfnyOq1agjlueeeQ48ePfDiiy+ipKREU/7FF19o/u3t7Y3k5GRMmTKlNV2hpqYGAwcOxNy5czFt2jSdx3ft2oWEhAQkJydj5MiReP/99xEVFYXc3Fz4+/sDABQKBWpra3V+98svv0R9fT0yMzORk5MDT09PPPLIIwgLC8PDDz/cZDy1tbVabVVVVbXq+RER2aXMTN2RpcaEAAoL1fXGjm2zsIjuVqvnnKZOnYpJkybhiy++wL/+9S+UlJRApVKhZ8+eGDlyJKKjo02yFigqKgpRUVHNPr527Vo8++yzmD9/PgBg/fr1OHz4MN59910kJSUBALKzs5v9fV9fX4SFhcHPzw8AMHHiROTk5DSbMCUlJWH58uXGPh0iovahtNS09YgsxOiEKSMjAwAwYMAAdOvWDZMnT8bkyZNNFtjdqKurQ3Z2NhYvXqxVHhkZiaysLIPaCAsLw8WLF3HlyhW4uroiIyMDCxYsaLb+kiVLkJiYqLlfVVWlSbaIiOh/vL1NW4/IQoxOmMaOHQsHBwdcunTJlPEYpaKiAiqVSmeNlJeXF8rKygxqw9HREStXrsTo0aMhhEBkZCQeffTRZuvL5XLI5XIolUoolUqdCw4TERGA8HDA11e9wLupdUySpH78fyfcEFkroxMmV1dXyGSyZi+4awmSJGndF0LolLVE37RfU+Lj4xEfH6+52jERETUikwEbNqjPhpMk7aSp4f15/Xp1PSIrZvRZcvfeey+uXbvW5CLqtubh4QGZTKYzmlReXt6qM/MMoVQqERwcjLCwMLP2Q0Rks2JigL17gV69tMt9fdXlMTGWiYvoLhidMD355JOor6/H7t27TRmPUZycnKBQKJCamqpVnpqaihEjRpi17/j4eOTm5uL48eNm7YeIyKbFxAD5+UBaGrBjh/pnXh6TJbIZRk/JLVy4ECkpKXjxxRfh7u6OiRMnmjIuHdXV1Th//rzmfl5eHnJycuDm5gZ/f38kJiZi5syZGDx4MIYPH45NmzahoKAAcXFxZo2La5iIiAwkk3HrALJZkhBNrcLTb8WKFbhx4waUSiVqamrQv39/jBw5Ep6enpC1MBf95ptvGhVoeno6IiIidMpnz56NrVu3AlBvXLlq1SqUlpYiJCQE69atw+jRo43q7241rGGqrKyEi4tLm/RJRERErWPo57fRCZODgwMkSULjXzdkgbW9jsQwYSIiIrI9hn5+Gz0lN3r06Ls6A81ecUqOiIjI/hk9wkTaOMJERERkewz9/Db6LDkiIiKi9oIJUytxHyYiIiL7xyk5E+GUHBERke3hlBwRERGRiTBhIiIiItKDCVMrcQ0TERGR/eMaJhPhGiYiIiLbwzVMRERERCbChImIiIhIDyZMRERERHowYWolLvomIiKyf1z0bSJc9E1ERGR7uOibiIiIyESYMBERERHpwYSJiIiISA8mTERERER6MGEiIiIi0oMJUytxWwEiIiL7x20FTITbChAREdkebitAREREZCJMmIiIiIj0YMJEREREpAcTJiIiIiI9mDARERER6cGE6X/OnDmDQYMGaW4dO3bEgQMHLB0WEVHbUamA9HRg5071T5XK0hERWQ1HSwdgLfr27YucnBwAQHV1NQIDA/Hwww9bNigioraSkgIsXAgUFf1W5usLbNgAxMRYLi4iK8ERpiYcPHgQ48ePR+fOnS0dChGR+aWkALGx2skSABQXq8tTUiwTF5EVsZmEKSMjA9HR0fDx8YEkSU1OlyUnJyMoKAjOzs5QKBTIzMw0qq/du3fjiSeeaGXERGQUc08L2Xr7pqZSqUeWmtrDuKEsIcH6nweRmdlMwlRTU4OBAwdi48aNTT6+a9cuJCQkYOnSpTh58iTCw8MRFRWFgoICTR2FQoGQkBCdW0lJiaZOVVUVvv32W0ycOLHFeGpra1FVVaV1I6JWSkkBAgOBiAhgxgz1z8BA041w2Hr75pCZqTuy1JgQQGGhuh5ReyZsEACxf/9+rbIhQ4aIuLg4rbJ+/fqJxYsX31Xb27ZtE08//bTeesuWLRMAdG6VlZV31R8R/c++fUJIkhDqj+jfbpKkvu3b177bN5cdO3Rjbuq2Y4elIyUyi8rKSoM+v21mhKkldXV1yM7ORmRkpFZ5ZGQksrKy7qotQ6fjlixZgsrKSs2tsLDwrvohokbMPS1k6+2bk7e3aesR2Sm7SJgqKiqgUqng5eWlVe7l5YWysjKD26msrMSxY8cwYcIEvXXlcjlcXFywfft2DBs2DOPHj7/ruInof8w9LWTr7ZtTeLj6bDhJavpxSQL8/NT1iNoxu0iYGkh3/MELIXTKWuLq6oqLFy/CycnJ4N+Jj49Hbm4ujh8/bvDvENEdSktNW8/e2jcnmUy9dQCgmzQ13F+/Xl2PqB2zi4TJw8MDMplMZzSpvLxcZ9TJ1JRKJYKDgxEWFmbWfojsmrmnhWy9fXOLiQH27gV69dIu9/VVl3MfJiJIQjQ16W7dJEnC/v37MXXqVE3Z0KFDoVAokJycrCkLDg7GlClTkJSUZPaYqqqq4OrqisrKSri4uJi9PyK7olKpzyYrLm56HZAkqT+88/KMG+mw9fbbikqlnjYsLVUnd+Hh1h0vkQkY+vltMzt9V1dX4/z585r7eXl5yMnJgZubG/z9/ZGYmIiZM2di8ODBGD58ODZt2oSCggLExcWZNS6lUgmlUgmVNS7mJLIVDdNCsbHq5KJx0mGKaSFbb7+tyGTA2LGWjoLIOrXFKXumkJaW1uRp/LNnz9bUUSqVIiAgQDg5OYnQ0FBx5MiRNovP0NMSiagF+/YJ4eurfTq7n5/pTsm39faJyOQM/fy2ySk5a8QpOSITMfe0kK23T0QmZejnNxOmVmo8JXf27FkmTERERDaECVMb4wgTERGR7TH089suthUgIiIiMicmTK3EfZiIiIjsH6fkTIRTckRERLaHU3JEREREJsKEiYiIiEgPJkytxDVMRERE9o9rmEyEa5iIiIhsD9cwEREREZkIEyYiIiIiPZgwEREREenBhKmVuOibiIjI/nHRt4lw0TcREZHt4aJvIiIiIhNhwkRERESkBxMmIiIiIj2YMBERERHpwYSJiIiISA8mTK3EbQWIiIjsH7cVMBFuK0BERGR7uK0AERERkYkwYSIiIiLSgwkTERERkR5MmIiIiIj0YMJEREREpAcTpkbWrVuH/v37Izg4GC+99BJ4AiEREREBTJg0Ll26hI0bNyI7OxunTp1CdnY2vvvuO0uHRURERFbA0dIBWJNbt27h5s2bAID6+np4enpaOCIiIiKyBjYzwpSRkYHo6Gj4+PhAkiQcOHBAp05ycjKCgoLg7OwMhUKBzMxMg9vv0aMHFi1aBH9/f/j4+OChhx5C7969TfgMiIiIyFbZzAhTTU0NBg4ciLlz52LatGk6j+/atQsJCQlITk7GyJEj8f777yMqKgq5ubnw9/cHACgUCtTW1ur87pdffomOHTvi0KFDyM/PR8eOHREVFYWMjAyMHj26yXhqa2u12qqqqjLRMyUiIiJrY5OXRpEkCfv378fUqVM1ZUOHDkVoaCjeffddTdn999+PqVOnIikpSW+be/bsQXp6OpRKJQBg9erVEELglVdeabL+W2+9heXLl+uU89IoREREtqNdXRqlrq4O2dnZiIyM1CqPjIxEVlaWQW34+fkhKysLN2/ehEqlQnp6Ovr27dts/SVLlqCyslJzKywsbNVzICIiIutlM1NyLamoqIBKpYKXl5dWuZeXF8rKygxqY9iwYZg4cSIefPBBODg4YPz48Zg8eXKz9eVyOeRyOZRKJZRKJVQqVaueAxEREVkvu0iYGkiSpHVfCKFT1pJ33nkH77zzzl31GR8fj/j4eM2QHhEREdkfu5iS8/DwgEwm0xlNKi8v1xl1MjWlUong4GCEhYWZtR8iIiKyHLtImJycnKBQKJCamqpVnpqaihEjRpi17/j4eOTm5uL48eNm7YeIiIgsx2am5Kqrq3H+/HnN/by8POTk5MDNzQ3+/v5ITEzEzJkzMXjwYAwfPhybNm1CQUEB4uLizBoX1zARERHZP5vZViA9PR0RERE65bNnz8bWrVsBqDeuXLVqFUpLSxESEoJ169Y1u4+SqRl6WiIRERFZD0M/v20mYbJ2TJiIiIhsT7vah8mSuOibiIjI/nGEyUQ4wkRERGR7DP38tplF30RE7Z5KBWRmAqWlgLc3EB4OyGSWjoqoXeCUXCtxSo6I2kRKChAYCEREADNmqH8GBqrLicjsOCVnIpySIyKzSUkBYmOBO9+uG65ksHcvEBPT9nER2QEu+iYisgcqFbBwoW6yBPxWlpCgrkdEZsOEiYjImmVmAkVFzT8uBFBYqK5HRGbDhKmVuIaJiMyqtNS09YjIKEyYWonXkiMis/L2Nm09IjIKEyYiImsWHg74+v62wPtOkgT4+anrEZHZMGEiIrJmMhmwYYP633cmTQ3316/nfkxEZsaEiYjI2sXEqLcO6NVLu9zXl1sKELUR7vTdSkqlEkqlEiqe0ktE5hQTA0yZwp2+iSyEG1eaCDeuJCIisj3cuJKIiIjIRJgwEREREenBhImIiIhIDyZMRERERHowYSIiIiLSgwlTK/FackRERPaP2wqYCLcVICIisj3cVoCIiIjIRJgwEREREenBS6OYSMPMZlVVlYUjISIiIkM1fG7rW6HEhMlErl27BgDw8/OzcCRERER0t65duwZXV9dmH+eibxO5ffs2SkpK0LVrV0iSZOlwrEZVVRX8/PxQWFjIxfB34LFpHo9N83hsWsbj0zwem6YJIXDt2jX4+PjAwaH5lUocYTIRBwcH+Pr6WjoMq+Xi4sI/0Gbw2DSPx6Z5PDYt4/FpHo+NrpZGlhpw0TcRERGRHkyYiIiIiPRgwkRmJZfLsWzZMsjlckuHYnV4bJrHY9M8HpuW8fg0j8emdbjom4iIiEgPjjARERER6cGEiYiIiEgPJkxEREREejBhIiIiItKDCRMRERGRHkyYyKIee+wxdO/eHbGxsVrlZ86cwaBBgzS3jh074sCBA5YJ0kKaOzYNrl+/joCAACxatKiNI7MOLR0ffceuPVm3bh369++P4OBgvPTSS3ovMNqeODo6at5j5s+fb+lwrEZhYSHGjh2L4OBgDBgwAHv27LF0SFaB2wqQRaWlpaG6uhoff/wx9u7d22Sd6upqBAYG4sKFC+jcuXMbR2g5+o7N0qVLce7cOfj7+2PNmjUWiNCyWjo+hryu2oNLly5h2LBh+Omnn9ChQweMHj0aa9aswfDhwy0dmlXw8PBARUWFpcOwOqWlpbh48SIGDRqE8vJyhIaG4syZM+3q/bcpHGEii4qIiEDXrl1brHPw4EGMHz++3f2xtnRszp07h19++QUTJ05s46isR0vHx5DXVXtx69Yt3Lx5E/X19aivr4enp6elQyIr5+3tjUGDBgEAPD094ebmhv/+97+WDcoKMGGiZmVkZCA6Oho+Pj6QJKnJKbHk5GQEBQXB2dkZCoUCmZmZJo9j9+7deOKJJ0zebmtY+tgsWrQISUlJJmvP1Cx9fGyFuY9Tjx49sGjRIvj7+8PHxwcPPfQQevfubcJnYD5t8RqqqqqCQqHAqFGjcOTIERNFbn5t+fd14sQJ3L59G35+fq2M2vYxYaJm1dTUYODAgdi4cWOTj+/atQsJCQlYunQpTp48ifDwcERFRaGgoEBTR6FQICQkROdWUlJiUAxVVVX49ttvrW4kxZLH5rPPPkOfPn3Qp08fkz4nU7KG144tMPdxunLlCg4dOoT8/HwUFxcjKysLGRkZbfX0WqUtXkP5+fnIzs7Ge++9h1mzZqGqqqpNnltrtdXf1+XLlzFr1ixs2rTJ7M/JJggiAwAQ+/fv1yobMmSIiIuL0yrr16+fWLx48V21nZaWJqZNm9bkY9u2bRNPP/30XbXX1tr62CxevFj4+vqKgIAA4e7uLlxcXMTy5cuNir0tWOq109Jj1sgcx2n37t3ihRde0NxftWqV+NOf/tTqWNuaOV9DDR555BFx/PhxY0O0GHMdm5s3b4rw8HCxbds2U4RpFzjCREapq6tDdnY2IiMjtcojIyORlZVlsn6scTpOH3Mfm6SkJBQWFiI/Px9r1qzBc889hzfffLPV7baVtnrt2DpTHCc/Pz9kZWXh5s2bUKlUSE9PR9++fc0RbpsyxbG5cuUKamtrAQBFRUXIzc3FPffcY/JY25opjo0QAnPmzMG4ceMwc+ZMc4RpkxwtHQDZpoqKCqhUKnh5eWmVe3l5oayszOB2JkyYgB9++AE1NTXw9fXF/v37ERYWBgCorKzEsWPHsG/fPpPGbm5tcWxsWVscH3s4dqY4TsOGDcPEiRPx4IMPwsHBAePHj8fkyZPNEW6bMsWx+fnnn7FgwQI4ODhAkiRs2LABbm5u5gi3TZni2Hz77bfYtWsXBgwYoFkftX37djzwwAOmDtemMGGiVpEkSeu+EEKnrCWHDx9u9jFXV1dcvHjR6NgszZzHpsGcOXPuNiyrYc7jY8ixsxWtPU7vvPMO3nnnHVOHZRVac2xGjBiBU6dOmSMsq9CaYzNq1Cjcvn3bHGHZNE7JkVE8PDwgk8l0vrGUl5frfLNpb3hsWsbjYxgep+bx2DSPx8Z8mDCRUZycnKBQKJCamqpVnpqaihEjRlgoKuvAY9MyHh/D8Dg1j8emeTw25sMpOWpWdXU1zp8/r7mfl5eHnJwcuLm5wd/fH4mJiZg5cyYGDx6M4cOHY9OmTSgoKEBcXJwFo24bPDYt4/ExDI9T83hsmsdjYyEWPEOPrFxaWpoAoHObPXu2po5SqRQBAQHCyclJhIaGiiNHjlgu4DbEY9MyHh/D8Dg1j8emeTw2lsFryRERERHpwTVMRERERHowYSIiIiLSgwkTERERkR5MmIiIiIj0YMJEREREpAcTJiIiIiI9mDARERER6cGEiYiIiEgPJkxEREREejBhIiJqxtKlSyFJEiIjIy0dChFZGBMmIqJm/PDDDwCA0NBQC0dCRJbGhImIqBknT54EwISJiJgwERE1qbi4GBcvXgTAhImImDARETWpYTrOxcUFvXv3tnA0RGRpTJiIiJrQkDA9+OCDkCQJX375JWJiYuDj44POnTtjwIABUCqVuH37toUjJaK24GjpAIiIrFFDwvTAAw/gmWeewSeffAKZTIauXbvi+vXrOHXqFF588UVcuXIFr7/+uoWjJSJzk4QQwtJBEBFZGz8/PxQVFcHd3R1OTk5Ys2YNHnvsMXTs2BFFRUWYMWMGMjMz0aVLF1RVVUGSJEuHTERmxISJiOgOFRUV6NGjBwDA3d0d2dnZCAgI0Kpz/PhxDBkyBABQWFgIX1/fNo+TiNoO1zAREd2hYToOAD744AOdZAlQJ1INZDJZm8RFRJbDhImI6A4NCZO3tzemTJnSZJ1z584BADp37oyePXu2WWxEZBlMmIiI7tCQMEVHR8PBoem3yZycHADAgAEDuH6JqB1gwkREdIeGhCksLExvnQcffLBNYiIiy2LCRETUSGVlJX799VcALe/wzevMEbUvTJiIiBo5efIkhBBwcnJCSEhIk3Wqqqrwn//8BwBHmIjaCyZMRESNNIwc9e/fH05OTk3WycnJgRACHTp0aDapIiL7woSJiKgRQ6baDEmqiMi+MGEiImrkbhImrl8iaj+YMBER/c/169dx9uxZAC0nQydPngTA9UtE7QkvjUJERESkB0eYiIiIiPRgwkRERESkBxMmIiIiIj2YMBERERHpwYSJiIiISA8mTERERER6MGEiIiIi0oMJExEREZEeTJiIiIiI9GDCRERERKQHEyYiIiIiPZgwEREREenBhImIiIhIj/8PGZWvSNQLPv4AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "plt.title(\"Accuracy of the numerical derivative\")\n", "plt.xlabel(\"${h}$\", fontsize=18)\n", "plt.ylabel(\"relative error\", fontsize=18)\n", "plt.xscale('log')\n", "plt.yscale('log')\n", "plt.scatter(arr_h, arr_err, color=\"red\")\n", "\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.16" } }, "nbformat": 4, "nbformat_minor": 4 }