{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In lecture 21 of mcs 320, we work with functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Expressions in Many Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the definition of the following sequence of polynomials."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "   f_i = x_i^3 + \\left(\n",
    "     \\sum_{ \\scriptsize \\begin{array}{c}\n",
    "                j=1 \\\\\n",
    "                j \\not=i\n",
    "            \\end{array} }^n (i + j + 1) x_j \\right) - 1,\n",
    "   \\quad i=1,2,\\ldots, n.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us look at a specific value for $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As it should work for any number of variables, we define a list of variables ``X``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[x1, x2, x3, x4]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = [var('x%d' % k) for k in range(1, n+1)]\n",
    "X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Keep in mind that lists start with index 0, so ``X[0]`` is ``x1``."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us define a function for the first expression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function depends on ``i`` and ``n``.  Instead of defining the function with ``f(i, n) = ``, we must us ``lambda`` to prevent evaluation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "x1^3 + 4*x2 + 5*x3 + 6*x4 - 1"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f = lambda i, n: X[i-1]^3 + sum((i+j+1)*X[j-1] for j in range(1, n+1)) - (2*i+1)*X[i-1] - 1\n",
    "f(1, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[x1^3 + 4*x2 + 5*x3 + 6*x4 - 1,\n",
       " x2^3 + 4*x1 + 6*x3 + 7*x4 - 1,\n",
       " x3^3 + 5*x1 + 6*x2 + 8*x4 - 1,\n",
       " x4^3 + 6*x1 + 7*x2 + 8*x3 - 1]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F = [f(i, 4) for i in range(1, 5)]\n",
    "F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Polynomials in Several Variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A polynomial in $n$ variables $x_1, x_2, \\ldots, x_n$ with $m$ terms is defined by\n",
    "\n",
    "1. a list of $m$ coefficients, let us take rational numbers,\n",
    "\n",
    "2. a list of $m$ exponent tuples, each tuple is defined by $n$ natural numbers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the example below, we take $n = 5$ and $m = 8$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "(n, m) = (5, 8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[x1, x2, x3, x4, x5]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = [var('x%d' % k) for k in range(1, n+1)]\n",
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "set_random_seed(20230707)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coefficients are stored in the list ``C``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, -2, 1/2, -1/2, -6, 0, -3, -7/3]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C = [QQ.random_element() for _ in range(m)]\n",
    "C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(1, 0, 1, 1, 2),\n",
       " (1, 0, 4, 1, 1),\n",
       " (4, 0, 2, 1, 1),\n",
       " (3, 7, 0, 1, 2),\n",
       " (0, 1, 0, 36, 1),\n",
       " (4, 1, 0, 1, 36),\n",
       " (0, 1, 0, 0, 1),\n",
       " (1, 8, 1, 2, 13)]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E = [tuple([abs(ZZ.random_element()) for _ in range(n)]) for _ in range(m)]\n",
    "E"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use ``prod`` to make the first monomial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "x1*x3*x4*x5^2"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prod(X[k]^E[0][k] for k in range(n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is good to have a helper function to define a monomial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "x1*x3*x4*x5^2"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "monomial = lambda x, e: prod(x[k]^e[k] for k in range(len(x)))\n",
    "monomial(X, E[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = lambda x, c, e: sum(c[k]*monomial(x, e[k]) for k in range(len(c)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-6*x2*x4^36*x5 - 7/3*x1*x2^8*x3*x4^2*x5^13 - 1/2*x1^3*x2^7*x4*x5^2 + 1/2*x1^4*x3^2*x4*x5 - 2*x1*x3^4*x4*x5 - 3*x2*x5"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p(X, C, E)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To verify, let us look at the coefficients and the exponents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0, -2, 1/2, -1/2, -6, 0, -3, -7/3]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(1, 0, 1, 1, 2),\n",
       " (1, 0, 4, 1, 1),\n",
       " (4, 0, 2, 1, 1),\n",
       " (3, 7, 0, 1, 2),\n",
       " (0, 1, 0, 36, 1),\n",
       " (4, 1, 0, 1, 36),\n",
       " (0, 1, 0, 0, 1),\n",
       " (1, 8, 1, 2, 13)]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Function Composition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the trajectory of a projectile in space modeled by a parabola,\n",
    "subject to the following constraints.  At time $t = 0$ it is launched from\n",
    "the ground and at time $t = 45$ it hits the ground 120 miles further.\n",
    "Assuming constant horizontal speed we create a function $f(t)$ to give\n",
    "the altitude of the projectile in function of time.\n",
    "First we model the shape of the trajectory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAADHCAYAAAAJUBeIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA2G0lEQVR4nO3deXiU1dnH8e9NIGGPqCAoWNQioGBBCgaKgIAbCoK2LihV64ra2mqVorVVW6UUUQqiFaoiuGDdEGkRFMENROStCqKgghoEWQQnrAmE8/5xnpDZyJ6ZTPL7XNdzxZznzMyZY0juOct9zDmHiIiIiMRXK9kNEBEREanKFCyJiIiIFEHBkoiIiEgRFCyJiIiIFEHBkoiIiEgRFCyJiIiIFEHBkoiIiEgRqm2wZF5jM7Nkt0VERERSV+1kN6CMis2kGQqFyMzMJBQKJaI9IiIikjpKNZCSqsGSiFRnzkEoBDk5sG2bv7ZvL/y6fTvk5cGePbB3b+TX/HxIS4PatWO/pqdDw4bxr8aN4ZBDoE6dZL97EaliFCyJSOLk5sLatZCdXXht2ACbNvlr40b/dfNmH/gkQ2YmNG0Khx5a+LVZM2jZElq1KrwOPRQ0yy9SIyhYEpGK45wPdD7/HFat8tfnn8NXX8E33/hgqKoLhfz1xRdF16tb1wdQRx4JP/4xHHsstGnjvx59tB/FEpFqwUpzkK6ZDQeGA62Dok+Au51zs4P7U4BLox622DmXFfYcGcB9wEVAPWAecJ1zbm1YnSbAeGBQUDQT+LVz7ofg+wM2euLEiUycOJH8/HxWrVpFKBSicePGJX6PIlJCGzbAxx8XXp9+6gOjH35IdsuSr1YtaN3aB08dOsAJJ/irfXvIyEh260SklGuWShssDQTygYKPXJcCtwCdnXOfBMHSYcDlYQ/Lc85tCXuOh4GBwGXA98BY4GCgi3MuP6gzG2gJXB08bBLwlXNuYPB9sY3OycnZv8BbwZJIOezb54Og99+HDz8sDI4SOUpUv75fV9Sokb8aNvRBR+3afo1R+Ne0NL9uKT/fr2Mq+Lp3L+zeDTt2FK57Krj27k3M+0hLg3btCoOnn/7UXwcdlJjXF5EClRcsxX0Csy3ALc65R4Ng6SDn3OAD1M0ENgHDnHPPBmWHA9nAAOfcHDNrD6wAspxzi4M6WcAioJ1zbiUKlkQqz4YNPjB6/31YvBiWLKn40aJGjfz0VatW0KKFXxsU72rSxAdGaWkV+/rR8vJg61b4/vvCNVPhX9evL1xjtX69DyArUrt20K0bnHSS/3rCCZrGE6lcidkNZ2ZpwC+ABvhApkAfM9sI/AC8CdzunCv4CNoFqAPMLajsnFtnZsuBHsAcoDsQKgiUgjrvmVkoqLMyXntyc3PJzc3d/31OTk5Z35pIzeGcX0/05pv+eustWL26/M9bpw4cc0zhOp4f/9gHRwUBUmZm+V+jIqWnw2GH+as4e/dGBk+rV/uRt4J1WmUZcfvsM39Nneq/z8jwgVPv3tCrF3TvDg0alP55RaRClDpYMrOO+OCoLrAdGOKcWxHcng08B3wNHAX8BXjDzLo453KB5vhpua1RT7shuEfwNd5vm41hdWKMGjWKu+66q7RvR6RmcQ6+/BLmzy8MkNauLf5xB1K/fuGanA4d/AhJmzY+KKpdTfeP1K5duCMunlCoMHhasQKWLfPTlmvWlPw1cnN94PrWW4Wv2bWrD5wKAigFTyIJU+ppODNLB44EDgLOA64EeocFTOF1W+ADpwudcy+a2VDgcedcRlS914AvnXPXmtltwKXOubZRdT4HHnXO/Y0403DxRpZatWqlaTiRbdt8cPTqqzBnTtlHjlq08H+wO3UqXHNz9NGVP0VWXeTkwPLlhWu+li71a8Dy8kr/XOnpcPLJcPrp/urYUWkMREqncqfhnHN5FC7w/sDMugI3AtfEqbvezL4G2gRF3wHpZtYkanSpGbAwrE68sfCm+BGouDIyMsjQLhMRP3r00UeFwdG775Y+Z1GDBj4w6tatcC3NEUfoD3J5NG4MPXr4q0Burg+cFi8uXCO2alXxz5WXB/Pm+evWW30ge9ppPnA64wy/1ktEKkxFLPCeB2Q75y6Lc+8Q4Fvgaufc1LAF3pc45/4d1GkBrCV2gfdJzrn3gzonAe+hBd4i8e3ZA2+/DTNmwMsv+5xGpdG0aeQUT4cOGjFKlq1bYeHCwmnSpUv9jr6SSkuDPn3gnHP8deSRldZUkRRWqakD7sWvS8oGGgEXAn8AzsCvY7oTeAFYj8/FdC9+yq69c25b8BwPA2fjUwdswedcOoTY1AGHUzhaNQn4uiSpA5RnSWqMHTv8yNGMGTBrlv8jW1KHHQannFIYHLVvr1Gjqmr7dh88vfUWLFjgR59Kk+rgxBN90DR4sKbrRApVarD0KNAPaAGEgI+B0c6518ysHjAD6Ixfz7QemA/c4ZzLDnuOusAYYCiRSSnD6xxMbFLKG0qSlLKARpakWtq1C2bPhunTfYC0a1fJHlenDvTsWThNc8IJ+qOZqnJy4I03CqdZv/qq5I899li48EJ/tW9faU0USQGJzbOUJAqWpObYswdee80HSDNm+AXbJXHUUXDmmT44OuUUn69Iqhfn/K67OXN88DRvnl8HVRInnOCDpgsu8Av1RWoWBUugYElSnHN+umXKFHjuOdiypdiHANCli59uGTwYjj9eo0c1zfbtMHeuX7f2yisln5rt1g2GDYOLLoJDDqncNopUDTU7WNKaJUlp69b5xIRTpsDKuPlXIxUs5h082K9LOVDuH6l59u6Fd97xo5EzZsDXXxf/mPR0GDQILrvMT9lW11xZIjU9WCqgkSVJGbt3w8yZPkCaM6f4ozTMfIB00UVw7rkaCZDiOeePrZk+HZ591gflxWne3I82XXYZHHdcpTdRJMEULIGCJUkBn38OjzwCjz9esmm2rCy/xuT8831eHZGy2LfPp5mYPh2ef96ff1ecHj1g+HD4+c+hbt3Kb6NI5VOwBAqWpIrau9evJXn4Yb9ouzjHHOM/2V98sV+wLVKR9uzxi8KfeAJeeqn4xeGHHgqXXw7XXON/NkVSl4IlULAkVcy6dTB5sr++/bboug0b+tGjyy+Hn/1Mi7QlMbZu9VN0jz/us4kX5/TT/WjT2WcrgamkopodLGmBt1QpH3wADzwA//538YkE+/TxAdJ55+mQVEmuTz7xo01Tp8KGA54y5bVuDTfeCL/6lT/SRSQ11OxgqYBGliRp8vP9VNsDDxSeGn8gTZoUTmsce2xi2idSUnv2+DQEDz/sE2EWpVEjuPJK+M1vfAAlUrUpWAIFS5IE27f7KYx//AO+/LLoullZcO21frqtXr3EtE+kPFauhH/+0+/a/OGHA9erVQuGDIGbboLu3TWNLFWVgiVQsCQJtHkzjBsHEycW/UekXj2/Ffvaa6Fz50S1TqRi7dzpp5UnTvTTzEXp0QNuuw0GDFDQJFVNqX4ga5Xqmc2Gm9nHZpYTXIvM7Myw+2Zmd5rZOjPbZWYLzOz4qOfIMLMJZrbZzHaY2UwzaxlVp4mZTTOzUHBNM7ODStLGiRMnctxxx9G1a9fSvDWR0lu7Fn73O/jRj+Ceew4cKB1+ONx7r6//yCMKlCS11a/vd2i+/75PQXDuuX40KZ6FC/0C8M6d/eLx/PyENlWkopT2IN2BQD7wRVB0KXAL0Nk594mZjQBuBy4DVgF/BHoBbZ1z24LneBgYGNT5HhgLHAx0cc7lB3VmAy2Bq4PXmQR85ZwbGHyvkSVJni++gL//3U9H7Nlz4HqdO/upiPPP95mRRaqr1ath/Hh49FE/HX0gbdrAiBF+hFX/JiS5EjsNZ2Zb8AHTY8A6YJxzbnRwLwPYAIxwzj1iZpnAJmCYc+7ZoM7hQDYwwDk3x8zaAyuALOfc4qBOFrAIaOecW4mCJUmGTz7xI0TTpx84y7YZDBzog6RevTT1IDVLKOQDpnHjIDv7wPVatoRbb4WrrlKSS0mWypuGi3gVszQzuxBogA9kjgKaA3ML6jjncoE3gR5BURegTlSddcDysDrdgVBBoBTUeQ8IhdWJkZubS05OTsQlUiFWrYKhQ6FjR3j66fiBUloa/PKXPqB6+WXo3VuBktQ8mZn+g8KXX/qR17Zt49dbu9bvmmvTxi8az8tLaDNFSqvUwZKZdTSz7UAu8E9giHNuBT5QAj+SFG5D2L3mQJ5zLvoo7Og6G+O89MawOjFGjRpFZmbm/quVDhSV8lqzxm/rb98ennnGn68VLSPDJ+b74gufl6Z9+8S3U6SqqVMHLr3Uf3h4/nk48cT49dau9f9+2rSBf/2r6GltkSQqy8jSSqATkAU8DDxhZuGnLEb/RbE4ZdGi68SrX+TzjBw5klAotP/KLmoIWKQo2dmFeY+mTIk/ktSgAfz+9z6geugh5ZURiSctzSdZ/eADf0h0797x633zjZ+Sa9fOf+goLoGrSIKVOlhyzuU5575wzn3gnBsJfATcCHwXVIke/WlG4WjTd0C6mTUpps5hcV66KbGjVvtlZGTQuHHjiEukVDZu9JmIf/xjmDQp/i/sRo3gjjv8L/cxY3SgrUhJmMFpp8GCBX4HXd++8eutXu132nXo4M+qS83UNlINlXnNUhgDMoA1+EDn1P03zNKB3sDCoGgpsCeqTgugQ1idRUCmmXULq3MSkBlWR6Ti7Nzpt/7/+Md+R0+89RP168Mf/uBHku6+Gw4+OPHtFKkOevb0h/fOn+//O56VK31Kgp494d13E9s+kThKm2fpXjM72cxaB2uX7gH6AE85v61uHHCbmQ0xsw7AFGAn8DSAcy4EPAqMNbN+ZtYZeBJYBrwe1PkUeBWYbGZZwU64ycCsYCdckZRnSUosP9/v3GnTBv74R9i2LbZO3bp+weqaNTBqFBxySOLbKVId9enjjwOaOxdOOil+nYULfcA0ZAh89llCmycSrrR5lh4F+gEt8LvTPgZGO+deC+4b8GfgGqAJsBi43jm3POw56gJjgKFAPWAecJ1zLjuszsHAeGBQUDQTuME590PwvVIHSNk5B7Nn+63Ln3wSv06dOnD11T778OGHJ7Z9IjVNwb/JO+6A//u/+HXS0uCKK+DOOzX9LRVBx52AgiU5gP/9D26+2U8BxGPmE+bdfbfPzC0iibNvHzz3nP+Qsnp1/Dr16/v7N92kcxWlPBKTZ0kkpWza5He4dely4EDptNP8p9onnlCgJJIMtWrBBRfAp5/6A6kPPTS2zs6dftr8uOPghRe0CFwSotoFS1qzJBH27PHZhNu08Tvc4v1i/clP/LbmOXOgU6dEt1BEoqWn+6SVX34Jt98efwTpq6/g5z/3O+s++ijhTZSaRdNwUn3NnQu//a3/lBpPq1Z+F9zFFx/4IFARSb5vv/VrlR59NP4Hnlq1/BrDu++Gpk0T3jxJSZqGkxpu9WoYPBhOPz1+oFS/Pvz1r3578rBhCpREqrojjoDJk/00ea9esff37fPHprRpAw8+6He6ilQg/ZWQ6iM31wdBxx/vz2eLZ+hQHyQdaGhfRKquTp18Ysvnnou/rjAUgl//2qciWLIk0a2TakzBklQP8+f7tUd33AG7d8feP/FEeOcdeOopf+K5iKQmM79W6dNP/bRbvA89S5f6gOm66+CHHxLeRKl+ql2wpAXeNcyGDX4qrW9fP2IUrWlTf0Dn++/Dz36W+PaJSOWoV89/OFq50o8YR3MOHn4Y2raFJ5/UrjkpFy3wltS0b5/f3TZyZPxPjrVqwQ03wF13wUEHJbp1IpJob74Jw4cfeENHnz7+0Ov27RPaLKmytMBbqrnly6FHD/+LMV6g1LWrX6/wj38oUBKpKXr3hg8/9McSxZuaW7DAr3n6y1/in/8oUgQFS5I68vL8SNGJJ8LixbH3MzP9J8dFi3wdEalZ0tP9gdcrVsDAgbH38/LgT3+Cn/5UC8ClVEp7kO5IM1tiZtvMbKOZzTCztlF1ppiZi7rei6qTYWYTzGyzme0ws5lm1jKqThMzm2ZmoeCaZmYHFddGrVmqppYs8b/g7rzTJ5qMNnSoP2hz+HB/hpSI1FytW8PMmX5X7JFHxt5ftgyysuD3v/cZwUWKUdqDdF8FpgNLgNrAPUBH4Djn3I6gzhTgMODysIfmOee2hD3Pw8BA4DLge2AscDDQxTmXH9SZDbQErg4eNgn4yjk3EK1Zqjl27oQ//xnuv9+vU4rWpo1fxNmvX+LbJiJV344dfkR67Nj4v0OOPtpvAjnllMS3TZIpcQfpmllTYCPQ2zn3VlA2BTjIOTf4AI/JBDYBw5xzzwZlhwPZwADn3Bwzaw+sALKcc4uDOlnAIqCdc+6z4tqmYKkaePNNuPJK+OKL2HtpaXDLLX5IXfmSRKQ4H3wAV1wBH38c//5VV8GYMX46X2qChC7wLvip2hJV3ieYpltlZpPNrFnYvS5AHWBuQYFzbh2wHOgRFHUHQgWBUlDnPSAUVidCbm4uOTk5EZekqB07/E62Pn3iB0o/+Ylfs3SghZwiItF++lMfMP31r35tU7TJk6FjR5g3L/FtkyqvzMGSmRlwP/COc2552K3ZwMVAX+BmoCvwhpllBPeb46fltkY95YbgXkGdjXFedmNYnQijRo0iMzNz/9WqVauyvC1JtkWL/I6ViRNj76Wn+7PcliyBLl0S3jQRSXF16vjs/R995HfURsvOhv79fRZwrWWSMOUZWXoQOAG4KLzQOfesc+4/zrnlzrlXgDOBY4Gzink+I3ItUrz5weg6+40cOZJQKLT/ys7OLun7kKogL8//EuvZM/5oUvfuflvwbbf5X3giImXVrh28/TZMmAANG8bef/BB/6Htvfdi70mNVKZgycwmAIOAU5xza4uq65xbD3wNtAmKvgPSzaxJVNVm+NGlgjqHxXm6pmF1ImRkZNC4ceOIS1LEsmXQrRvce2/sAsx69WDcOP+LTcnkRKSiFCSuXb7cnwAQ7fPPfdb/227z505KjVba1AFmZg8C5wJ9nXNrSvCYQ4BWwPqgaCmwBzg1rE4LoAOwMChaBGSaWbewOifh10gV1JFUl58Pf/+7X0vw0Uex97Oy/GjSjTcqHYCIVI4f/Qheew3Gj49dA7lvn18b2a1b/N9RUmOUdmRpInAJMBTYZmbNg6segJk1NLP7zKy7mbU2sz7AK8Bm4CUA51wIeBQYa2b9zKwz8CSwDHg9qPMp8Cow2cyygp1wk4FZzrk4B4CFNVB5llLD11/7BdwjRsRm061Tx69NevttOPbYpDRPRGqQWrX8OqUPP/QH8Eb7+GMfMD3wQPz0A1LtlTbP0oEqX+6cmxIETTOAzsBB+NGk+cAdzrn9i4jMrC4wBh901QPmAddF1TkYGI+f7gOYCdzgnPsB5VlKbc8+C9dcA6FQ7L0OHWDaNL9eQEQk0fbu9SkE/vzn+AlwTz8dpkyB5nH3GknqSFyepSRSsJSKtm+H3/wGHn889p6Zz5t0992QkRF7X0QkkT76CIYN82sqozVt6gOmAQMS3iypMDpIV6qgpUv9eW3xAqWjjoK33oLRoxUoiUjV8JOf+DQlt94ae2/TJjjrLL+ecvfuxLdNEk7BklSuffvgvvv81v/PP4+9f8klfp1Az54Jb5qISJEyMvyHuNdegxYtYu+PH+/XMn3ySeLbJglV7YIlLfCuQr77Ds44w0+vRc/9N2zo1yZNmwaaJhWRqqx/f7/Ie9Cg2HvLlvkdvY88Aqm5rEVKQGuWpHK88QZcdBFsjJOIvVs3ePppOOaYxLdLRKSsnIN//hNuuin+9NtFF/mgqVGjxLdNSktrliSJ8vP9Iu3+/WMDJTMYORLeeUeBkoikHjMYPtyfMdexY+z9Z57xo0zxFoVLSlOwJBVn40Y480y/5TZ6xPLww+H1132Wbh1XIiKp7Pjj4f33fW6maKtW+dHzxx7TtFw1Uu2CJa1ZSpK334bOnf1CyGgDBvhtuPGOFBARSUV16/oF3i++CJmZkfd274YrroDLLoMdO5LSPKlYWrMk5bNvn0/gdvvtfgouXFqaz8R9yy0+Q66ISHW0ejWcf75PkRLtuOPguef8V6lKtGZJEmTrVjjnHPjDH2IDpRYt/CLvESMUKIlI9Xb00fDuu/5g3mgrVkDXrjB9euLbJRVGf8WkbJYt878AZs2Kvde/v8+d1KtXwpslIpIUGRkwYQL8+9+xu+F27vQ75W6+2R+nIimnVMGSmY00syVmts3MNprZDDNrG1XHzOxOM1tnZrvMbIGZHR9VJ8PMJpjZZjPbYWYzzaxlVJ0mZjbNzELBNc3MDiqujVqzlADPPgtZWfDll5HlZnDnnfDqq9CsWVKaJiKSVL/4hZ+Oi3e+5f33w6mnxk+pIlVaaQ/SfRWYDiwBagP3AB2B45xzO4I6I4DbgcuAVcAfgV5AW+fctqDOw8DAoM73wFjgYKCLcy4/qDMbaAlcHbz8JOAr59xAtGYpOfbu9VNuY8fG3mva1OdO6t8/8e0SEalqdu/2Z2FOnhx7r2VLeOEFv2tOkiVxB+maWVNgI9DbOfeWmRmwDhjnnBsd1MkANgAjnHOPmFkmsAkY5px7NqhzOJANDHDOzTGz9sAKIMs5tziokwUsAto55z4rrm0KlirYpk1w4YV+HVK0bt38P/yWLWPviYjUZP/6F1x/PeTlRZanp8NDD/ldc5IMCV3gXbBfckvw9SigOTC3oIJzLhd4E+gRFHUB6kTVWQcsD6vTHQgVBEpBnfeAUFidCLm5ueTk5ERcUkGWLvWJ1uIFSlde6Q/BVaAkIhLryit9apXo35F5ef7eNddAbm5y2iYlVuZgKRhFuh94xzm3PChuHnzdEFV9Q9i95kCec25rMXXiTepuDKsTYdSoUWRmZu6/WrVqVfI3Iwf2xBPws5/BN99Elqenw6RJfog5IyM5bRMRSQXduvkPnX36xN6bNAl694b16xPeLCm58owsPQicAFwU51703J7FKYsWXSde/QM+z8iRIwmFQvuv7OzsYl5OipSfD7//vU+qFv2p54gj4M034aqrktI0EZGU06yZT9p7002x9xYv9ruL4+VpkiqhTMGSmU0ABgGnOOfWht36LvgaPfrTjMLRpu+AdDNrUkydw+K8dFNiR60AyMjIoHHjxhGXlFFOjj9dO95C7pNP9v+gs7IS3y4RkVRWu7b/vfrMM1C/fuS9b7/1v1+ffTY5bZMilTZ1gJnZg8C5QF/n3JqoKmvwgc6pYY9JB3oDC4OipcCeqDotgA5hdRYBmWbWLazOSfg1UgV1pDJ8+SV07w7//W/svd/8BubNg8PixbEiIlIiF14I773nk1mG27XL3/vTn/zpCFJllDZ1wEPAUOAcYGXYrZBzbldQZwQwErgc+By4DehDbOqAs/GpA7YA9wGHEJs64HDgmuA1JgFfK3VAJVqwAM47D7ZsiSyvUwceeQQuvzwpzRIRqZY2b/Z5mRYsiL137rkwdSo0aJDwZtUQlbobbjh+dGcBsD7suiCszt+BccBDwAfAEcBpBYFS4HfADODfwLvATmBgQaAUuBhYht81Nxf4GBhWXAOVlLKMJk3yydKiA6VDD/W74BQoiYhUrEMPhblz/Y64aC++6DfXfP114tslMXSQbk23d69fcDhhQuy9jh1h5kxo3TrhzRIRqTGc8zmXbrwx9pzNZs3gpZegR9ysOVJ2OkhXSignB84+O36gNGiQPxhSgZKISOUy84krX30VDjoo8t7GjdC3rxZ+J5mCpZoqOxt69oQ5c2Lv/eEP/pNM9GGQIiJSefr3h/ffh7ZtI8tzc/3C77/9zY9CScJVu2BJa5ZK4MMP/db/ZcsiyzMyYNo0GDUKalW7Hw0RkaqvTRu/U+7002PvjRwJV18Ne/Ykvl01nNYs1TSzZ8P558P27ZHlTZvCyy/7tAEiIpJce/f6dC0PPxx777TT4LnnQH/bykNrluQA/vlPGDgwNlBq29Z/klGgJCJSNdSuDRMnwn33+TVN4ebO9csodFJFwihYqgn27YNbb4Xhw2N3WvTuDQsXxiZHExGR5DKDm2/2o0h160beW7YMTjoJ/u//ktO2GqbaBUtasxSlICPsmDGx9y65xC/wPvjgxLdLRERK5rzzYP58v1wi3Pr10KsX/Oc/yWlXDaI1S9XZli2FKQCi3XEH3HVX7PCuiIhUTatXw1lnwWefRZanpfnEwr/6VXLalZq0ZkmAtWv9oYzRgVLt2vDYY3D33QqURERSydFH+2UTffpElufnwxVXwD33KLVAJVGwVB19+qnP9rpiRWR548Z+N5yOLhERSU1NmvjlE8PinP71xz/Cr38duzZVyq3aBUs1fs3SokXxd0kccQS8845PeiYiIqkrPR2eeMInEI42caJfp5qbm/h2VWOlDpbMrJeZvWJm68zMmdngqPtTgvLw672oOhlmNsHMNpvZDjObaWYto+o0MbNpZhYKrmlmdlBx7bv++utZsWIFS5YsKe1bS32zZkG/frGH4bZr54duO3ZMTrtERKRimfkEwuPGxd57/nk44wwIhRLerOqqLCNLDYCPgBuKqPMq0CLsGhB1fxwwBLgQ6Ak0BGaZWVpYnaeBTsAZwdUJmFaG9tYMjz8Ogwf73W/hsrL8iNKRRyalWSIiUoluvBGeeQbq1IksX7DAp4ZZvz4pzapuyrUbzswcMMQ5NyOsbApwkHNu8AEekwlsAoY5554Nyg4HsoEBzrk5ZtYeWAFkOecWB3WygEVAO+fcZ/GeO1yN2Q3nHIwe7dPgRzvrLH/4YoMGiW+XiIgkzrx5MGQIbNsWWd66tV/jdOyxSWlWFVYldsP1MbONZrbKzCabWbOwe12AOsDcggLn3DpgOdAjKOoOhAoCpaDOe0AorE6E3NxccnJyIq5qb98+uOmm+IHSpZf6w3AVKImIVH/9+vnRpGbNIsu/+sqvY/3f/5LRqmqjMoKl2cDFQF/gZqAr8IaZZQT3mwN5zrmtUY/bENwrqLMxznNvDKsTYdSoUWRmZu6/WrVqVc63UcXt3eu3isabrx4xwk/LRQ/LiohI9XXiiX596jHHRJZv2gSnnBI/556USIUHS865Z51z/3HOLXfOvQKcCRwLnFXMQ43IZJPx5gej6+w3cuRIQqHQ/iu7Op+Zk5cHF10EU6bE3nvgAfjb35RDSUSkJjrmGB8wdekSWR4Kwamn+ik5KbVKTx3gnFsPfA20CYq+A9LNrElU1Wb40aWCOofFebqmYXUiZGRk0Lhx44irWtq5E845x+92CFe7Njz1FPz2t0lploiIVBHNmvnjUU45JbJ81y5/mPoLLySnXSms0oMlMzsEaAUULMlfCuwBTg2r0wLoACwMihYBmWbWLazOSUBmWJ2aJyfHbwd99dXI8owMvz5p6NDktEtERKqWRo38mXEDB0aW79kD55/v8zRJiZUlz1JDM+tkZp2CoqOC748M7t1nZt3NrLWZ9QFeATYDLwE450LAo8BYM+tnZp2BJ4FlwOtBnU/x6Qcmm1lWsBNuMjDLObeyqPZV26SUmzdD377w9tuR5Q0b+qzcZ5+dnHaJiEjVVK+eH0W66KLI8n374LLLYPz4pDQrFZU6dUAQAM2Pc+sJYDgwA+gMHIQfTZoP3OGc27+IyMzqAmOAoUA9YB5wXVSdg4HxwKCgaCZwg3PuB2raQbrr1vm55ujjS5o08YHSSSclp10iIlL15efD9dfDI4/E3rv7bn9MSs1b51qqN1yuPEtJVHOCpTVr/BElq1dHlh92GLz2mrJyi4hI8Zzzx6P8/e+x926+GcaMqWkBU5XIsyQVYeVKOPnk2EDpyCP9dJwCJRERKQkzv1P63ntj740d6w/g3bcv8e1KEdUuWKo2a5ZWrPCp6r/9NrL82GP98SVt2sR/nIiISDxmPonxgw/G3ps4EYYPV8B0AJqGq4qWLfPZWDdtiiz/yU98jozD4mVVEBERKaFp0/wi7+jg6PLLYfJkSEuL+7BqRNNwKe3DD31ujOhAqVs3nzdDgZKIiJTXsGEwfXpsUPT44z6I2rs3Kc2qqhQsVSVLl/r0AN9/H1nevTvMnet3v4mIiFSEX/wCnnsu9misJ5/0wZQCpv2qXbCUsmuW3n/fT71tjToy7+ST/dRbZmZy2iUiItXXkCE+F1N6emT59Olw4YU+iaVozVKVsGiRz8ydkxNZ3qcPzJoFDRokpVkiIlJDzJ7tA6fc3Mjyc86BZ5/1J0VUL1qzlFLefhtOOy02UOrXz6eqV6AkIiKV7cwz/YfzevUiy19+Gc47LzaIqmEULCXTO+/4H9Dt2yPLTz8dXnkF6tdPTrtERKTm6d8f/vvf2L89//kP/PznkJeXnHZVAQqWkmXRIh8o7dgRWT5gAMyYERvdi4iIVLY+ffxh7Q0bRpbPmgUXXFBj1zCV5SDdXmb2ipmtMzNnZoOj7puZ3Rnc32VmC8zs+Kg6GWY2wcw2m9kOM5tpZi2j6jQxs2lmFgquaWZ2UHHtS4kF3u+/70ePokeUBg2CF1+EunWT0y4REZGTT/Y7sBs1iiyfMcMfylsDA6ayjCw1AD4CbjjA/VuBm4L7XYHvgNfMLLzXxwFDgAuBnkBDYJaZhSd8eBroBJwRXJ2AacU17vrrr2fFihUsWbKkxG8ooZYu9WuUtm2LLB840G/hrH6L6EREJNV07+4XfUevm33hhRqZVqBcu+HMzAFDnHMzgu8NWAeMc86NDsoygA3ACOfcI2aWCWwChjnnng3qHA5kAwOcc3PMrD2wAshyzi0O6mQBi4B2zrnPimtbldwN9+GHPo9SdHqAAQP8iJICJRERqUreessvGdm5M7J86FCYOjWVM30ndTfcUUBzYG5BgXMuF3gT6BEUdQHqRNVZBywPq9MdCBUESkGd94BQWJ0Iubm55OTkRFxVyrJlfvFcdKB02mk+UlegJCIiVU2vXn6Bd/Q62qefhl/9CvLzk9OuBKvoYKl58HVDVPmGsHvNgTzn3NZi6myM8/wbw+pEGDVqFJmZmfuvVq1albrxlWbFCp8KIDozd79+fg5Ya5RERKSq6tMHZs6M/Vs1dSpcfXWNOHy3snbDRc/tWZyyaNF14tU/4POMHDmSUCi0/8rOzi5xYyvVypV+6i36rLfevf0Pn3a9iYhIVde/v/9wH53p+7HH4Nprq33AVNHB0nfB1+jRn2YUjjZ9B6SbWfRBZ9F14p0Y25TYUSsAMjIyaNy4ccSVdF9+6QOlDVFN7tnTb8NUHiUREUkVp58OL70Ue5bc5Mnw299Cap4IUiIVHSytwQc6pxYUmFk60BtYGBQtBfZE1WkBdAirswjINLNuYXVOAjLD6lRta9f6abZ16yLLs7J80q/oHBYiIiJV3YABfp1t7dqR5RMmwB13JKdNCVCWPEsNzayTmXUKio4Kvj/S+a1144DbzGyImXUApgA78akAcM6FgEeBsWbWz8w6A08Cy4DXgzqfAq8Ck80sK9gJNxmY5ZxbWVT7qkSepY0b/ZDl119Hlnft6pN9ReeuEBERSRUDB8K//x0bMN1zD4wenZw2VbJSpw4wsz7A/Di3nnDOXRakD/gzcA3QBFgMXO+cWx72HHWBMcBQoB4wD7jOOZcdVudgYDwwKCiaCdzgnPuBqnyQ7tatcMop8NFHkeWdOsEbb0CT6NlHERGRFPTMM3DxxbHTbxMnwnXXJadNJVeq1AHlyrOURFUzWNq+3acCWLQosrxdO5+romnTxLRDREQkEf71L7jqqtjyqVN98sqqK6l5lmqu3bth8ODYQOmoo+D11xUoiYhI9XPllXD//bHll13mky1XE9UuWErKmqU9e/wBg/PmRZa3aOEDpSOOSFxbREREEul3v4O77oos27cPLrwQ5sxJTpsqmKbhyis/H375S5/NNNwhh/ipt+OOq7zXFhERqQqcg1tugbFjI8vr1fOH8vbsmZx2HZim4RLGOb+ILTpQatzYR9MKlEREpCYwgzFjfEbvcLt2wVln+UPkU5iCpfIYORImTYosq1fPn6PTpUty2iQiIpIMZvDQQ/6Q3XA5OXDGGf5EixSlYKms7r8/Np9EerpPB1/1hhtFREQqX1oaTJkC55wTWb55s98t/u23SWlWeVW7YCkhC7ynTYObb44sS0uD6dP9D4OIiEhNVaeO/3vYv39k+Tff+CNTtmxJTrvKQQu8S2v2bBg0CPbujSx//HG/VVJERER87sG+fWHJksjyHj3gtdeSfT6qFnhXmkWL4LzzYgOl0aMVKImIiIRr2NCv4W3bNrJ84UI4/3yfdidFKFgqqRUr/Ir+Xbsiy2+6yW+XFBERkUhNm/rd4dH5Bv/zH5/Qct++5LSrlCo8WDKzO83MRV3fhd23oM46M9tlZgvM7Pio58gwswlmttnMdpjZTDNrWZLXr5Q1S9nZfp5169bI8ksu8VslrVSjeSIiIjXHj37kA6bos1GnToURI5LTplKq8DVLZnYn8HMgfGVXvnNuU3B/BHA7cBmwCvgj0Ato65zbFtR5GBgY1PkeGAscDHRxzuWTyDVL33/vd7d99llk+Zlnwssv+4VsIiIiUrSFC/2i7+gZmr//PRkzNFVizdJe59x3YVdBoGTAb4F7nHMvOueWA5cC9YGhQZ1M4ArgZufc6865/wGXAB2JDMAq344dfuotOlDKyoLnnlOgJCIiUlI9evi/nWlpkeW33urTDVRhlRUstQmm2daY2XQzOzooPwpoDswtqOicywXeBHoERV2AOlF11gHLw+rEyM3NJScnJ+Iqlz174Oc/h8WLI8vbt4dZs6BBg/I9v4iISE1z1lnw2GOx5VdeCf/9b+LbU0KVESwtBn4JnA5chQ+OFprZIcF/A2yIesyGsHvNgTzn3NYi6sQYNWoUmZmZ+69WrVqV/R04B9dcA6++GlnesqWfdz3kkLI/t4iISE32y1/CffdFluXn+x1yVfRYlAoPlpxzs51zLzjnljnnXgfOCm5dGl4t6mEWpyxakXVGjhxJKBTaf2VnZ5e67fvdfbfPmxTu4IP9YYDlCcJERETEJ3aOXqdUsPRlzZrktKkIlZ46wDm3A1gGtAEKdsVFjxA1o3C06Tsg3cyaFFEnRkZGBo0bN464yuTxx+HOOyPL6tXzU2/t25ftOUVERCTS6NFw8cWRZRs2+A1U33+fnDYdQKUHS2aWAbQH1gNr8MHQqWH304HewMKgaCmwJ6pOC6BDWJ3KMXdu7InJtWrBM89A9+6V+tIiIiI1iplfv9S3b2T5ypX+bLndu5PTrjgqI8/SfWbW28yOMrOTgOeBxsATzucpGAfcZmZDzKwDMAXYCTwN4JwLAY8CY82sn5l1Bp7Ej069XtzrlznP0ocfxs/OPX587IGAIiIiUn7p6fDii9ChQ2T5u+/CsGFVJmllZeRZmo7Pm3QosAl4D7jDObciuG/An4FrgCb4BeHXB2kECp6jLjAGn06gHjAPuM45V7AQqWLzLGVn+3QA69ZFlt9yi8//ICIiIpUnO9vP4Hz7bWT5734H999fGa9YqjxLOkj3hx980slPPoksv+ACePppPw0nIiIilevjj+HkkyE69c8DD8Bvf1vRr1YlklKmhtxcOPfc2ECpVy+fIEuBkoiISGKccIKfkqtdO7L8ppvghReS06ZAzY0GnIMrroD58yPL27WDl16CunWT0y4REZGaql+/2KSVzvldc+++m5w2UQ2DpRIv8L7jDnjqqciy5s1h9myfU0lEREQSb9gw+OtfI8tyc/1mqy+/TEqTauaapSlT4PLLI8saNIC33oITT6zAZoqIiEipOQfXXguTJkWWt20LixZBk+hUjKWmNUtFevPN2FxKaWn+cD8FSiIiIslnBhMnwoABkeUrV/o0P3l5CW1OzQqWPv/cL+jesyey/KGHfMZQERERqRpq14bp06Fjx8jy+fPhuuv86FOCVLtg6YBrlrZs8WfObNkSWX7zzbEjTSIiIpJ8jRr548aaR52S9uijMGZMwppRM9Ys1a0Lp58OCxZEVho0yG9TTEurpGaKiIhIuS1ZAr17w65dkeUvvOBnjEpPa5YiOAfDh8cGSp06+d1wCpRERESqtq5d4cknY8svucQHUpWs+gdL//hHbM6Gww+HV16Bhg2T0yYREREpnXPPhdGjI8t27fKzRN98U6kvXe2m4SZOnMjEiRPJz89n1apVhPCn+O5Xv75PEdClS6U3UkRERCqQc3DVVX7NUriOHeGdd6C4s2AL6Ww4gJw33ySzT5/IYMnMz28OGVKpjRMREZFKkpfnd7C/8UZk+YABMHNmSZfXaM0Sa9f6g3Cj/e1vCpRERERSWXo6PP+8T1AZ7r//hVtvrZSXTLmRJTOzUCi0L7o8NzeX3Nxc2LGDBueey44VK2gFZBOMLF1yCTz4oB9dEhERkdS2ejX07Qtbt0aWjx8Pl15a5EMzMzMzgW2uhEFQKgZLjYFQstshIiIiKS3TOZdTkoqpGCwVPbKUn0/GvfeSO2kSrYCvWremyfz5FXY4bteuXVlSCdsUK+N5K/o5c3JyaNWqFdnZ2bHn7ZVDKrz3ynpe9an6NBWeU32aGj+jkBrvv8Kf86mnfEZvYF+zZtSaPr3YTVylHVmqXe5GJliJ3tgjj7Dh6KPhD3+gztSpNG7dusJePy0trcJ/uCvreSurrY0bN67Q502l964+VZ9WpFR5zgLq04pV0f0JqfP+K/Q5hw8n97PP+GT8eNosWECj9u2LfUhJR5TCH1AtL+AI/K65Iyr4ea+vpPZW+PNW9HPil385oHFVbmcq/X9Sn6pPU+Q51acp0J+p8v5TrU+dc6k3DVdSYWubMl1pI0iJS31a8dSnFU99WvHUpxVL/VnxKrtPq2fqAC8XuCv4KhVDfVrx1KcVT31a8dSnFUv9WfEqtU+r7ciSiIiISEWoziNLIiIiIuWmYElERESkCAqWRERERIqgYElERESkCAqWRERERIqgYEkimNlIM1tiZtvMbKOZzTCztlF1zMzuNLN1ZrbLzBaY2fHJanOqCfrYmdm4sDL1aSmZ2RFm9qSZfW9mO83sQzPrEnZffVoKZlbbzP5qZmuC/lptZn8ys1phddSnRTCzXmb2StA/zswGR90vtv/MLMPMJpjZZjPbYWYzzaxlQt9IFVJUn5pZHTMbbWbLgr5aZ2ZTzezwqOcod58qWJJovYGJQBZwKv5InLlm1iCszq3ATcANQFfgO+A1M2uU4LamHDPrClwNfBx1S31aCmbWBHgX2AOcCRwH3Az8EFZNfVo6I4Br8f3VHt9/twC/DqujPi1aA+AjfP/EU5L+GwcMAS4EegINgVlmllZJba7qiurT+sCJwF+Cr+cCxwIzo+qNo7x9WhlpwXVVnwtoik8h3yv43oD1wIiwOhn4P1LXJLu9VfkK/oGuAvoDC4Bx6tMy9+XfgLeLuK8+LX2fzgIejSp7AZimPi1TfzpgcNj3xfYfkAnkAReE1TkcyAdOT/Z7SvYV3acHqNM1qHdkRfapRpakOJnB1y3B16OA5sDcggrOuVzgTaBHYpuWciYC/3HOvR5Vrj4tvUHAB2b2XDBd/D8zuyrsvvq09N4B+pnZsQBm9hP8p/D/BvfVp+VTkv7rAtSJqrMOWI76uKQy8cHSD8H3FdKntSuufVLdmJkB9wPvOOeWB8XNg68boqpvAH6UqLalGjO7ED9M3DXObfVp6R0NDMf/fN4LdAPGm1muc24q6tOyGI3/Q/OZmeUDacDtzrlngvvq0/IpSf81B/Kcc1vj1GmOFMnM6uJHnZ92hefDVUifKliSojwInID/dBkt+pwci1MmgJm1Av4BnOac211EVfVpydUCPnDO3RZ8/79goexwYGpYPfVpyV0AXAIMBT4BOgHjzGydc+6JsHrq0/IpS/+pj4thZnWA6fjfDdeV5CGUok81DSdxmdkE/FTHKc65tWG3vgu+RkfkzYj9xCReF3z/LDWzvWa2F7+Q/jfBfxf0m/q05NYDK6LKPgWODP5bP6elNwb4m3NuunNumXNuGvAAMDK4rz4tn5L033dAerCB4UB1JEoQKP0bP9V5atioElRQnypYkgjB1tYH8bsK+jrn1kRVWYP/4Ts17DHp+D/+CxPW0NQyD+iI/6RecH0APBX892rUp6X1LtA2quxY4Ovgv/VzWnr1gX1RZfkU/p1Qn5ZPSfpvKX6HZ3idFkAH1MdxhQVKbYD+zrnvo6pUSJ9qGk6iTcQPw58DbDOzgk9BIefcLudcQX6g28zsc+Bz4DZgJ/B0Mhpc1TnntuEXE+5nZjuA7wvWgqlPS+0BYKGZ3Yb/RdkNn5LhagD9nJbJK8DtZvYNfhquM36b+2OgPi0JM2sI/Dis6Cgz6wRscc59U1z/OedCZvYoMNbMvsdvrLkPWAZEbwypEYrqU2Ad8Dx+PejZQFrY36wtzrm8CuvTZG8F1FW1LvwcbrzrsrA6BtyJnwrZjd/N0SHZbU+li7DUAerTMvfh2cEvvN34Kbirou6rT0vXn43w+Wi+BnYBXwJ/BdLVpyXuwz4H+P05paT9B9QFJgDf4wOpV4BWyX5vVbFPgdZF/M3qU5F9asETiYiIiEgcWrMkIiIiUgQFSyIiIiJFULAkIiIiUgQFSyIiIiJFULAkIiIiUgQFSyIiIiJFULAkIiIiUgQFSyIiIiJFULAkIiIiUgQFSyIiIiJFULAkIiIiUoT/B84TjPEw7YmQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "Graphics object consisting of 1 graphics primitive"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y(x) = x*(120 - x)\n",
    "plot(y(x), (x, 0, 120), aspect_ratio=1/100, thickness=3, color='red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want a function in time $t$.\n",
    "The assumption of constant horizontal speed implies that $x$ is just a rescaling of $t$.\n",
    "This means that when $t = 45$, x must be 120."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x at 0 : 0\n",
      "x at 45 : 120\n"
     ]
    }
   ],
   "source": [
    "x(t) = 120/45*t\n",
    "print('x at 0 :', x(0))\n",
    "print('x at 45 :', x(45))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then the altitude is given as the composition of $y$ after $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "altitude at 0 : 0\n",
      "altitude at 22.5 : 3600.00000000000\n",
      "altitude at 45 : 0\n"
     ]
    }
   ],
   "source": [
    "f(t) = y(x(t))\n",
    "print('altitude at 0 :', f(0))\n",
    "print('altitude at 22.5 :', f(22.5))\n",
    "print('altitude at 45 :', f(45))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we plot the altitude as a function of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAD9CAYAAABObSg6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4VklEQVR4nO3dd3xUVfrH8c8jQkQwoaxSFFddKzZcDMUGoqyNqrsuq6LYC5ZFRA3qT7GhLioq0V1FEbEhK6KiCAqKhW5BEAULq1GqggkCBgjn98e5gWEITMpM7pTv+/W6rzH3nrnzXG9Inpx7znPMOYeIiIiIbNsOYQcgIiIikuyUMImIiIjEoIRJREREJAYlTCIiIiIxKGESERERiUEJk4iIiEgMSphEREREYlDCJCIiIhJD2iZM5mWbmYUdi4iIiKS2HcMOoJJilicvLCwkJyeHwsLC6ohHREREUlfMzpW07WESERERiRclTCIiIiIxpOojORHJBCUlsHIl/PLL5u3nn2HFCvjtN1izZstt7VrYuBFKFxWPXFw8Kwtq14add978uvPOUK8eNGgADRtuuWVng4ZAikgg7RKm/Px88vPzKSkpCTsUEYll3Tr49lv4+mv4/nv44QcoKNi8LVrkk6Zou+zit9Kkp3TbaSeoUcMnOqXJjplPnFav9snW2rWbk6s1a3xCVly89WfUqQN77LHl1qwZ7Lsv7L8/7L477KBOepFMYc7FHD+9ubHZ5cDlwF7Bri+A251z44LjTwPnRb1tunOuTcQ5soBBwD+A2sBE4Arn3I8RbeoDDwNdgl2vAVc5534Nvo4ZdFFR0aZB39nZ2eW+RhFJgDVrYM4cmD0bvvoK5s/328KFvkcIfA9QaVJSuu2xBzRqBH/4w+aenwYNoFat+MXmnI9vxYote7F++gl+/HHzVprAlcZbuzbst59PnvbfHw49FA4/3P93jRrxi09EqkPM7uSKJkydgRLgm2DXeUA/4Ajn3BdBwtQIOD/ibeuccysizvEY0BnoBfwC3A80AFo650qCNuOAPYBLgrc9DvzPOdc5+FoJk0iyKiyEGTPgk0/gs8/8tmCBTzRq1IA//QkOOMAnFqWv++8PjRsn/yOw9evhf//z1xO5ffWVT6bAJ1KHHAItWvitVSufSNWsGWLgIhJDfBOmMk9gtgLo55x7MkiY6jnnum2jbQ6wHOjpnBsZ7GsKFACnOufGm9lBwDygjXNuetCmDTAVONA5Nx8lTCLJwTn/OG3qVJgyxW9ffOH31627OWko3Q4+2D82S0e//AKff+4TxNmz/fbFFz7Jql0bjjwS2rbdvDVqFHbEIrJZzISp0mOYzKwG8DegDj6ZKdXezJYBvwKTgZucc8uCYy2BmsCE0sbOuUVmNhc4ChgPtAUKS5OloM00MysM2syvbMwiEgfffQcTJ8I778C778Ly5X7/wQf7RODaa/3r/vtn1hifhg3h+OP9Vqq42Pe0TZ3qt+eeg/vu88cOPBBOOAE6dPDvqV8/nLhFpFwqnDCZ2aH4BGkn4Degu3NuXnB4HDAK+B7YG7gDmGRmLZ1zxUBj/CO6lVGnXRocI3hdxtaWRbTZSnFxMcURAzeLiooqeGUiUqaVK2HCBJ8gTZzoxx3tsAPk5sLFF8Oxx0Lr1vqFX5asrM09SqUKCuCjj3yy+dZbkJ/vH0X++c8+gTr5ZDjmGD3CE0kylelhmg+0AOoBZwDDzaydc25e6WO2wFwzm4VPnk4DRm/nnMaWj9nKeuQW3WYLAwcOZMCAAeW6ABGJ4Ztv4PXX4bXX4IMP/Ey1gw6C006DE0+Edu38dHypuGbNoEcPv4GfHThxIkyaBM8843ugcnJ84tS5M5xyih/oLiKhiscYpneAb51zl27j+NfAUOfcvWbWAT8rrkFkL5OZzQbGOOduNbMLgAecc/WizvMr0Mc5N4wyEqeyepiaNWumMUwi5eEczJoFL7/sk6Qvv/S9Ix06+F/anTr5X/SSWBs3wqef+mT19df947waNeDoo6FbN/jrX3UfRBKjWgZ9TwQKnHO9yjjWEPgJuMQ590zEoO9znHMvBW2aAD+y9aDv1s65GUGb1sA0NOhbJH6c8wOTR46El17yY5N23dUnSJ07+56kunXDjjKz/fQTvPGGT54mTPB1q44+Gs480ydPTZuGHaFIuoh7WYG78eOUCoBdgB7AjcDJ+HFNtwEvA4vxtZruBvYEDnLOrQrO8RjQCV9WYAW+JlNDti4r0BQo7bV6HPheZQVE4mD+fD/4eORIPyW+QQM44wz4+9/9o7Yd066ebXooLPS9fyNH+uRpwwY/fqz08Z7GkIlURdwTpieBE4AmQCHwOXCvc+5tM6sNjAGOwI9vWgy8C9zinCuIOMdOwL+As9iycGVkmwZsXbjyShWuFKmkoiLfizRsmJ/6n5MD3bv7JOmEEzTAONWsXAljxvjk6Z13/GO7bt2gVy/o2FFJr0jFJf6RXEi2GXTk0igLFixQwiSZa+NGmDzZJ0n//S/8/jv85S9w/vnQtWv61kPKNEuW+B7DYcN83aemTaFnT588HXhg2NGJpIrMS5hKqYdJMtaKFf6X52OP+XXa9t3XJ0nnnuuXGpH05Bx8/DE8/TQ8/7zvherQAa64Arp0US+iyPYpYVLCJBnjk098TZ8XXvDjW848Ey691Nf0SfYlRyS+fv8dRo+GRx/1NZ+aNoVLLvF1szRQXKQsSpiUMElaW7cORo2CIUNg2jTfg3TZZXDRRVp6Q7zZs31v47PP+kTq9NN9NfY2bWK/VyRzxEyYMmjdApE0UlQEgwbBPvvAOedAnTq+R2HhQrjpJiVLstnhh8O//+1LFAwe7BOotm19z+Mrr/iipCISkxImkVSyaBHccIMvXti/v58RNWeOnynVvbtmR8m25eTAlVf6oqSvvuqXtzn9dD8w/LHHYM2asCMUSWpplzDl5+fTvHlzcnNzww5FJH7mz/cDt/fay/cWXHqp700aNgwOOSTs6CSV7LCDHwT+/vv+Me4RR/hEas894a67fL0nEdmKxjCJJLMvv4Q77/QDuZs0gT59/MDdnJywI5N08t13cP/9MHQo7Lwz/POfcPXVKoYpmURjmERS0rx58I9/wMEH+56A/Hz/S+2665QsSfzts8/m77HzzoN77vG9mTffDL/8EnZ0IklBCZNIMpk3zy9zccghfjr4o4/CN9/A5Zf7xXBFEmn33f3A8IULfU/mgw/6xOmmm+DXX0MOTiRcSphEkkFBAVxwARx6KEyd6gfhfv21LxGgREmqW+PGfhbm//4HvXv7xGmffeBf/4K1a8OOTiQUSphEwrRiBfTrB/vtB2PHwkMP+UTp0kuVKEn4dt3VP5779lvf89m/v/9efeIJXxxVJIOkXcKkWXKSEtas8b+I9tnHz3rLy/O/lK68EmrVCjs6kS01aeIfD3/5JRx3nK8afvDBvmhqak4cEqkwzZITqU4bN/qKy/37w7Jl/pHbzTfDbruFHZlI+X36qf8efustOPpoP+7pyCPDjkqkKjRLTiRpTJvmKyyfd57/JfPll/Dww0qWJPUccQSMGwcTJvjB4Lm5vk7YokVhRyaSMEqYRBLtp5+gZ0+fLK1f78sEjBwJf/pT2JGJVE3HjvDZZ36SwtixsP/+vvilBoZLGlLCJJIoa9f6opP77w/jx/uBsjNnwrHHhh2ZSPzsuKN/tFw6q3PAAL/cyiuvaHyTpJW0S5g06FuSwltv+UGxt98OV1zhf5lcdBHUqBF2ZCKJUa+eL0XwxRe+jtjpp0Pnzr4Ypkga0KBvkXhatMgvKzFqFJxwgq+efMABYUclUr2c8wv8Xn01LF/uC1/266dSGZLM4jvo28wuN7PPzawo2Kaa2SkRx83MbjOzRWa21szeM7ODo86RZWaPmNnPZrbazF4zsz2i2tQ3sxFmVhhsI8ysXkViFalWJSXwyCP+UcTkyfDcc/D220qWJDOZQbduvnL91Vf7x3SHHQbvvBN2ZCKVVtFHcj8CNwJHBtsk4NWIpOh64FrgSiAXWAK8bWa7RJxjMNAd6AEcA9QFxppZ5LOK54EWwMnB1gIYUcFYRarHxx9D69ZwzTVw9tnw1Vdw1ln+l4ZIJqtbF+691w8Mb9zYDxI/6yzf6ySSYqr8SM7MVgD9gKeARcBg59y9wbEsYClwg3PuP2aWAywHejrnRgZtmgIFwKnOufFmdhAwD2jjnJsetGkDTAUOdM7NR4/kJBmsWQO33OJr0Bx6qC9A2aZN2FGJJCfnYMQI6NMHdtjBl9To0UN/WEiySFwdJjOrYWY9gDr4ZGZvoDEwobSNc64YmAwcFexqCdSMarMImBvRpi1QWJosBW2mAYURbbZSXFxMUVHRFptIwnzwARx+uB+jNHAgzJqlZElke8zg3HN9/bEOHXxPU9euvuyGSAqocMJkZoea2W9AMfBvoLtzbh4+WQLfoxRpacSxxsA659zKGG2WlfHRyyLabGXgwIHk5ORs2po1a1buaxIpt9Wr/ZiMdu18wcnZs+H66/3UahGJbbfdfB2yV17xZTaaN/clN1JzApJkkMr0MM3HjylqAzwGDDez5hHHo7/rrYx90aLblNV+u+fJy8ujsLBw01ZQUBDjI0Uq6N13/aO3oUP96u3vv69B3SKVVToo/Iwz/Np0J54I//tf2FGJbFOFEybn3Drn3DfOuVnOuTxgNnANfoA3bN0LtBube52WALXMrH6MNo3K+Ohd2br3apOsrCyys7O32ETi4rfffC2lDh1gzz1hzhw/wFs1lUSqpn59eOopv8TKN9/4mXRPP63eJklK8ShcaUAWsBCf7HTcdMCsFtAOmBLs+hhYH9WmCXBIRJupQI6ZtYpo0xrIiWgjUj2mT/frZg0fDkOGwKRJWtJEJN46dvR/iJxxhl+T7vTTNZNOkk5F6zDdbWbHmtlewVimu4D2wHPOT7cbDPQ3s+5mdgjwNLAGXyYA51wh8CRwv5mdYGZHAM8Cc4B3gjZfAm8BT5hZm2CG3BPA2GCGnEjibdjgq3QffTQ0aOBXZ+/d28/uEZH4y86GYcNg9Gj48ENfLXzs2LCjEtmkoj/9G+HrIc0HJgKtgZOdc28Hx+/DJ02PArOA3YG/OOdWRZyjDzAGeAn4CJ9QdXbOlUS0ORufRE0Its+BnhWMVaRyvvnGr/d2++2+QvGHH/r14EQk8bp3971NrVr5pVUuvhhWrYr9PpEES7ulUfLz88nPz6ekpIQFCxaoDpOUn3Pw5JN+aZPGjeHZZ1UqQCQszvkJFn36QKNG8MILPokSSYyYdZjSLmEqpcKVUiErV8KFF/qpzhdd5GfB1a0bdlQi8u23vmbTJ5/APfdsLnwpEl+JK1wpkjamTfMDu997z4+feOIJJUsiyeJPf/KFYvv0geuug06dYFlZpfpEEksJk2SujRth0CA/XqlpUz+wu3v3sKMSkWi1asF998G4cb6qfosWfsaqSDVSwiSZ6eefoUsX6NcPrr0WJk+GP/4x7KhEZHtOPtlX1z/oIF/o8pZb/IxWkWqghEkyzwcf+L9Qp0+HN9/0q6nXrBl2VCJSHk2a+EKXd97p13E8/nhYvDjsqCQDpF3ClJ+fT/PmzcnNzQ07FEk2zvnkqH17Py7is8/glFPCjkpEKqpGDejf3/cMf/edH4P4/vthRyVpTrPkJDMUFUGvXn4WXP/+MGCAFswVSQdLl0KPHr7n+J57oG9fsJgTnkSiaZacCF9+6eu3TJwIY8bAXXcpWRJJF40awdtv+xl0/frB3/7m/0ASiTMlTJLe/vtfnyztuCPMnAldu4YdkYjE2447+t6lV17xyVNuLnzxRdhRSZpRwiTpacMGuP56/9fmaaf5Wkta3kQkvXXr5ssO1Krl/1AaNSrsiCSNpF3CpEHfwvLl8Je/wAMP+O2FF1SIUiRT7Lef/wOpSxc480xfemDjxrCjkjSgQd+SXubM8T8oV6/2f122axd2RCIShtJZsf37+58JI0bALruEHZUkLw36lgzy2mtw1FGQk+PHKylZEslcZnDjjf7nwqRJ0LatL0EgUklKmCT1lf4l2a0bdOwIH36oqt0i4nXq5B/RFRf7weBaUkUqSQmTpLbff4dzz/V/Sd58s58Vp/FKIhKpeXNf2f/Pf/bjG4cM8X9oiVSAitFI6lqyxPcqzZ7tB3b36BF2RCKSrBo08Iv39usHV10F8+fDgw+qJpuUW9p9p+Tn55Ofn09JSUnYoUgizZnjywVs2OCXRNCsSBGJZccdfZJ04IHQuzcsXAgvvqheaSmXCj2SM7M8M5tpZqvMbJmZjTGzA6LaPG1mLmqbFtUmy8weMbOfzWy1mb1mZntEtalvZiPMrDDYRphZvVgx9u7dm3nz5jFz5syKXJqkkrffhqOPhoYN/eBuJUsiUhGXXgpvvOH/2Dr2WPjpp7AjkhRQ0TFM7YB8oA3QEd9DNcHM6kS1ewtoErGdGnV8MNAd6AEcA9QFxppZjYg2zwMtgJODrQUwooLxSroZNgxOPRWOOcb/sNt997AjEpFUdNJJfoLIzz9D69b+0b7IdlSpDpOZ7QosA9o5594P9j0N1HPOddvGe3KA5UBP59zIYF9ToAA41Tk33swOAuYBbZxz04M2bYCpwIHOua9ixaY6TGnGObj1VrjjDv/X4ZAhGnsgIlW3aBF07gwLFsDIkf4PMslECa/DlBO8roja3z54ZLfAzJ4ws90ijrUEagITSnc45xYBc4Gjgl1tgcLSZCloMw0ojGgjmWLdOjjvPJ8s3XMPPPaYkiURiY+mTX1v9fHH+8Tp0UfDjkiSVKV/65iZAQ8AHzrn5kYcGgeMAr4H9gbuACaZWUvnXDHQGFjnnFsZdcqlwTGC12VlfOyyiDZbKC4upri4eNPXRVqtOj2sXAmnnw5TpmgmnIgkRp06fuHea6/1g8F//BHuussXvxQJVOXP9CHAYfgxSJuUPmYLzDWzWfjk6TRg9HbOZ2y55ElZzwqj22wycOBABgwYUI6wJWX8+KMfZ7B4Mbzzjh+cKSKSCDVqwEMPwZ57wnXX+bIljz+u3mzZpFKP5MzsEaALcLxz7sfttXXOLcYnTPsFu5YAtcysflTT3fC9TKVtGpVxul0j2mwhLy+PwsLCTVtBQUH5LkaS01df+WVOfvvN9y4pWRKR6tC3r193bsQI6N4d1qwJOyJJEhUtK2BmNgQ4HejgnFtYjvc0BJoBi4NdHwPr8bPsSts0AQ4BpgS7pgI5ZtYqok1r/Jip0jZbyMrKIjs7e4tNUtSMGX4WXHa2T5YOPDDsiEQkk5xzDrz+Orz7Lpx4IvzyS9gRSRKo0Cw5M3sUOAvoCsyPOFTonFtrZnWB24CX8QnSXsDdwJ7AQc65VcF5HgM6Ab3wA8YHAQ2Bls65kqDNOKApcGnwGY8D3zvnOrONx3KRNEsuRU2Y4McsHX64/4HVoEHYEYlIppoxwxfI3XVXGD8emjULOyJJnLjPkrsc38vzHj4hKt3+HhwvAQ4FXgUWAMOD17alyVKgDzAGeAn4CFgDdC5NlgJnA3Pws+kmAJ8DPSsYr6SSF17wC2W2b++LUypZEpEwtWoFH30Ea9dC27bwxRdhRyQhqlIdphBtM+jIpVEWLFigHqZU8cgjcM010LMnDB0KNWuGHZGIiLd4MZxyCvzwg1+PrnXrsCOS+IvZw5R2CVMpPZJLEZEFKfv2hfvugx2qWh5MRCTOCgv947nZs/1wgfbtw45I4ivhhStFKs856NNnc0HKQYOULIlIcsrJ8eOY2rb1vU1vvhl2RFLN9NtJwrFxI1x2ma97kp8PN9wQdkQiIttXp47vXTr5ZOjWDUaNCjsiqUZKmKT6bdjglzoZOtQvpnvFFWFHJCJSPllZ8NJLcOaZfuWBYcPCjkiqiUqYSvVatw7OOgtefRWefx7+/vfY7xERSSY1a8Izz0DdunDBBb7A7lVXhR2VJFjaJUyRs+QkyaxdC3/9q1/mZPRov9CliEgq2mEHvxD4LrvA1VfDqlXQv3/YUUkCaZacVI/ffoOuXWHqVN+71LFj7PeIiCQ75+DOO+H//s/P+L3ttrAjksqJOUsu7XqYJAkVFflZJXPm+FkmWhdORNKFGdxyi39Ml5fnE6jbbvP7Ja0oYZLEKiryM0rmzfOP4lq1iv0eEZFUc+ON/jHdDTdASYkvl6KkKa0oYZLEWbXK9yzNm+eXOsnNDTsiEZHEuf56nzT16+dLp9x1l5KmNKKESRJj1SrfszR3rpIlEckc113nk6a+fX3SNHCgkqY0kXYJk2bJJYHoZEmP4UQkk1x7rU+a+vTxSdO99yppSgOaJSfxVfoYbs4cJUsiktkeecSXHLj2Wr/0k5KmZKZZclKNVq2CU0/1ydKECUqWRCSzXXWV72m68krYcUe/ZqaSppSlhEni47fffLI0e7ZPllq3DjsiEZHw9e4N69f7x3O1a6tOUwpTwiRVt3YtdOmyOVlq0ybsiEREksc//+l/Tvbv75MmLTaekpQwSdWsW+eXO5k2zRelVLIkIrK1vDyfNN14I+y0E1xzTdgRSQWlXcKkWXLVaMMGv5DuO+/A2LGq4C0isj0DBsDvv/sep9q14ZJLwo5IKmCHijQ2szwzm2lmq8xsmZmNMbMDotqYmd1mZovMbK2ZvWdmB0e1yTKzR8zsZzNbbWavmdkeUW3qm9kIMysMthFmVi9WjL1792bevHnMnDmzIpcmFbVxo1+le8wYGDVKa8OJiMRi5ksMXHUVXHYZPPNM2BFJBVQoYQLaAflAG6AjvodqgpnViWhzPXAtcCWQCywB3jazXSLaDAa6Az2AY4C6wFgzqxHR5nmgBXBysLUARlQwXkkE5/xAxmef9VuXLmFHJCKSGsxg8GC48EI4/3wYOTLsiKScqlSHycx2BZYB7Zxz75uZAYuAwc65e4M2WcBS4Abn3H/MLAdYDvR0zo0M2jQFCoBTnXPjzewgYB7Qxjk3PWjTBpgKHOic+ypWbKrDlCDO+bL/998PTz7pe5lERKRiSkqgVy944QUYPVp/eIYvZr2HivYwRcsJXlcEr3sDjYEJpQ2cc8XAZOCoYFdLoGZUm0XA3Ig2bYHC0mQpaDMNKIxos4Xi4mKKioq22CQBBgzwydLDDytZEhGprBo1YNgw6NoVzjwTJk8OOyKJodIJU9Cb9ADwoXNubrC7cfC6NKr50ohjjYF1zrmVMdosK+Njl0W02cLAgQPJycnZtDVr1qz8FyPlc//9PmEaONA/gxcRkcrbcUd4/nk45hjfw/Tpp2FHJNtRlR6mIcBhwD/KOBb9nM/K2Bctuk1Z7bd5nry8PAoLCzdtBQUFMT5OKmT4cL+oZF6enxYrIiJVl5UFr7wCBxwAJ50ECxaEHZFsQ6USJjN7BOgCHO+c+zHi0JLgNboXaDc29zotAWqZWf0YbRqV8dG7snXvFQBZWVlkZ2dvsUmcvP66H6B40UVw111hRyMikl522QXefBP+8Ac/4/jHH2O/R6pdRcsKmJkNAU4HOjjnFkY1WYhPdjpGvKcWfnbdlGDXx8D6qDZNgEMi2kwFcsysVUSb1vgxU6VtpDp8+KF/vt6lCzz2mNZBEhFJhD/8wa+UAL6n6Zdfwo1HtlKhWXJm9ihwFtAVmB9xqNA5tzZocwOQB5wPfA30B9oDBzjnVgVtHgM6Ab3wA8YHAQ2Bls65kqDNOKApcGnwGY8D3zvnOhP78Z5mycXD55/DccdBixbw1lu+Oq2IiCTOggV+TNPee8PEiVC3btgRZYqYvQEVTZi21fh859zTQRsDbsUnOvWB6UDviIHhmNlOwL/wyVdtYCJwhXOuIKJNA+Bh/KM/gNeAK51zv6KEKfEWLoSjj4ZGjeC99yAnJ+ZbREQkDj75BNq394uYjx3rxzlJosU3YUoi2ww6cmmUBQsWKGGqjGXLfLLkHHz0kU+aRESk+kye7B/NdekCL74IO1S1CpDEkHkJUyn1MFVSUREcfzwsWuSTpX32CTsiEZHM9OqrcPrpcPXV8OCDYUeT7hJeuFLSSXExdO8O334L48crWRIRCVPXrjBkiF9K5YEHwo4m4+0YdgCSJJzzlbs//BDefhsOOyzsiERE5PLLoaAA+vaFpk2hR4+wI8pYSpjEu+kmX3H2xRf9zDgREUkOd93lazOddx40buwHhEu10yM5gf/8xy938q9/wd//HnY0IiISyQyGDvV/zHbrBnPmhB1RRkq7hCk/P5/mzZuTm5sbdiip4Y034IoroHdv3+UrIiLJp1YtePllX5/plFNUDTwEmiWXyWbNgnbtfCn+l1/2q2eLiEjyWrQI2raF7Gz44AOoVy/siNKFZsnJNixcCJ06wSGH+LFLSpZERJJf06Z+5YWffvKzmtetCzuijKGEKROtWAGnngp16viFdXfeOeyIRESkvA46CF57DaZMgUsu8bOcJeGUMGWa33/3gwaXL4dx42C33cKOSEREKuqYY2DYMBg+3E/akYRLu7ICkUujSBTn4KKLYMYMmDQJ9t8/7IhERKSyzjoLvv7al4XZd18488ywI0prGvSdSe68E265xddaUvkAEZHU5xz07An//a9fKL1Nm7AjSlUa9C2Bl17yydKAAUqWRETShRk8+STk5vqFehcuDDuitKUepkwwY4YvH9C9Ozz3nP8HJiIi6ePnn33vUlaWXzhd5QYqSj1MGe+HH/xfHUccAU89pWRJRCQd/eEPvhDxokV+LNP69WFHlHaUMKWzVaugc2fYaSd45RX/KiIi6emAA2D0aHj3XbjySpUbiLO0S5i0NEqgpATOPts/zx47Fho1CjsiERFJtOOPh8cf99tDD4UdTVqpcMJkZseZ2etmtsjMnJl1izr+dLA/cpsW1SbLzB4xs5/NbLWZvWZme0S1qW9mI8ysMNhGmFm9WPH17t2befPmMXPmzIpeWnq58UbfPTtypK/mLSIimeH886FfP78+6IQJYUeTNirTw1QHmA1cuZ02bwFNIrZTo44PBroDPYBjgLrAWDOLXJ/jeaAFcHKwtQBGVCLezDN0KAwaBIMH+0UaRUQkswwcCCed5GdFL1gQdjRpoUqz5MzMAd2dc2Mi9j0N1HPOddvGe3KA5UBP59zIYF9ToAA41Tk33swOAuYBbZxz04M2bYCpwIHOua9ixZaxs+Q+/BA6dPAFKvPzNchbRCRTFRZC69b+v6dPh5yccONJbqHNkmtvZsvMbIGZPWFmketvtARqApv6CZ1zi4C5wFHBrrZAYWmyFLSZBhRGtJFoBQVwxhlw1FH+2bWSJRGRzJWT49ecW7oU/vEPP7ZVKi0RCdM44GygA9AXyAUmmVlWcLwxsM45tzLqfUuDY6VtlpVx7mURbbZQXFxMUVHRFltGWbvWrxG3004wahTUrBl2RCIiErb99/djWceP92NbpdLinjA550Y6595wzs11zr0OnALsD5wW463GlgUpy3pWGN1mk4EDB5KTk7Npa9asWWXCT03OwcUXw5dfwquvwq67hh2RiIgki7/8BR54wI9tfeaZsKNJWQkvK+CcWwx8D+wX7FoC1DKz+lFNd8P3MpW2KWse/K4RbbaQl5dHYWHhpq2goKDqwaeK++/3FbyHDYMWLcKORkREks3VV8MFF/g/rqdNi91etpLwhMnMGgLNgMXBro+B9UDHiDZNgEOAKcGuqUCOmbWKaNMayIlos4WsrCyys7O32DLC+PFwww2Ql6c14kREpGxm8OijcOSRfpmsn34KO6KUU+FZcmZWF9g3+PJT4FrgXWBFsN0GvIxPkPYC7gb2BA5yzq0KzvEY0AnoFbxnENAQaOmcKwnajAOaApcGn/U48L1zrjNaS877+mto1QqOPto/iqtRI/Z7REQkcy1d6hfqbdwY3n9fK0BslpBZckfiE6VPg68fCP77dqAEOBR4FVgADA9e25YmS4E+wBjgJeAjYA3QuTRZCpwNzMHPppsAfA70rES86amoCLp29RW8n3tOyZKIiMTWqJFfPuXzz/1jOim3KtVhCtE2g87Pzyc/P5+SkhIWLFiQnj1MGzf6GXGTJ8OMGX79IBERkfIaNsyPaXr8cT+uSWL2MKVdwlQqrR/JDRjgt7Fj4dToIuoiIiLlcMUV8OST/tFcaYHLzKWEKe0SpjffhE6d4Pbb4eabw45GRERS1bp10K6dL3r88ceZvki7Eqa0Spi++w5atoRjjvGDvHdI+CRHERFJZ4sWwZ//DAceCG+/nclFj0NbGkXibe1av+xJw4YwYoSSJRERqbqmTf3qEB99BNdfH3Y0SU2/dVOBc3D55TB/Prz8MtSrF3ZEIiKSLo491lcCHzwYnn8+7GiSVtolTPn5+TRv3pzc3NywQ4mfxx+H4cP96+GHhx2NiIikmyuvhJ494aKLfMkB2YrGMCW76dN99n/JJTBkSNjRiIhIulqzxhdCLiqCWbOgfvQKZmlNY5hS2vLl8Ne/+oHeDzwQdjQiIpLOdt7ZF7VcuRJ69fLDQWQTJUzJasMG6NHDT/scNQpq1Qo7IhERSXd77w3PPAOvvQaDBoUdTVJRwpSsbrkF3nsPXnwR9tgj7GhERCRTdOq0eVH3Dz4IO5qkoTFMyWjsWOjcGe69V9M8RUSk+m3YACec4Bd5//TTTChqmXljmFJ+ltz338O55/qEqV+/sKMREZFMtOOO/glHSQmcdZZ/zXDqYUom69bBccfBkiXwySfQoEHYEYmISCabNAk6doSbbvJLcqWvzOthSml5eT5RGjlSyZKIiISvQwefKN15J4wfH3Y0oVLClCxefdWXDrjvPq0aLSIiySMvD04+Gc4+2y/Um6H0SC4ZLFzoFz88/ni/9InF7BkUERGpPr/8AkccAbvvDpMnp2OpGz2SS3rr1sHf/+4rqj71lJIlERFJPg0bwksvwccfw403hh1NKNIuYUq5WXLXXw+zZ/tvRC2qKyIiyapNGz9s5MEH4fXXw46m2lU4YTKz48zsdTNbZGbOzLpFHTczuy04vtbM3jOzg6PaZJnZI2b2s5mtNrPXzGyPqDb1zWyEmRUG2wgzqxcrvt69ezNv3jxmzpxZ0UurfqNHw0MPwf33w5FHhh2NiIjI9l1zjS9sef758OOPYUdTrSrTw1QHmA1cuY3j1wPXBsdzgSXA22a2S0SbwUB3oAdwDFAXGGtmNSLaPA+0AE4OthbAiErEm5y++w4uuAD+9jfo3TvsaERERGIzg2HDYKed/CDwDRvCjqjaVGnQt5k5oLtzbkzwtQGLgMHOuXuDfVnAUuAG59x/zCwHWA70dM6NDNo0BQqAU51z483sIGAe0MY5Nz1o0waYChzonPsqVmxJPeh73To45hhYscI/D87JCTsiERGR8nv/fT9R6ZZb4Lbbwo4mHqp90PfeQGNgQukO51wxMBk4KtjVEqgZ1WYRMDeiTVugsDRZCtpMAwoj2myhuLiYoqKiLbakdcst8Nlnvt6SkiUREUk1xx0Ht94Kd9zh1z3NAPFOmBoHr0uj9i+NONYYWOecWxmjzbIyzr8sos0WBg4cSE5OzqatWbNmFQ6+WkyY4AfNDRwILVuGHY2IiEjl3HSTT5zOPhuWLw87moRL1Cy56Od8Vsa+aNFtymq/zfPk5eVRWFi4aStIxuJay5b5deJOOgn69Ak7GhERkcqrUQOee84PM+nVC1KzrmO5xTthWhK8RvcC7cbmXqclQC0zqx+jTVlLI+/K1r1XAGRlZZGdnb3FllQ2boTzzvPfUMOHww5pV9FBREQyTdOm/nfam2/C4MFhR5NQ8f6tvRCf7HQs3WFmtYB2wJRg18fA+qg2TYBDItpMBXLMrFVEm9ZATkSb1DJ4MLz1FjzzDDQqKxcUERFJQaeeCn37wg03wKxZYUeTMBWeJWdmdYF9gy8/xZcQeBdY4Zz7wcxuAPKA84Gvgf5Ae+AA59yq4ByPAZ2AXsAKYBDQEGjpnCsJ2owDmgKXBp/1OPC9c64zqbY0yscfQ9u2vn7Fv/4VbiwiIiLxFjn7+5NPIOzfuxUXc5ZcZRKm9vgEKdpw51yvoLTArfhEpz4wHejtnJsbcY6dgH8BZwG1gYnAFc65gog2DYCHgS7BrteAK51zv5JKCdOqVX6duJwcmDIlHdffERER8fUFjzgCunSBESlXNjH+CVOS2GbQ+fn55OfnU1JSwoIFC8JPmHr18gvqfvIJ7LdfeHGIiIgk2nPPwTnn+Nezzgo7morIvISpVFL0MJV+4zzzDPTsGU4MIiIi1enss2HsWL9O6l57hR1NeSlhCi1h+vZb3zXZtWsqdk2KiIhUzq+/QosW0KyZL2pZo0aMNySFaq/0LeDX1jnnHNh1V8jPDzsaERGR6lOvnu8omDIF7rkn7GjiRglTIgwcCDNmwLPPpuJMARERkao59ljo398vnzJjRtjRxIUeycXbjBlw1FG+ZPyAAdX3uSIiIslk/frNpQY+/RTq1g07ou3JvEdy+fn5NG/enNzc3Or/8N9+84PdWraEm2+u/s8XERFJFjVr+slPixfDP/8ZdjRVph6meLr0Uv8Y7rPPVEJAREQE4Kmn4MIL4b//hTPOCDuabcm8HqbQvP46PP44PPigkiUREZFS55/vE6WLL4Yffww7mkpTD1M8LF0Khx4KbdrAq6+CxUxURUREMseKFXDYYXDAAfD228m4AL16mBLOOd/VaAZDhypZEhERidaggS/iPGkSPPJI2NFUStolTNU+6Ps//4E33oAnn4TddquezxQREUk1HTr4RehvvBG+/DLsaCpMj+SqYv58X837vPPgsccS8xkiIiLpYu1avyB93bq+sGXNmmFHVEqP5BJm/XpfzXuPPWDQoLCjERERSX61a/tHc59+CnffHXY0FaKEqbLuvNPf8Oeegzp1wo5GREQkNeTm+lqFd9wBs2aFHU256ZFcZXz8MbRu7W/4bbfF99wiIiLpbv16aNsWVq+GTz7xPU/h0iO5uPv9dzj3XD898qabwo5GREQk9dSs6R/NLVzo15xLAWmXMCV8ltytt8I338Dw4ck0WE1ERCS1NG/uF6sfPBjefTfsaGKK+yM5M7sNuDVq91LnXOPguAXHLwHqA9OB3s65LyLOkQUMAv4B1AYmAlc450pLhIbzSG7KFL+Q4MCBcMMN8TmniIhIptq40ZcbWLgQPv8ccnLCiiS0R3JfAE0itkMjjl0PXAtcCeQCS4C3zWyXiDaDge5AD+AYoC4w1sxqJCje2Fav9uUDWreG664LLQwREZG0scMO8PTTsHJl0i/Qm6iEaYNzbknEthw29S79E7jLOTfaOTcXOA/YGTgraJMDXAj0dc6945z7FDgHn3SdmKB4Y8vLg59+8o/iaoSXt4mIiKSVvfaChx7yidOYMSEHs22JSpj2M7NFZrbQzF40s32C/XsDjYEJpQ2dc8XAZOCoYFdLoGZUm0XA3Ig21au0lPvAgbD//qGEICIikrZ69YIuXeCSS+Dnn8OOpkyJSJimA+cCJwEX4xOkKWbWMPhvgKVR71kacawxsM45t3I7bbZSXFxMUVHRFltcFBXBBRdAu3Zw1VXxOaeIiIhsZuaXGispSdrftXFPmJxz45xzLzvn5jjn3gFOCw6dF9ks6m1Wxr5o220zcOBAcnJyNm3NmjWrcOxl6tsXfvkFhg1LxtWVRURE0kPjxv5pzosvwujRYUezlYRnAM651cAcYD/8AG/YuqdoNzb3Oi0BaplZ/e202UpeXh6FhYWbtoKCgqoHP24cDB0K998Pe+9d9fOJiIjItv3jH9C1K1x+edI9mkt4whSUCDgIWAwsxCdEHSOO1wLaAVOCXR8D66PaNAEOiWizlaysLLKzs7fYqmTlSrjoIjjpJLj44qqdS0RERGIz84vZr18PV18ddjRbiHvCZGaDzKydme1tZq2B/wLZwHDniz4NBvqbWXczOwR4GlgDPA/gnCsEngTuN7MTzOwI4Fl8L9U78Y53m/r08aUEhg71N1BEREQSr0kTePhheOEFeOWVsKPZZMcEnHMP4AXgD8ByYBrQxjn3fXD8PnwxykfZXLjyL865VRHn6ANsAF5ic+HKXs65kgTEu7Vx43z5gKeegj32qJaPFBERkcDZZ8OoUf7R3HHHQcOGYUeUfovv5ufnk5+fT0lJCQsWLKh4pe+iIjj4YF+y/a231LskIiIShsWL/e/i006DZ59N9KfF/GWfdglTqUovjXLZZfDcczB3Lvzxj1WJUURERKpixAi/4P2YMX4weOKEtjRKapo0ydeBuPdeJUsiIiJhO+cc6NTJd2asWBFqKOphKrV6NRx6KDRr5ldNVs0lERGR8C1a5IfKdOrke5wSQz1M5XbTTbBkCTz5pJIlERGRZNG0KQwe7McxvfZaaGEoMwD46CM/hfHOO2HffcOORkRERCKde64f/H3ZZfDrr6GEkHYJU35+Ps2bNyc3N7d8b/j9d7jwQmjVCq65JrHBiYiISMWZwb//Db/9Bv36hRNCxo9hysuDBx6ATz/10xdFREQkOT32GFxxhZ+kdfzx8TyzygpsN2GaNQvatIHbb4f+/eMdo4iIiMTTxo3Qvr2v0fT551C7drzOrEHf27RuHVxwARx2WGjdeyIiIlIBO+wATzwBP/wAAwZU70dX66clk3vvhS+/9Muf1KwZdjQiIiJSHgccAP/3fzBoEHzySbV9bGY+kps/3/cs9e0Ld9+dqBhFREQkEdavhyOPhBo1YMYM2LHKS+Nm3iO5mLPkNm6ESy6BPfeEW26p3uBERESk6mrWhKFDYfZsP3GrGmReD9MTT/iEaeJE6NAhkTGKiIhIIvXtC48+6geA77dfVc6kWXJbJEyLF8NBB8Hpp/uxSyIiIpK6Spc1++MffakBi5n3bEvmPZLbrmuugVq1/EAxERERSW116sDjj8N77/mlzRIocxKm11+HUaPgoYegQYOwoxEREZF4OPFE6NULrrvOL9SbIJnxSM7MV/E+9FB4442qdNmJiIhIslmxwv+eP+ooGD26MmfIvEdyZc6Su/lm/z/z0UeVLImIiKSbBg3gkUfglVdgzJiEfET69zBNnEj2iSfC/fdDnz7VEZuIiIhUN+egc2f47DOYNw+2t47s1jRLrrB5c7Jr14Zp0+JR2EpERESS1fff+0dzF14IDz9ckXdm3iO5rXz1la+9pGRJREQkvf3xj3DHHTBkiK8AHkcp18NkZlZYWLgxen9xcTHFxcWbvl49dy4HnnIKBZdeSvZ991VrjCIiIhKSDRt8YeqSEl9uoBzrxebk5OQAq9x2kqJUTJiygcKw4xAREZG0kuOcK9rWwVRMmMrVw7R48WJatWrFvHnz2H333ct17tzcXGbOnBlq24q0LyoqolmzZhQUFGy9wHAcYtE1Vs+5dY3VF0ciz61rrL44EnluXWP1xZHIcxcVFXFws2Z8Uc5rLE8PU8oN7NnexZRll112Kfc3RI0aNUJvW5n22dnZocedCdeY6HPrGqsvDl1j1c6ta6y+OHSNlW9bRPmvcXs9S6XSf9B3BfTu3Tv0tpVpn6hz6xqr99zJEoeusWrnTpY4dI1VO3eyxKFrrNq54ynlHsmVl5ntDvwI7OGc+ynseBIhYjzXdp+7pjJdY3rQNaYHXWN60DVWTjr3MP0M3BO8pqtiYEDwmq50jelB15gedI3pQddYCWnbwyQiIiISL+ncwyQiIiISF0qYRERERGJQwiQiIiISgxImERERkRiUMKUgM7vNzFzUtiTsuKrCzI4zs9fNbFFwPd2ijltw3YvMbK2ZvWdmB4cUbqWU4xqfLuO+Tgsp3Aozszwzm2lmq8xsmZmNMbMDotqk9H0s5zWm9H0EMLPLzexzMysKtqlmdkrE8ZS+j1Cua0z5+xgp+N51ZjY4Yl/K38dI27jGuN1HJUyp6wugScR2aLjhVFkdYDZw5TaOXw9cGxzPBZYAb5vZLtUTXlzEukaAt9jyvp5aDXHFSzsgH2gDdMSvJDDBzOpEtEn1+1iea4TUvo/ga9jdCBwZbJOAVyN+mab6fYTY1wipfx8BMLNc4BLg86hD6XAfge1eI8TrPjrntKXYBtwGfBZ2HAm8Pgd0i/jagMXADRH7soBfgUvDjjce1xjsexoYE3ZscbzGXYPrPC6N7+MW15iO9zHiulYAF6bjfYy+xnS6j0BdYAFwIvAeMDjYnzb3cVvXGO/7qB6m1LVf0I260MxeNLN9wg4ogfYGGgMTSnc454qBycBRYQWVIO2DRz0LzOwJM9st7ICqICd4XRG8puN9jL7GUmlzH82shpn1wPeQTiUN72MZ11gqHe5jPvCGc+6dqP3pdB+3dY2l4nIfU27xXQFgOnAuPqNuBNwMTDGzg51zv4QaWWI0Dl6XRu1fCvyxmmNJpHHAKOB7/A+zO4BJZtYy+EGWMszMgAeAD51zc4PdaXUft3GNkCb30cwOxScPOwG/Ad2dc/PMrPSXacrfx21dY3A45e9jkAT+Gf+4LVpa/HuMcY0Qx/uohCkFOefGRXw5x8ymAt8C5+F/gKer6LL0Vsa+lOWcGxnx5Vwzm4X/R34aMDqcqCptCHAYcEwZx9LlPpZ5jWl0H+cDLYB6wBnAcDNrF3E8He5jmdfonJuX6vfRzJoBDwF/cc79vp2mKXsfy3ON8byPeiSXBpxzq4E5wH5hx5IgpTMAG0ft342t/zpKG865xfh/2Cl1X83sEaALcLxz7seIQ2lzH7dzjVtJ1fvonFvnnPvGOTfLOZeHn7BwDWl0H7dzjWW1TbX72BJ/Tz42sw1mtgE/aeHq4L9L71Uq38ftXqOZ1Yh+Q1XuoxKmNGBmWcBB+AF86Wgh/od0x9IdZlYL/w9jSlhBJZqZNQSakSL3NZiiPAQ4HejgnFsY1STl72M5rrGs96TUfdwOww8KTvn7uB2l17j1gdS7jxPxs6dbRGyzgOeC//6O1L+P271G51xJ9Buqch/1SC4Fmdkg4HXgB3x2fTOQDQwPM66qMLO6wL4Ru/Y2sxbACufcD0Fdjf5m9jXwNdAfWAM8X92xVtb2rjHYbgNexv9D3gu4G/gZeKU646yCfOAsoCuwysxK/3ItdM6tdc65NLiP273G4B7fRmrfR8zsbvzYjwJgF6AH0B44OU3u43avMR3uo3NuFRA5tg4zWw38UjrmLtXvY6xrjPt9DHs6oLZKTaF8EVgErAN+Cr4ZmocdVxWvqT3+uXn09nRw3IJv/MXA7/iZHIeEHXe8rhGoDYwHlgX39ftgf7Ow467A9ZV1bQ7oFdEmpe9jrGtMh/sYXMeTwP+A4uBa3gE6pst9jHWN6XIfy7jm99hyyn3K38ftXWO876MFJxURERGRbdAYJhEREZEYlDCJiIiIxKCESURERCQGJUwiIiIiMShhEhEREYlBCZOIiIhIDEqYRERERGJQwiQiIiISgxImERERkRiUMImIiIjEoIRJREREJAYlTCIiIiIx/D+CD26LqZ8NAgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "Graphics object consisting of 1 graphics primitive"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(f(t), (t, 0, 45), aspect_ratio=1/200, ticks=5, color='red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. The Newton Operator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Newton operator defined below takes on input a function and returns a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def newton_step(f, x):\n",
    "    \"\"\"\n",
    "    Returns the function that will perform one Newton step\n",
    "    to solve the equation f(x) = 0.\n",
    "    \"\"\"\n",
    "    n(x) = x - f(x)/diff(f,x)\n",
    "    return n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, consider the computation of the square root of 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "x |--> x - 1/2*(x^2 - 2)/x"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sqr(x) = x^2 - 2\n",
    "our_sqrt = newton_step(sqr, x)\n",
    "our_sqrt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3/2"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s1 = our_sqrt(2)\n",
    "s1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17/12"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s2 = our_sqrt(s1)\n",
    "s2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "577/408"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "s3 = our_sqrt(s2)\n",
    "s3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.1234155e-6"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d3 = sqrt(2) - s3\n",
    "d3.n(digits=8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "665857/470832 has error -1.59161572810e-12\n"
     ]
    }
   ],
   "source": [
    "s4 = our_sqrt(s3)\n",
    "d4 = sqrt(2) - s4\n",
    "print(s4, 'has error', d4.n(digits=12))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting from an integer approximation, we obtain a sequence of quadratically converging rational approximations.  The quadratic convergence means that the accuracy doubles in each step."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 10.3",
   "language": "sage",
   "name": "sagemath"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}