path: root/src/Lecture7/slides/.ipynb_checkpoints/X2-StudentsRequests-notebook-checkpoint.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Diffie-Hellman key exchange\n",
    "\n",
    "The following is a simple implementation of the classic [Diffie-Hellman key exchange](https://en.wikipedia.org/wiki/Diffie%E2%80%93Hellman_key_exchange) cryptographic protocol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Public key: p = 75521 and g = 58258 \n",
      "\n",
      "[[ Alice's secret key: a = 22794 ]]\n",
      "[[ Bob's secret key: b = 69773 ]] \n",
      "\n",
      "Alice sends h1 = 31067 to Bob\n",
      "Bob sends h2 = 54398 to Alice \n",
      "\n",
      "Alice computed 30031 using h2 and her secret a\n",
      "Bob computed 30031 using h1 and his secret b\n"
     ]
    }
   ],
   "source": [
    "# Public information:\n",
    "p = Primes()[10^3 + randint(1,10000)] # random prime\n",
    "g = randint(2, p-1) # random integer\n",
    "\n",
    "print(\"Public key: p =\", p, \"and g =\", g, \"\\n\")\n",
    "\n",
    "a = randint(2, p-1) # Only Alice knows this\n",
    "b = randint(2, p-1) # Only Bob knows this\n",
    "\n",
    "print(\"[[ Alice's secret key: a =\", a, \"]]\")\n",
    "print(\"[[ Bob's secret key: b =\", b, \"]]\", \"\\n\")\n",
    "\n",
    "h1 = (g^a) % p # Alice sends this to Bob\n",
    "h2 = (g^b) % p # Bob sends this to Alice\n",
    "\n",
    "print(\"Alice sends h1 =\", h1, \"to Bob\")\n",
    "print(\"Bob sends h2 =\", h2, \"to Alice\", \"\\n\")\n",
    "\n",
    "secret_a = (h2^a) % p # Alice can compute this because she knows a\n",
    "secret_b = (h1^b) % p # Bob can compute this because he knows b\n",
    "\n",
    "print(\"Alice computed\", secret_a, \"using h2 and her secret a\")\n",
    "print(\"Bob computed\", secret_b, \"using h1 and his secret b\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## General Diffie-Hellman\n",
    "\n",
    "The following code is an implementation of a generic Diffie-Hellman key exchange protocol that uses a group $G$ instead of $(\\mathbb Z/p \\mathbb Z)^\\times$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Public key:\n",
      "G = Additive abelian group isomorphic to Z/171 embedded in Abelian group of points on Elliptic Curve defined by y^2 = x^3 + x + 156 over Finite Field of size 157 \n",
      "g = (53 : 90 : 1) \n",
      "\n",
      "[[ Alice's secret key: a = 145 ]]\n",
      "[[ Bob's secret key: b = 65 ]] \n",
      "\n",
      "Alice sends h1 = (150 : 80 : 1) to Bob\n",
      "Bob sends h2 = (4 : 58 : 1) to Alice \n",
      "\n",
      "Alice computed (28 : 28 : 1) using h2 and her secret a\n",
      "Bob computed (28 : 28 : 1) using h1 and his secret b\n"
     ]
    }
   ],
   "source": [
    "def genericDH(G):\n",
    "    if G.cardinality() == 1:\n",
    "        print(\"Group is trivial, can't do anything\")\n",
    "        return\n",
    "    g = G.random_element()\n",
    "    while g == G.identity(): # Make sure g is not trivial\n",
    "        g = G.random_element()\n",
    "        \n",
    "    print(\"Public key:\\nG =\", G, \"\\ng =\", g, \"\\n\")\n",
    "    \n",
    "    a = randint(2, G.exponent()-1) # Only Alice knows this\n",
    "    b = randint(2, G.exponent()-1) # Only Bob knows this\n",
    "\n",
    "    print(\"[[ Alice's secret key: a =\", a, \"]]\")\n",
    "    print(\"[[ Bob's secret key: b =\", b, \"]]\", \"\\n\")\n",
    "    \n",
    "    # \"Ternary operator\", I did not explain this\n",
    "    # https://docs.python.org/3/reference/expressions.html#conditional-expressions\n",
    "    h1 = g^a if G.is_multiplicative() else a*g # Alice sends this to Bob\n",
    "    h2 = g^b if G.is_multiplicative() else b*g # Bob sends this to Alice\n",
    "\n",
    "    print(\"Alice sends h1 =\", h1, \"to Bob\")\n",
    "    print(\"Bob sends h2 =\", h2, \"to Alice\", \"\\n\")\n",
    "    \n",
    "    secret_a = h2^a if G.is_multiplicative() else a*h2 # Alice can compute this because she knows a\n",
    "    secret_b = h1^b if G.is_multiplicative() else b*h1 # Bob can compute this because he knows b\n",
    "\n",
    "    print(\"Alice computed\", secret_a, \"using h2 and her secret a\")\n",
    "    print(\"Bob computed\", secret_b, \"using h1 and his secret b\")\n",
    "    \n",
    "E = EllipticCurve(GF(157), [1,-1])\n",
    "G = E.abelian_group()\n",
    "genericDH(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical methods for differential equations\n",
    "\n",
    "## Euler's method (ODE)\n",
    "\n",
    "In sage you can use [`ode_solver()`](https://doc.sagemath.org/html/en/reference/calculus/sage/calculus/ode.html) to solve any ordinary differential equation by hand, but Euler's method is very simple to implement by hand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAGGCAYAAACNCg6xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5yOdf7H8fcYzFCMs3EYzCLHQowIlYTVQee2ItRurUUl20ltRduGTmtpapeKWp1+5dhBpeQUyjGHEFKImRznhBlzz/X747NjHIa5cc993fd9vZ6Px/24m9M1n7lNeff9fq7PN8pxHEcAAAAoUgm3CwAAAAgXBCcAAAA/EZwAAAD8RHACAADwE8EJAADATwQnAAAAPxGcAAAA/ERwAgAA8BPBCfCT4zhKT08XM2MBwLsIToCfMjIyFBcXp4yMDLdLAQC4hOAEAADgJ4ITAACAnwhOAAAAfiI4AQAA+IngBAAAIs6XX0qzZwf+ugQnoAjJyclq2rSpkpKS3C4FAOCHnBypf3/phRcCf+0oh6E0gF/S09MVFxentLQ0lS9f3u1yAAAnMXasNHiw9P33UvPmgb02K04AACBipKVJTz8t9esX+NAkEZwAAEAEGTVKysqy8FQcCE4AACAibN8u/fOf0pAhUq1axfM9CE4AACAiPPmkVK6c9PDDxfc9ShbfpQEAAIJj1Spp4kRrDC/O+3e4qw7wE3fVAUBochypa1dp2zZpzRqpVKni+16sOAEAgLD26afSV19J06cXb2iS6HECisQATAAIXYcPSw8+KF1+uXTNNcX//diqA/zEVh0AhJ6XX5buu09avlxq2bL4vx8rTgAAICzt2ycNGybdeWdwQpNEcAIAAGHqH/+QDh2SnnkmeN+T4AQAAMLO5s3SmDHSo49KNWoE7/sSnOAJ9erVU1RU1AmPgQMHul0aAOAMPPKIVL26TQkPJsYRwBOWLFkin8935O01a9aoa9euuvnmm12sCgBwJubNkyZPliZNksqWDe735q46eNLgwYP18ccfa+PGjYqKivLra7irDgDcl5cntW0rlSghLV5sz8HEihM8JycnR5MmTdKQIUNOGZqys7OVnZ195O309PRglAcAOIW335aWLZPmzw9+aJLocYIHTZs2Tfv371e/fv1O+XkjRoxQXFzckUdCQkJwCgQAFCojw3qbbr5Z6tjRnRrYqoPndO/eXaVLl9ZHH310ys8rbMUpISGBrToAcMnQodLo0dL69VLduu7UwFYdPOWXX37Rl19+qSlTphT5uTExMYqJiQlCVQCAomzcKL30kvTYY+6FJomtOnjMhAkTVK1aNV111VVulwIAOA1Dhti8pocfdrcOVpzgGXl5eZowYYL69u2rkiX51QeAcDFzpvTxx9KHH0plyrhbCytO8Iwvv/xSW7du1V133eV2KQAAP+XkSIMHS5dfLt1wg9vVsOIED+nWrZu4FwIAwsuYMXa8yuTJkp9j94oVK04AACAk7dwpDR8uDRggNW/udjWG4AQAAELS0KFSTIyFp1BBcAKKkJycrKZNmyopKcntUgDAM779VnrzTenZZ6WKFd2upgADMAE/cVYdAARHXp7Urp2UmystWSJFR7tdUQGawwEAQEh5800LTPPnh1ZoktiqAwAAISQtTXr0Uen22907j+5UCE4AACBkPPGElJUljRrldiWFY6sOAACEhBUrpORk6bnnpNq13a6mcDSHA36iORwAik9ennTxxVJmpgWoUqXcrqhwrDgBAADXvfGGjSCYNy90Q5NEjxMAAHDZ7t3SI49IffpInTq5Xc2pEZyAIjAAEwCK19Chks9nvU2hjh4nwE/0OAFA4C1aZL1Nycl2Jl2oIzgBfiI4AUBg5eZKSUk25PLbb0Nv2GVhaA4HAACueOUV6fvvpcWLwyM0SfQ4AQAAF+zYYcMu77lHatvW7Wr8R3ACAABBd999Umys9OyzbldyetiqAwAAQfXRR9LkydK770qVKrldzemhORzwE83hAHD2MjKkZs3s8emnUlSU2xWdHrbqAABA0DzxhLRnjzWGh1tokghOQJEYgAkAgbFkiTR2rDR8uJSY6HY1Z4atOsBPbNUBwJnLn9kkWYAqGaZd1mFaNgAACCejR0urVtnMpnANTRJbdQAAoJht2SI99ZR0770Fq07hiuAEAACKjePYGXSVK0t//7vb1Zy9MF4sAwAAoe7996XPPpNmzJDKlXO7mrPHihMAACgW+/ZJ998v3XSTdM01blcTGAQnAABQLB5+WDp0SPrXv9yuJHDYqgMAAAH31VfSa69Jr74q1azpdjWBw4oTUAQGYALA6cnMlO6+W7rsMumee9yuJrAYgAn4iQGYAOCf+++Xxo+XVq+W6td3u5rAYqsOAAAEzIIFdqzKiy9GXmiSWHEC/MaKEwCc2sGDUsuWUqVKFqCio92uKPBYcQIAAAExfLj088/StGmRGZokmsPhIb/++qt69+6typUrq2zZsmrZsqWWLVvmdlkAEBGWLJGef14aNkxq0sTtaooPK07whH379qlDhw7q3LmzZs6cqWrVqmnz5s2qUKGC26UBQNjLyZHuuktq0UJ68EG3qyleBCd4wqhRo5SQkKAJEyYceV+9evXcKwgAIsizz0rr19uqU6lSbldTvNiqgyfMmDFDbdq00c0336xq1aqpVatWGj9+vNtlAUDYW7VK+sc/pKFDrTE80nFXHTwhNjZWkjRkyBDdfPPN+u677zR48GD95z//UZ8+fQr9muzsbGVnZx95Oz09XQkJCdxVBwD/k5MjtW0r+XzS0qVSTIzbFRU/ghM8oXTp0mrTpo0WLlx45H333XeflixZokWLFhX6NcOGDdPw4cNPeD/BCQDME09II0dK330ntWrldjXBwVYdPKFGjRpq2rTpMe9r0qSJtm7detKvGTp0qNLS0o48tm3bVtxlAkDY+O47acQI6cknvROaJJrD4REdOnTQhg0bjnnfjz/+qLp16570a2JiYhTjhXVnADhNBw9KfftaYHr0UberCS6CEzzhgQce0MUXX6xnn31Wt9xyi7777juNGzdO48aNc7s0AAg7f/ubtGWLtHx55N9Fdzx6nOAZH3/8sYYOHaqNGzcqMTFRQ4YM0d133+3313PkCgBI8+dLl14qPfdc5M9sKgzBCfATwQmA12Vm2pDLGjWkuXMj91iVU2GrDgAA+OXhh6WUFOnzz70ZmiSCEwAA8MOsWdKrr0rJyVKDBm5X4x626gA/sVUHwKv27ZMuuEBq3NhWm0p4eJiRh390wD/Jyclq2rSpkpKS3C4FAILOcaT+/a2/6Y03vB2aJFacAL+x4gTAi956y2Y2vf++dMstblfjPo/nRgAAcDJbtkiDBkl9+hCa8hGcAADACXJzpTvukCpXlsaOdbua0MFddQAA4AQjR0qLFknz5kl0JxRgxQkAABzju++kYcOkxx6TOnRwu5rQQnM44CeawwF4QWamHd5bsaL0zTfeO4uuKGzVAQCAIx54QNqxQ/rkE0JTYQhOAABAkjR1qvTaa9K4cdJ557ldTWiixwkoAgMwAXjBzp3S3XdL114r/elPblcTuuhxAvxEjxOASJWXJ3XrJq1dK61aJVWt6nZFoYutOgAAPG7UKGn2bOmLLwhNRWGrDgAAD1u4UHriCWnoUOmKK9yuJvSxVQf4ia06AJFm3z6pZUupdm1pzhzuovMHK04AAHiQ41gTeHq69M47hCZ/0eMEAIAH/fvf0pQp0uTJUt26blcTPlhxAgDAY1atskGXAwZIN9zgdjXhhR4nwE/0OAGIBFlZUps2UunS0rffSrGxblcUXtiqA4qQnJys5ORk+Xw+t0sBgLN2333S1q3S0qWEpjPBihPgJ1acAIS7t9+WeveW3nhDuvNOt6sJT/Q4AQDgAT/8IN1zjwWnfv3criZ8seIE+IkVJwDhKjNTattWKlHC+prOOcftisIXPU4AAEQwx5H697e+piVLCE1ni+AEAEAEGzfOepvefVdq0sTtasIfPU4AAESoZcvsLroBA6Rbb3W7mshAjxPgJ3qcAISTffuk1q2lypWlBQukmBi3K4oMbNUBABBhHMfGDezbJ331FaEpkNiqA4qQnJyspk2bKikpye1SAMAvL74oTZ8uvfWWlJjodjWRha06wE9s1QEIB3PnSl26SA8+KI0c6XY1kYfgBPiJ4AQg1G3fbn1NzZtLn38ulaQhJ+DYqgMAIAJkZ0s33mjnz733HqGpuPCyAgAQAe69V/r+e7uDrmpVt6uJXKw4wROGDRumqKioYx7x8fFulwUAATF+vD1efVVq08btaiIbK07wjGbNmunLL7888nZ0dLSL1QBAYHz7rTRokPSXv9gIAhQvghM8o2TJkqwyAYgoqanW19S6tTR6tNvVeANbdfCMjRs3qmbNmkpMTNStt96qn3766ZSfn52drfT09GMeABAqDh+WbrlF8vmkDz+USpd2uyJvIDjBEy666CK99dZb+vzzzzV+/HilpKTo4osv1p49e076NSNGjFBcXNyRR0JCQhArBoBTe+ghaeFC6YMPpJo13a7GO5jjBE/KyspS/fr19fDDD2vIkCGFfk52drays7OPvJ2enq6EhATmOAFw3ZtvSv36SWPHWn8TgoceJ3jSOeeco/PPP18bN2486efExMQohgOeAISYRYuke+6R/vhHaeBAt6vxHrbq4EnZ2dlat26datSo4XYpAOC37dul66+X2raVXnlFiopyuyLvITjBEx588EHNnTtXW7Zs0bfffqubbrpJ6enp6tu3r9ulAYBfDhyQrr1WiomRJk+mGdwtbNXBE7Zv367bbrtNu3fvVtWqVdWuXTstXrxYdevWdbs0ACiS40h33SWtXy99841UrZrbFXkXwQme8N5777ldAgCcsWefld5/38YOtGzpdjXexlYdAAAhbPp06W9/k556yoZdwl0EJwAAQtTq1VKvXhaYnnzS7WogEZyAIiUnJ6tp06ZKSkpyuxQAHpKaKvXsKTVoYHObSvA3dkhgACbgp/T0dMXFxTEAE0CxO3hQ6txZ+uUXO8S3Th23K0I+msMBAAgheXlS377SqlXSvHmEplBDcAIAIIQ88YTdPTd5stSmjdvV4HgEJwAAQsTEiTZ64LnnbEI4Qg+tZgAAhIA5c+wMurvvlh580O1qcDIEJwAAXLZhg3TDDdKll0rJyZxBF8oITgAAuGj3bumqq6QaNaQPPpBKlXK7IpwKPU4AALjk0CHrZUpPt7EDFSq4XRGKwooTUAQGYAIoDnl5Up8+0tKl0owZUmKi2xXBHwzABPzEAEwAgfTAA9KYMTZ24Lrr3K4G/mKrDgCAIHvpJWn0aGsEJzSFF7bqAAAIovffl/76V+nRR6UBA9yuBqeL4AQAQJDMmWN9Tb1726BLhB+CEwAAQbBmjW3Ldeokvf46s5rCFcEJAIBitn271KOHVLeuNYOXLu12RThTBCcAAIrRnj1St25SiRLSp59KcXFuV4SzwV11AAAUk8xMmwq+a5e0YIFUq5bbFeFsseIEFIEBmADORE6OdOON0tq10mefSY0auV0RAoEBmICfGIAJwF8+n9SrlzR1qjRzpnT55W5XhEBhqw4AgAByHOm+++zA3g8+IDRFGoITAAABNGyY9Mor0vjx0g03uF0NAo0eJwAAAmTsWOnpp6URI6Q//cntalAcCE4AAATApEm2RffXv0qPPOJ2NSguBCcAAM7SlClSv37SnXdKzz3HVPBIRnACAOAsfPqpdOut0s03W19TCf5mjWj88QIAcIZmz7YG8CuvlN56S4qOdrsiFDeCE1AEBmACKMw330g9e0qXXiq9/75UqpTbFSEYGIAJ+IkBmADyLV0qdekitWxpAy7LlnW7IgQLK04AAJyG1aul7t2lJk2kjz8mNHkNwQkAAD9t2CB17SrVqWPnz5Ur53ZFCDaCEwAAftiyxbbnKleWvvhCqlDB7YrgBoITPGnEiBGKiorS4MGD3S4FQBj4+Wepc2epTBnpyy+lqlXdrghuITjBc5YsWaJx48bpggsucLsUAGHg55+lyy6TSpa08QM1arhdEdxEcIKnZGZmqlevXho/frwqVqzodjkAQtyWLTZuoGRJac4cKSHB7YrgNoITPGXgwIG66qqrdMUVV7hdCoAQt2WLrTSVLm2hqXZttytCKCjpdgFAsLz33ntavny5lixZ4tfnZ2dnKzs7+8jb6enpxVUagBDz00/W01S6tPT114QmFGDFCZ6wbds23X///Zo0aZJiY2P9+poRI0YoLi7uyCOBNXrAE376iZUmnByTw+EJ06ZN0/XXX6/oow6S8vl8ioqKUokSJZSdnX3Mx6TCV5wSEhKYHA5EsM2bbaUpNtZWmmrVcrsihBq26uAJXbp00erVq49535133qnGjRvrkUceOSE0SVJMTIxiYmKCVSIAlxGa4A+CEzyhXLlyat68+THvO+ecc1S5cuUT3g/AezZskK64wuY0zZkj1azpdkUIVfQ4AQA8bdUq6ZJLpPLlCU0RZcYM6YEH7DmA6HEC/JSenq64uDh6nIAI8t130u9/L9WrZ8eoVKnidkUIiGnTpOuvl0qUkPLypOnTpZ49A3JpVpwAAJ40b56dPdekiU0EJzRFiCVLpP797Z/z8qToaFtKDBCCEwDAcz7/3Faa2ra1f+bA3giwZ48Fposusg5/yUKTz2fzJQKE4AQA8JSpU6VrrrHVpk8+kc491+2KcFby8qTXXpMaNZLefVcaPVratMm25+67L6DbdBI9TkCRkpOTlZycLJ/Ppx9//JEeJyCMTZok9esn3Xij/XOpUm5XhLOybJk0YIA1q/XpI40aJcXHF+u3JDgBfqI5HAhv48bZTk6/ftL48baLgzC1d6/0+OPSf/4jNW8uJSdLnToF5VuzVQcAiHgvvCD9+c/SoEG2q0NoClN5edLrr0vnnSe98470z39Ky5cHLTRJBCcAQATLy5MefFB66CFboPjXv+wOdYSh5culiy+W/vQn6corbWrp/fdLJYM7y5tfHwBARDp82LblXnxRGjNGeuYZKSrK7apw2vbutT6mNm2krCxp7lzprbeKvZfpZDhyBQAQcbKypFtukWbNshutbr3V7Ypw2vLypIkTpUcekbKzpZdekgYOdL2jn+AEAIgoe/ZIV18trV5t4wa6dnW7Ipy25cstJC1eLPXqJT3/vFSjhttVSWKrDgAQQbZulTp2tDE+c+YQmsLOvn3WwZ+UJGVk2B/ipEkhE5okVpwAABHihx+k7t2tV/ibb+zGK4SJvDzpzTdtW+7QIbsNctAg17flCsOKE1CE5ORkNW3aVElJSW6XAuAkFi2ylaZKlQhNYWfFCvvDu+suWyJcv1564IGQDE0SAzABvzEAEwhN06ZJt99uuzvTp3PuXNjYv1964gnplVekxo1tiGUAz5QrLqw4AQDC1ujR0g03WDM4h/WGify75c47z56ff15auTIsQpNEcAIAhCGfz85vfeAB6eGHpffek2Jj3a4KRVq50qZ833mndMUVNsRyyJCQ3ZYrDM3hAICwkpUl3Xab9OmndlTZPfe4XRGKtH+/9OSTth3XqJE0e7bUubPbVZ0RghMAIGzs3Cldc40tVHz8sfT737tdEU7JcaT//tfOvDlwQBo1yo5JCaMVpuMRnAAAYWHtWjuizOeTFiyQWrRwuyKc0vff20iBBQtsdPsLL0i1arld1VmjxwkAEPK+/NLOd61Y0YZJE5pCWFqarSpdeKGNcf/qKzv3JgJCk0RwAgCEuAkTpB49LDjNny/Vru12RShU/rZco0bS669LI0daM/jll7tdWUARnIAiMAATcIfPZ4Ok77pL+uMfpY8+ksqVc7sqFGrVKumSS6Q+fWyswPr11tdUurTblQUcAzABPzEAEwie9HQbajlzpvTii7bzExXldlU4QVqa9NRT0ssvSw0b2nOXLm5XVaxoDgcAhJSffrI753791UYOdO/udkU4geNIb78tPfiglJkpjRhh6TYCV5iOx1YdACBkzJkjtW0r5eRYEzihKQStXi1deql0xx32HMHbcoUhOAEAQsK4cXbGa8uW0rff2vFlCCHp6Tblu1Ur6bffpFmzpPff91y3PsEJAOCq3Fzp3nulP/9Z6t/f+poqVXK7KhyRvy3XqJGNav/HP6wZ/Ior3K7MFfQ4AQBcs3u3zUacO1f6978tPCGErFkjDRwozZsn3XyzdeonJLhdlasITgAAVyxfLt1wg53EMWuW3cWOEJGeLg0fLv3rX1L9+tIXX9g+KtiqAwAE33//K3XoIFWpIi1dSmgKGY4jvfOONZj9+9/SM8/Ythyh6QiCE1AEBmACgXP4sN213qePbdHNny/VqeN2VZBkhwF27iz16mWpdt066dFHpZgYtysLKQzABPzEAEzg7KSmSrfcIi1caDtAf/kLQy1DQkZGwbbc734njR0rdevmdlUhix4nAECxW7zYeosPH5a+/lrq2NHtiiDHsXECf/2rtG+f9PTTNm6AFaZTYqsOAFBsHMcWMjp1spuxli0jNIWEH36wo1Fuu01q186GWA4dSmjyA8EJAFAs0tJslWnwYJvTNHeuVKuW21V5XEaGTflu0ULavl367DNp8mQazU4DwQme8Oqrr+qCCy5Q+fLlVb58ebVv314zZ850uywgYn3/vdSmjY0ZmDxZeuklqVQpt6vysPxtucaNpeRk62lavZozbc4AwQmeULt2bY0cOVJLly7V0qVLdfnll+vaa6/V2rVr3S4NiCiOI73+uu3+nHuubc3dcIPbVXncunU25fvWW6WLLrK3H3uMbbkzxF118KxKlSrp+eef1x//+Ee/Pp+76oBTO3BAGjBAevNN6e67rbepTBm3q/KwzExr+P7nP6V69aQxY6QePdyuKuxxVx08x+fz6YMPPlBWVpbat29/0s/Lzs5Wdnb2kbfT09ODUR4QljZskG66Sdq82YJTnz5uV+RhjiN98IHdIbd3rzRsmN05FxvrdmURga06eMbq1at17rnnKiYmRv3799fUqVPVtGnTk37+iBEjFBcXd+SR4PHzmYCTefdd62c6fFj67jtCk6vWr7cp33/4g5SUZHfPPf44oSmA2KqDZ+Tk5Gjr1q3av3+/Jk+erNdee01z5849aXgqbMUpISGBrTrgfzIz7W65iROl22+3EzrKlXO7Ko/KzLTjUV56ye6QGzNGuvJKt6uKSAQneNYVV1yh+vXr6z//+Y9fn0+PE1Bg+XLrNd6xw27S6tOHKeCucBzpww9tW273bmv6fughVpiKEVt18CzHcY5ZUQJQtLw8W9Ro185Wl5Yvl/r2JTS5Yv16Oxrllluk1q1tW+6JJwhNxYzmcHjCY489ph49eighIUEZGRl67733NGfOHH322WdulwaEjdRUqV8/m5n4179K//gHd7S7IivLtuVefNHGsX/8sXTVVW5X5RkEJ3hCamqq7rjjDu3cuVNxcXG64IIL9Nlnn6lr165ulwaEhS++sO04x5FmzpR+/3u3K/Igx7FpokOGSLt2SX/7m/Tww6wwBRk9ToCf6HGCF2Vn2+7P88/brtCbb0rx8W5X5UEbNlgn/qxZUs+e0ujRUmKi21V5Ej1OAIBCrVljg6ZHj5aee85WmghNQZaVZQ3f558vbdokffSRNH06oclFBCcAwDHyG8Bbt5Zyc20200MPSSX4GyN4HEeaMkVq0sT+MB5/XFq7Vrr6arcr8zz+NQCKkJycrKZNmyopKcntUoBit3Wr1KWLNX8PGiQtXSq1bOl2VR7z4492NMqNN0otWtjdck89xfk1IYIeJ8BP9DghkjmO9Pbb0sCBUlycDbW8/HK3q/KYrCzp2WelF16QatWyw/6uucbtqnAcVpwAwON277YTOu64w/qOV60iNAWV40hTp0pNm9qIgaFDbVuO0BSSGEcAAB42darUv7+dM/d//yfdfLPbFXnMxo3SfffZcKwrr5Rmz5bq13e7KpwCK04A4EF79tj5cjfcILVvb200hKYgOnDA5jA1b24TwKdPt0GWhKaQx4oTAHjMtGm2ypSTI02aZAGKI1OCxHEsJA0eLKWkSI8+ag8av8MGK04A4BF790q9e0vXX2/zmdaulXr1IjQFzaZNdjTK9ddbP9OaNdLw4YSmMENwAgAPmDpVatZM+vRT6b//tVWnGjXcrsojDhyQnnzS/gB++MFe/E8+kRo0cLsynAGCEwBEsF9/tQWOG26Q2ra1RY7evVllCor8bbmmTaVRo+xcuR9+kK69lj+AMEZwAorAAEyEo7w86dVX7e/sxYulDz+0hY6aNd2uzCM2b7Yp39ddZ9O/16yR/v53qWxZtyvDWWIAJuAnBmAiXPzwg3TPPdI339jzqFFShQpuV+URBw9KI0fai169ug2xZIUporDiBAARIjtbGjbMjkjZtUuaO1f6z38ITUHz0Ue2xDdypPTgg9K6dbbiRGiKKIwjAIAIsGCBrS5t3Gh3tz/+uBQb63ZVHrF5s3T//dbw3b279MUXUsOGbleFYsKKEwCEsbQ06S9/kTp1sjPmVqywVhpCUxAcPGhLfM2aSatXS1OmSDNnEpoiHCtOABCGHMeOSHngASkjQxo71gJUdLTblXnExx/bUSnbt0sPPSQ99ph0zjluV4UgYMUJAMLM+vVS167SrbdK7dpZM/igQYSmoPjpJzsJ+ZprbGVpzRrpH/8gNHkIK04AECaysqRnnpFefFGqU8eGWfbo4XZVHjBjhjRrlrR7t00SrVZNmjzZBmTR+O05BCcACHGOYzOYBg+WUlPtbNiHH6aPKSj++1+pT5+Ct2+6SZo4kRUmDyM4AUVITk5WcnKyfD6f26XAgzZvlu6913qOr7pK+vpr6Xe/c7uqCHf4sC3nTZxok7/zRUdLCQmEJo9jACbgJwZgIpgyMqQRI6SXXpLi422OYs+e7AwVq1WrLCxNmmSDsC68UGrdWho/3kKTz2dBqmdPtyuFi1hxAoAQkpcnvfWWNHSotH+/9Mgjti3HIkcx2b1beucdC0wrVlj/0h13SH37ShdcYJ9z9dXSnDnSZZcRmsCKE+AvVpxQ3L75xvqYli61O+ZGjbImcATY4cPSZ59JEybYWPCxO0kAAB77SURBVAHHsbvk+vWzbvtSpdyuECGMFScAcNnWrbay9N57tjM0f77UsaPbVUWg1asLtuJ++01q1Up64QXpttukqlXdrg5hguAEAC45cEB67jl7xMVJb7xhO0QlmLAXOLt3S+++a4Fp+XILSL172wvdooXb1SEMEZwAIMgcx/4uf+QRW/gYMsQGT5cr53ZlEeLwYenzzy0szZhhL/jVV0tPPmlbcaVLu10hwhjBCQCCaP58a/ZevNjmJ77wAuMFAmbNmoKtuNRUqWVL6fnnpdtvZysOAUNwAoAgWLfO7pSbPt3ucp89W+rc2e2qIsCePQVbccuWSVWqSL16WaN3y5ZuV4cIRHACgGK0c6c0bJj02mt2h9w770h/+AN9TGclN/fYrbi8PJsO+re/SVdeyVYcihXBCSgCk8NxJjIybBvuhRfsaJQXXpAGDJBiYtyuLIytXWth6b//ta24Cy6wmQ23327zl4AgYI4T4CfmOMEfOTm2ujR8uJSWZnOZHn1UqlDB7crC1N69Nqdh4kRpyRKpcmXbirvzTrbi4ApWnAAgAHJzpbfftm25X36xc2GffpoBlmckN1f64ouCs+J8PtuKmzLFntmKg4sITgBwFvLypA8/lJ56Slq/XrrhBhtG3ayZ25WFoR9+KNiKS0mRzj9fGjnStuKqV3e7OkCSRHsiPGHEiBFKSkpSuXLlVK1aNV133XXasGGD22UhjDmOBaQLL7Rm78REOypl8mRC02nZt0969VXpoovshXv9denmm+0Oue+/lx54gNCEkEJwgifMnTtXAwcO1OLFizVr1izl5uaqW7duysrKcrs0hKGvvpIuvtiON6tQwWYzffqpHZcCP+TmSjNnWuKMj5fuvdfC0eTJ0o4d0pgxlkijotyuFDgBzeHwpF27dqlatWqaO3euLrnkEr++huZwLFokPf649PXXUtu20jPPSFdcwd/vflu3rmArbudOW2G68047AoVVJYQJepzgSWlpaZKkSpUquVwJwsHChdbo/fnndgf89Om22kRg8sO+fdL771tg+vZbqVIl61nq149VJYQlghM8x3EcDRkyRB07dlTz5s1P+nnZ2dnKzs4+8nZ6enowykMImTvXAtPs2bY48t571n7D8Moi+HzSrFkWlqZNs625Hj2si/7qqxlmhbBGcILnDBo0SKtWrdKCBQtO+XkjRozQ8OHDg1QVQoXjWA/T3/8uzZsntWhhrTfXXUdgKtL69QVbcTt2SE2b2n5mr15SjRpuVwcEBD1O8JR7771X06ZN07x585SYmHjKzy1sxSkhIYEepwjlONJnn9kK0+LFUps20pNP2gIJu0mnsH9/wVbc4sVSxYoFW3GtW/PiIeKw4gRPcBxH9957r6ZOnao5c+YUGZokKSYmRjFsKUS8vDwbK/D3v9s4gXbt7Iav7t35O/+kfD7pyy8tLE2dKh0+bFtxH3xgzV/8e4MIRnCCJwwcOFDvvPOOpk+frnLlyiklJUWSFBcXpzJlyrhcHdxw+LD1LI0aZUegdepkbTlduhCYTmrDBunNN6W33pJ+/VVq0sQSZ+/ebMXBM9iqgydEneRvwgkTJqhfv35+XYNxBJEhK8vOknvpJWnrVunKK6VHHpH8nErhPWlpBVtxixbZ4Kr8rbg2bUiZ8BxWnOAJ/P8Bdu+WXn5ZGjvWssBtt0kPP2yneuA4Pp91yOdvxeXk2N7l++9LPXtKsbFuVwi4huAEIKL98ov04ot2kofjSH/6kzRkiFSvntuVhaAffyzYitu+XWrcWBo+3LbiatZ0uzogJBCcAESk776TRo+W/u//pLg46aGHpEGDpCpV3K4sxKSl2Ys0caJN+qxQwZbj+vWTkpLYigOOQ3ACEDFyc21nafRoywC/+531Mv3xj9I557hdXQjx+ezcmAkTpClTbCuuWzfrlr/2WrbigFMgOAFFSE5OVnJysnw+n9ul4CT277eG77FjreH70kttYPXVV0vR0W5XF0I2bizYitu2TWrUSBo2zLbiatVyuzogLHBXHeAn7qoLPZs2SWPGSG+8YYsmt90mDR4stWrldmUhJD29YCvum29s3/LWW+1w3bZt2YoDThMrTgDCiuNIc+ZI//ynDa6sXNmavQcMkOLj3a4uROTl2VbcxIl2XsyhQ7YV9+67thXH7DLgjBGcAISFAwfs7/2xY6Xvv7dDd8ePt5FC5ID/2bTJtuLefNO24s47z86N6d1bql3b7eqAiEBwAhDSNm2SXn3V+pj377eBlS+8wITvIzIyCrbiFiyQype3rbh+/ez8GF4kIKAITgBCjs8nffqp9MordvBupUo2f6l/f7tTzvPy8my/Mn8r7uBBqWtX6Z13pOuuYwkOKEYEJwAhY+dOW1kaP176+Wc70WPCBOkPfyALSJI2by7Yitu6VWrYUHr8cemOO6SEBLerAzyB4ATAVT6f9PnnFpY++kgqXVq65RY73aNtW7erCwEZGdKHH9rq0rx5thX3hz/YVlz79mzFAUFGcALgim3bbIzA66/bP7doYaMFbr/dhld7Wl6eNHeuhaUPP7StuCuukN5+27biypZ1u0LAswhOQBEYgBk4ubnSJ5/Y6tLMmbb9dvvt0t1327ac5xdPfvqpYCvul1+kBg2kxx6zrbg6ddyuDoAYgAn4jQGYZ27LFltZeuMN62NKSrKwdOutUrlyblfnsszMgq24uXPtBcnfirv4YtIkEGJYcQJQLA4csHPjJk6UvvrK8kDv3haYWrZ0uzoXzZghzZ4tVa1qR6B8+KG9WJdfLk2aJF1/PVtxQAgjOAEIGMexUz0mTrTRQhkZ0iWX2ErTzTd7/KDd9HQ7cXj48IL3xcdLjz5qW3F167pXGwC/EZwAnLWtW+3c2DfftIGV9erZMSh9+nh47tKvv9pAyvzHqlXW9J2vRAnbq/zb39yrEcBpIzgBOCNpadKUKba79PXXtrt0003W+H3JJZYLPCMvT1q/3gLS/Pn2/PPP9rEGDaSOHaVBg+wk4gEDpOhom8PQubOrZQM4fQQnAH7Lzra74d5+22Yu5eRIl11mW3E33SSde67bFQZJdra0bFnBatI330h791ogatXKRgZ07Ch16HDiycO1atnU78suk3r2dKN6AGeBu+oAP3n1rrq8PJu7+Pbb1se8f781d/fqZTtNnjg7dv9+aeHCgqD03XcWns45x4ZQduxoj4su8lB6BLyJFScAJ3Aca8l5+23p3Xel7dutb2nAAAtMTZu6XWEx27bt2P6k1avtRale3QLSyJH23KKFVKqU29UCCCKCE1AELw3AXLtW+uADuyNu3TqpcmU7/qRXrwgeKZSXZz/40UFp61b7WKNGFpAeeMCe69eP0BcBgL/YqgP8FKlbdevWWVD6v/+TfvjBjkK79loLTN262dlxEeXQIWnp0oJG7oULbSuuZEnpwgsLtt06dJCqVXO7WgAhhhUnwIM2bCgIS2vW2HDKa6+1Hahu3aSYGLcrDKC9e4/tT1qyxLrazz3XltGGDLGg1LatxwdNAfAHK06An8J5xclxbGVpyhTbilu1ynJDz562stS9uxQb63aVAeA4dsbb0dtua9fax2rUsIDUqZM9n3++rTIBwGkgOAF+CrfglJdniytTp9rjxx9tQeXosFSmjNtVniWfz5bMjg5K27fbx5o0Kdh269hRSkykPwnAWeN/t4AIcviwnRM7dao0bZq0Y4dUpYqFpRdflK64IsxXlg4etFEA+SFp4UI7yqRUKalNG+m22ywkXXyx/eAAEGAEJyDMHTggff65haWPP5b27ZPq1LGz4a6/3nqcw3ZHavduC0f507iXLbN0WL68haOHH7aglJTEwbgAgiJc/3MKeNq2bdInn9jjq69sIaZZM2ngQAtLrVqF4a6U40hbthy77bZunX2sVi3rTerd24JS8+Y2pRsAgoweJ8BPbvY4+Xy2Q/XxxxaWvv/eckPHjtJVV9kJHw0bBrWks5eba13qRwelnTvtY82aFfQmdepkS2hhlwQBRCJWnIAiuDUAMy3NtuA++UT69FPbtapcWerRQxo61Jq7K1QIaklnJyvr2P6kRYukjAwbFJWUJPXpU9CfVKmS29UCQKFYcQL8VNwrTo5jd77lryrNn2+LMhdcYKtKV19tR6GFzQ7Vb7/Z4bf5QWn5cvuBKlSwxqv8FaU2bcK8Yx2Al7DiBLgoLU2aPdtWlr74wlp8YmOlyy+XxoyxwFSnjttV+sFxpM2bC6ZxL1hgKVCyH6BjR6lfP3tu1kwqUcLVcgHgTBGcgCDy+ey0jy++sLC0eLG9r2FDC0ndu1toCvkbxHJzpZUrj+1PSk21PqTzz7e5B8OG2cpSWCQ/APAPwQkoZtu2FQSlL7+0cQHly0tdukjJyXbESWKi21UWITNT+vbbgtWkxYutZykmxo4quesua+Ju3z7MGq8A4PTQ4wTPmDdvnp5//nktW7ZMO3fu1NSpU3Xdddf5/fX+9jjt3y/Nm2dbcF98YXfUlyhh/c/du1tQuuiiEJ+tlJJybH/SihW2NFax4rHTuFu3jrCD7QDg1EL5P91AQGVlZalFixa68847deONNwbwupYtZs+2x/LldtxJnTpS167S8OG2uhSyN4rld6Ufve22aZN9LDHRAtLdd9tz48b0JwHwNIITPKNHjx7q0aPHWV8nO9uONckPSt9+a8Os4+OtP6l/f6lz5xA+Gu3wYVtBOjoo7dplxbZoIf3+9xaSOnSQatd2u1oACCkEJ+AksrOzlZ2drZwcaeXKaM2adUiSlJBg4alSJQtI//ynBabGjUM0KKWnW09SfkhavNhGjcfGSu3aSX/+swWldu2kuDi3qwWAkEZwAo6TmWmzGZ95ZrHmzcuTdJGkspL2SZKefFK68kqbrxSSu1Y7dlh/Un4j9/ff295h5coWkJ5+2hq5W7Wy4ZMAAL8RnOB5u3YVjB+aP7+gD7pKlUt09dW5at/ep/btM5WYaNtvgwbZXXGumzHD9gobNrQAlL+i9NNP9vH69S0oDRhgz40aheiSGACED+6qg6c4jvTLL1Ji4h3q2vXv2ratntavt4/VrWsLMfmP47fe3DyrTjk5Nh3zxx/tMWuWzTfIFxVlK0idOhX0J9WoEdwaAcADCE6IaNafZG09ixbZgsz27faxOnXSdOWVcUeCUkLCqa9V7MHJ55O2bpU2brRwdPTzli223SbZdMyyZaU9eywJligh/eUv0ssvB74mAMAx2KpDRNmxwwJSflBatkw6dMhGDbVo4VOXLrvVqlWWBg9uo8GDn1Dnzp1VqVIlJSQEabq140g7dxYejjZtsqQnSaVK2VbbeedJ111nzw0b2nPNmtJHH0nXXmsH1/l8NhwKAFDsWHFC2MrOttWkRYsKwtLWrfaxOnVsiHW7dvbcsqW0aNEcde7c+YTr9O3bVxMnTizy+53WitOePSeGo/x/zsqyzylRwvYHjw5F+c916hQ9IXPGDGnOHOmyy6SePYusHwBw9ghOCAu5udIPP0hLlthZb0uWSKtW2UiimBipTZtjg1LNmoGv4YTglJl5YijKf967t+ALa9YsPBz97ndM3QaAMMNWHUJOXp5lj/yAtGSJ3el28KD1QDdtakHpzjvtGJOWLYvxrvpDh6TNm62gVavsfT16WM/Rzp0Fn1e5soWhRo2ka64pCEgNGkjnnltMxQEAgo3gBFc5jm2v5QekpUutLyktzT5ev76FoxtvtLB04YXFkENyc6Wffy687+iXX6xIqeAb16hhZ6jkh6OGDUP4PBUAQCCxVYeg8fmkDRts9WjlSntesaJgV6t2bQtHSUn2aN06gHkkL0/69dfC+45++snCk2RbZw0aHLOlNmXNGo3+5BOlSvpx40Z3xhEAAEICwQnF4uBBafXqY0PSqlX2fkmqV8+22Fq1skdSkp31dlYcx6ZZFtZ3tGlTwTePjrZJluedd2LvUULCSceBuzrHCQAQEghOOCv5d9evXm2P77+3kLR+va0wRUfbIMn8gNSqlQWmihXP4pumpZ24pZb/nL/HJ1kIKiwcJSba7f6nieAEACA4wW/p6dKaNfbID0qrVxdstZUtK51//rEhqXlzqUyZM/hmBw7YKlFh4ei33wo+r3r1E+9Wa9jQmqPKlg3Iz52P4AQAIDjhBIcPWy/S0eFo9Wrrk5ZsFalhQwtJRz8SE0/z0Nv8Y0QK6zvKH+8tSXFxBStHxwekIAYYghMAgODkYYcPW1/0unX2yF9JWr/ePiZJtWqdGJAaN5ZiY/38Jj6ftG1b4eHo55/t45ItSxU266hhQ6lKlZA4nJbgBAAgOHlAVpatIOUHpPzHpk0FASkuzrbVzj//2Ge/7mpzHCklpfCm7M2bbcS3VHCMSGHhqGbN01yuCj6CEwCA4BRBdu2y1aLjA1L+MSSSrSA1aWKrRk2aFDyqV/djUWfv3pNPys7MtM+JirJb5o4PR/4eIxLCCE4AAIJTmDl4sKBn+ugMs26dHY8mWQ9S/fonBqTGjU/REjRjhvT113ZmyXnnFd6Unf8NJE8eI0JwAgCE7//+R7CT9Uxv3GjtQvnKly/ILN26FQSkBg0KyS45OXY32sZUKTXVttbyn5cvlxYsOLGQ/GNEzjtPuvrqgnDksWNEkpOTlZycLF9+PxYAwLNYcXLJwYMWjjZvtgbt/OPQNm48sWc6/1SP4xd4qlbMVdTuXceGoJM9H33obL7KlW3qZHq63cXmONZndMstUnIyx4gchxUnAAArTsXEcWyBJz8Y5Yej/Oejz4eNibFb+Rs2lK7v6dP5NXarccUUJZZJVeXDKSqx63+rRFtTpCVHBaLduwvOUctXsaI1LMXH2/MFFxz7dv5ztWoFQyBnzJCuvdb2+Hw+6bbbCE0AABSCFaezkJ5us41++cVWiY4PSAcO2OdFKU+Nq+xRq5qpalYlReeVT1XdmBTFR1kwKpOeqqjU/4WhXbvsXLWjxcUVHn6Of1+1amfeXzRjhjRnjnTZZVLPnmfxqkQuVpwAAASnk3AcW9DJD0YnPH52FLV/r6orVfFKUa3oVDWplKIG56aqTukUVVeqKuak6JzMVJXc+5uiju+POffcE4NQYc/Vqp3h6G0EGsEJAODZrTqfz7bLjg9Duzen6cBPKfLtSFWF7JRjglHr2BTViE5VVV+Kyh/6TdE6fNQFJWWVlcrFS5Xyw0+7wleJqleXzjnHtZ8dAACcmYgMTrm51ga0fbv063ZHv23OUPqPKTr4c6p8O1JU4rdUxexLUVXHQlEDpapTiVRVd1JU2sk55lp5MbFyqlVXiZrxiqpeXYpvffLVIQ/daQYAgBeFXXDKzpZ27JB2/JipvT+kKGNTqg79kirfrymK2pWq2H0pKncgVdVlq0UtlKIyOnTMNXKjS+tgheryVYlXiZrVFVuvhUrXLryHqET58iFx3AcAAHBfyASn7Gwp5acD2vNDqtJ/TNGBLak6vM0apkvuTlFseqrKZ6WoUq6tEiXqwDFfnxtVUullqutQ+erK/V28omo0VUydy5VXv7qcxHhFxRcEopIVKqgcYQgAAJymwAen/AnUnTtLPXsqY9ch7VqTqv3rU5S5OVXZW1OVt8NusY/ZZ83TFbJTVDUvVXWVobpHXSpX0dpXqpoyylTXgbh4+RLPU1aNS7S9TnWd2yBeFRpV1zn1/xeGKlZUpRA/6wwAAIS3wN5V9795QI6kKEmZKqtzj1sZylOU9kRXU1pMdWWVq67sivFyqloPUcxRgahC43iVqFo55A9+ReQ7enL4jz/+yF11AOBhgQ1ODzygvNH/Ugk5ylOUdiR21Lbuf1LZ38Ur7rzqqtwsXufWq6KoktEB+5ZAsDCOAAAQ2K26zp1VYvRoKTpaJXw+1R79oGozTBEAAESIwAannj2l6dOZQA0AACISk8MBP7FVBwCg8xoAAMBPBCcAAAA/EZwAAAD8RI8T4CfHcZSRkaFy5copisnzAOBJBCcAAAA/sVUHAADgJ4ITAACAnwhOAAAAfiI4AQAA+IngBAAA4CeCEwAAgJ8ITgAAAH4iOAEAAPiJ4AQAAOAnghMAAICfSvrzSflndAEAAEQqf84i9Ss4ZWRkKC4uLiBFAQAAhKK0tDSVL1/+lJ/j1yG/p7PilJ6eroSEBG3btq3Ib16UpKQkLVmy5KyuEehrhdJ1AvlaB6qmSL5OKP5uh9p1AnWtUHytA3mtULoO/x0J7nVC8Xc71K4TqGudyWsdsBWnqKio0/4DLl++/Fn/UkRHRwfkX+RAXivUriMF5rWWQu9nC7Xr5Aul3+1Qu06grxVKr3UgrxVq15H470iwrpMvlH63Q+06gb5WoH6384V0c/jAgQND7lqhdp1ACrWfLdSuE0ih9rOF4r9rgRKKP1uoXSeQQu1nC7XrBFKo/Wyh+O9acfBrq+50pKenKy4uzq99QpwdXuvg4vUOHl7r4OG1Di5e7+Aprtc6etiwYcMCdrX8i0ZH67LLLlPJkn7tBOIs8FoHF6938PBaBw+vdXDxegdPcbzWAV9xAgAAiFQh3eMEAAAQSghOAAAAfiI4AQAA+IngBAAA4KeAB6cpU6aoe/fuqlKliqKiorRy5cpAfwtPcRxHw4YNU82aNVWmTBlddtllWrt27Sm/ZuLEiYqKijrhcejQoSBVHXleeeUVJSYmKjY2Vq1bt9b8+fPdLinsnc5rOmfOnEJ/p9evXx/EiiPLvHnzdM0116hmzZqKiorStGnT3C4p7J3ua8rvdfEYMWKEkpKSVK5cOVWrVk3XXXedNmzYELDrBzw4ZWVlqUOHDho5cmSgL+1Jzz33nF566SW9/PLLWrJkieLj49W1a9cij8ApX768du7cecwjNjY2SFVHlvfff1+DBw/W448/rhUrVqhTp07q0aOHtm7d6nZpYetMX9MNGzYc8zvdsGHDIFUcebKystSiRQu9/PLLbpcSMc70NeX3OrDmzp2rgQMHavHixZo1a5Zyc3PVrVs3ZWVlBeYbOMVky5YtjiRnxYoVxfUtIl5eXp4THx/vjBw58sj7Dh065MTFxTn//ve/T/p1EyZMcOLi4oJRoie0bdvW6d+//zHva9y4sfPoo4+6VFH4O93X9Ouvv3YkOfv27QtGeZ4jyZk6darbZUQUf15Tfq+D47fffnMkOXPnzg3I9ehxCmFbtmxRSkqKunXrduR9MTExuvTSS7Vw4cJTfm1mZqbq1q2r2rVr6+qrr9aKFSuKu9yIlJOTo2XLlh3zZyBJ3bp1K/LPAIU7m9e0VatWqlGjhrp06aKvv/66OMsEgobf6+KVlpYmSapUqVJArkdwCmEpKSmSpOrVqx/z/urVqx/5WGEaN26siRMnasaMGXr33XcVGxurDh06aOPGjcVabyTavXu3fD7faf8Z4OTO5DWtUaOGxo0bp8mTJ2vKlClq1KiRunTponnz5gWjZKBY8Htd/BzH0ZAhQ9SxY0c1b948INc8qxnkb7/9tv785z8feXvmzJnq1KnTWRflVce/np988okkKSoq6pjPcxznhPcdrV27dmrXrt2Rtzt06KALL7xQY8eO1ZgxYwJctTec7p8BinY6r2mjRo3UqFGjI2+3b99e27Zt0wsvvKBLLrmkWOsEigu/18Vv0KBBWrVqlRYsWBCwa55VcOrZs6cuuuiiI2/XqlXrrAvysuNfz+zsbEm28lSjRo0j7//tt99O+L/1UylRooSSkpJYcToDVapUUXR09AkrIaf7Z4ACgXpN27Vrp0mTJgW6PMBV/F4Hzr333qsZM2Zo3rx5ql27dsCue1ZbdeXKlVODBg2OPMqUKROoujzp+NezadOmio+P16xZs458Tk5OjubOnauLL77Y7+s6jqOVK1ceE77gn9KlS6t169bH/BlI0qxZs07rzwAFAvWarlixgt9pRBx+r8+e4zgaNGiQpkyZotmzZysxMTGg1w/40cx79+7V1q1btWPHDkk6MjshPj5e8fHxgf52ES0qKkqDBw/Ws88+q4YNG6phw4Z69tlnVbZsWd1+++1HPq9Pnz6qVauWRowYIUkaPny42rVrp4YNGyo9PV1jxozRypUrlZyc7NaPEtaGDBmiO+64Q23atFH79u01btw4bd26Vf3793e7tLBV1Gs6dOhQ/frrr3rrrbckSaNHj1a9evXUrFkz5eTkaNKkSZo8ebImT57s5o8R1jIzM7Vp06Yjb2/ZskUrV65UpUqVVKdOHRcrC19Fvab8XgfHwIED9c4772j69OkqV67ckdXtuLi4wCzwBOTevKNMmDDBkXTC46mnngr0t/KEvLw856mnnnLi4+OdmJgY55JLLnFWr159zOdceumlTt++fY+8PXjwYKdOnTpO6dKlnapVqzrdunVzFi5cGOTKI0tycrJTt25dp3Tp0s6FF14YsNtavexUr2nfvn2dSy+99Mjbo0aNcurXr+/ExsY6FStWdDp27Oh88sknLlQdOfJvhT/+cfR/S3B6inpN+b0OjsL+DCQ5EyZMCMj1o/73TQAAAFAExhEAAAD4ieAEAADgJ4ITAACAnwhOAAAAfiI4AQAA+IngBAAA4CeCEwAAgJ8ITgAAAH4iOAEAAPiJ4AQAAOAnghMAAICfCE4AAAB++n9rbrq56HnpOgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "Graphics object consisting of 2 graphics primitives"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var('y')\n",
    "\n",
    "def euler_desolve(f, x0, y0, x1):\n",
    "    n = 5\n",
    "    h = (x1-x0)/n\n",
    "    S = []\n",
    "    Y = [y0]\n",
    "    for i in range(n+1):\n",
    "        S.append(x0 + i*h)\n",
    "        Y.append(N(  Y[i] + h*f(S[i], Y[i])  ))\n",
    "    return S, Y\n",
    "\n",
    "f(x,y) = y\n",
    "x0 = -1\n",
    "x1 = 2\n",
    "y0 = e^(-1)\n",
    "\n",
    "S, Y = euler_desolve(f, x0, y0, x1)\n",
    "plot(e^x, -1, 2) + line([(S[i], Y[i]) for i in range(len(S))], color='red', marker='o', markersize=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sage also has an `eulers_method()` function \"for pedagogical purposes only\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         x                    y                  h*f(x,y)\n",
      "        -1    0.367879441171442   0.0367879441171442\n",
      "-0.900000000000000    0.404667385288587   0.0404667385288587\n",
      "-0.800000000000000    0.445134123817445   0.0445134123817445\n",
      "-0.700000000000000    0.489647536199190   0.0489647536199190\n",
      "-0.600000000000000    0.538612289819109   0.0538612289819109\n",
      "-0.500000000000000    0.592473518801020   0.0592473518801020\n",
      "-0.400000000000000    0.651720870681122   0.0651720870681122\n",
      "-0.300000000000000    0.716892957749234   0.0716892957749234\n",
      "-0.200000000000000    0.788582253524157   0.0788582253524157\n",
      "-0.100000000000000    0.867440478876573   0.0867440478876573\n",
      "-1.38777878078145e-16    0.954184526764230   0.0954184526764230\n",
      "0.0999999999999999     1.04960297944065    0.104960297944065\n",
      "0.200000000000000     1.15456327738472    0.115456327738472\n",
      "0.300000000000000     1.27001960512319    0.127001960512319\n",
      "0.400000000000000     1.39702156563551    0.139702156563551\n",
      "0.500000000000000     1.53672372219906    0.153672372219906\n",
      "0.600000000000000     1.69039609441897    0.169039609441897\n",
      "0.700000000000000     1.85943570386086    0.185943570386086\n",
      "0.800000000000000     2.04537927424695    0.204537927424695\n",
      "0.900000000000000     2.24991720167165    0.224991720167165\n",
      "1.00000000000000     2.47490892183881    0.247490892183881\n",
      "1.10000000000000     2.72239981402269    0.272239981402269\n",
      "1.20000000000000     2.99463979542496    0.299463979542496\n",
      "1.30000000000000     3.29410377496746    0.329410377496746\n",
      "1.40000000000000     3.62351415246420    0.362351415246420\n",
      "1.50000000000000     3.98586556771062    0.398586556771062\n",
      "1.60000000000000     4.38445212448168    0.438445212448168\n",
      "1.70000000000000     4.82289733692985    0.482289733692985\n",
      "1.80000000000000     5.30518707062284    0.530518707062284\n",
      "1.90000000000000     5.83570577768512    0.583570577768512\n",
      "2.00000000000000     6.41927635545363    0.641927635545363\n"
     ]
    }
   ],
   "source": [
    "# Usage: eulers_method(f, x0, y0, h, x1)\n",
    "eulers_method(f, -1, N(e^(-1)), 0.1, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving the heat equation with a finite difference method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def heat_fdm(u0j, u1j, ui0):\n",
    "    m, n = len(u0j)-1, len(ui0)-1\n",
    "    k, h = 1/m, 1/n\n",
    "    \n",
    "    u = [[0] * (m+1) for i in range(n+1)]\n",
    "    for j in range(m+1):\n",
    "        u[0][j] = u0j[j]\n",
    "    for j in range(m+1):\n",
    "        u[n][j] = u1j[j]\n",
    "    for i in range(n+1):\n",
    "        u[i][0] = ui0[i]\n",
    "        \n",
    "    for j in range(0,m):\n",
    "        for i in range(1,n):\n",
    "            u[i][j+1] = (k/(h*h)) * (u[i+1][j] - 2*u[i][j] + u[i-1][j]) + u[i][j]\n",
    "    \n",
    "    return u\n",
    "\n",
    "n, m = 20, 20\n",
    "u0j = [10 - (j/m)*10 for j in range(m+1)] # One extreme goes from hot to cold\n",
    "u1j = [(j/m)*10 for j in range(m+1)] # The other does the opposite\n",
    "ui0 = [10 - (i/m)*10 for i in range(0,n+1)]\n",
    "\n",
    "u = heat_fdm(u0j, u1j, ui0)\n",
    "for t in range(m+1):\n",
    "    show(line([(i/n, u[i][t]) for i in range(n+1)], ymin=-1, ymax =12))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.0",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}