{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<!-- dom:TITLE: Polynomial interpolation: Error theory -->\n",
    "# Polynomial interpolation: Error theory\n",
    "<!-- dom:AUTHOR: Anne Kværnø (modified by André Massing) -->\n",
    "<!-- Author: -->  \n",
    "**Anne Kværnø (modified by André Massing)**\n",
    "\n",
    "Date: **Jan 21, 2021**\n",
    "\n",
    "If you want to have a nicer theme for your jupyter notebook,\n",
    "download the [cascade stylesheet file tma4125.css](https://www.math.ntnu.no/emner/TMA4125/2020v/part_II/notebooks/tma4125.css)\n",
    "and execute the next cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'> <link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'> <link href=\"https://fonts.googleapis.com/css?family=News+Cycle\" rel=\"stylesheet\"> <link href='http://fonts.googleapis.com/css?family=Nixie+One' rel='stylesheet' type='text/css'> <link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'> <link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'> <link href=\"https://fonts.googleapis.com/css?family=Source+Code+Pro\" rel=\"stylesheet\"> <style> @font-face {\n",
       "    /*font-family: \"Computer Modern\";*/\n",
       "    /*src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');*/\n",
       "    font-family: 'News Cycle', sans-serif;\n",
       "}\n",
       "\n",
       "\n",
       "/* main background */\n",
       "\n",
       "\n",
       "/*#notebook_panel { \n",
       "        background: #C9D4E2;\n",
       "    }*/\n",
       "\n",
       "div.cell {\n",
       "    /*max-width:600px;*/\n",
       "    /*margin-left:16% !important;*/\n",
       "    /*margin-right:auto;*/\n",
       "}\n",
       "\n",
       "div.text_cell {\n",
       "    max-width: 120ex;\n",
       "    /* instead of 100%, */\n",
       "}\n",
       "\n",
       "div.text_cell_render {\n",
       "    font-family: 'Alegreya Sans', Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva,   sans-serif;\n",
       "    line-height: 135%;\n",
       "    /*font-size:    45pt;*/\n",
       "    font-size: 130%;\n",
       "    /*width:600px;*/\n",
       "    /*margin-left:auto;*/\n",
       "    /*margin-right:auto;*/\n",
       "}\n",
       "\n",
       ".CodeMirror {\n",
       "    /*font-family: 'Source Code Pro', monospace;*/\n",
       "    /*font-style: italic;*/\n",
       "    /*font-weight: bold;*/\n",
       "    /*font-size: 120%;*/\n",
       "}\n",
       "\n",
       "\n",
       "/*.prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Nixie One', serifV;\n",
       "    font-style: regular;\n",
       "    font-weight: 400;\n",
       "    font-size: 45pt;\n",
       "    line-height: 100%;\n",
       "        color: #497BB9;    \n",
       "    /*color: #5889C5;*/\n",
       "    /*color: rgb(0, 51, 102);*/\n",
       "    /*font-weight: 200;*/\n",
       "    /*font-size: 50pt;*/\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Nixie One', serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 30pt;\n",
       "    line-height: 100%;\n",
       "    /*color: rgb(0,51,102);*/\n",
       "    color: #497BB9;    \n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.3em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Nixie One', serif;\n",
       "    margin-top:16px;\n",
       "    font-size: 22pt;\n",
       "    font-weight: 600;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: rgb(102,102,0);\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Nixie One', serif;\n",
       "    font-size: 14pt;\n",
       "/*    text-align: center;*/\n",
       "    margin-top: 0em;\n",
       "    margin-bottom: 2em;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Nixie One', sans-serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 16pt;\n",
       "    color: rgb(163,0,0);\n",
       "    /*color: #5889C5;*/\n",
       "    /*margin-bottom: .5em;\n",
       "    margin-top: 0.5em;*/\n",
       "    font-style: italic;\n",
       "    margin-bottom: .1em;\n",
       "    margin-top: 0.8em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".warning {\n",
       "    color: rgb( 240, 20, 20)\n",
       "}\n",
       "\n",
       "</style> \n",
       "<script>\n",
       " MathJax.Hub.Config( {\n",
       "    TeX: {\n",
       "        Macros: {\n",
       "            NN: \"{\\\\bf N}\",\n",
       "            RR: \"{\\\\mathbb{R}}\",\n",
       "            mcO: \"{\\\\mathcal{O}}\"\n",
       "        },\n",
       "        equationNumbers: { autoNumber: \"AMS\" },\n",
       "        extensions: [\"AMSmath.js\", \"begingroup.js\"] \n",
       "    },\n",
       "    tex2jax: {\n",
       "        inlineMath: [ ['$', '$'], [\"\\\\(\", \"\\\\)\"]], displayMath: [ ['$$', '$$'], [\"\\\\[\", \"\\\\]\"]]\n",
       "    },\n",
       "    displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "    \"HTML-CSS\": {\n",
       "        styles: {\n",
       "            '.MathJax_Display': {\n",
       "                \"margin\": 4\n",
       "            }\n",
       "        }\n",
       "    }\n",
       "}\n",
       ");\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from IPython.core.display import HTML\n",
    "from numpy import pi\n",
    "from numpy.linalg import norm, solve  # Solve linear systems and compute norms\n",
    "\n",
    "\n",
    "def css_styling():\n",
    "    try:\n",
    "        with open(\"tma4125.css\") as f:\n",
    "            styles = f.read()\n",
    "            return HTML(styles)\n",
    "    except FileNotFoundError:\n",
    "        pass  # Do nothing\n",
    "\n",
    "\n",
    "# Comment out next line and execute this cell to restore the default notebook style\n",
    "css_styling()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "And of course we want to import the required Modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "\n",
    "newparams = {\n",
    "    \"figure.figsize\": (6.0, 6.0),\n",
    "    \"axes.grid\": True,\n",
    "    \"lines.markersize\": 8,\n",
    "    \"lines.linewidth\": 2,\n",
    "    \"font.size\": 14,\n",
    "}\n",
    "plt.rcParams.update(newparams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Finally, we also run the LagrangePolynomial.ipynb to import\n",
    "the function `cardinal` and `lagrange`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAFpCAYAAACS4uOlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABSpklEQVR4nO3deVhV1f7H8fc6wGE6gIiIgCKCMqgoOKQ2apb1y26zWdktK7PSyrI0y9K8WZaW5a20ObPR2zzfBqfMHDEBFQRRHJgUEJnn9fvjIBcUZObA4ft6nvPU2dP5LsAPm7XXXltprRFCCGG9DJYuQAghROuSoBdCCCsnQS+EEFZOgl4IIaycBL0QQlg5CXohhLBytpYu4HTdunXT/v7+Td4/Pz8fZ2fnliuonets7QVpc2chbW6cyMjIDK21Z23r2l3Q+/v7s2PHjibvv379ekaPHt1yBbVzna29IG3uLKTNjaOUOlTXOum6EUIIKydBL4QQVk6CXgghrJwEvRBCWDkJeiGEsHLtbtSNEK2poqKCjIwMsrOzKS8vt3Q5Debm5kZsbKyly2hT0ub/sbGxoUuXLnTr1g2DofHn5xL0olM5evQoSin8/f2xs7NDKWXpkhokNzcXFxcXS5fRpqTNZlprSktLSU9P5+jRo/j5+TX6uNJ1IzqV/Px8fH19MRqNHSbkReemlMJoNOLr60t+fn6TjmEVQZ+Xl8f8+fPx9PTk4osvxtPTk/nz55OXl2fp0kQ71JQ/fYWwtOb83DZoT6XUhUqp75RSyUoprZSa3IB9wpRSG5RShZX7zVOtcAqVl5fHyJEjWbx4MRkZGWitycjIYPHixYwcOVLCXgjR6TX0V4QJ2A3MAArr21gp5Qr8BqQDwyv3mwXMbFqZdVuyZAmJiYkUFRXVWF5UVERiYiJLlixp6Y8UQogOpUEXY7XWPwE/ASilVjZgl0mAE3C71roQ2K2UCgFmKqWW6hZ8UO3y5cvPCPlTioqKWPr66/S85x5cbW1xs7Ghi60t3ezs8DQacbWxkX5aIYTVa61RN6OAjZUhf8ovwDOAP3Cw+sZKqanAVAAvLy/Wr1/f4A/KzMw86/q8rCxm7Iqn0OnMdXZAV8Cz8uUF+FS+fIHutP+LGHl5eY36elmD5rTZzc2N3Nzcli2oDZSXl3fIuptD2nymoqKiJv3st1bQ9wCOnrYsvdq6GkGvtX4LeAtg2LBhujGzt3l4eJCRkVHn+i7ajS+vgrTBtiSMsiVqhIHowHLSy8vIKy8nvVphp3M0GAhydKS/szODnJ0ZbDIRbjLhbW/f4Ppam8zw1zixsbEdcsheWww1PHHiBCEhIfz1118EBga26medzYQJExg5ciRTp07tkN+r5qjv++zg4EBERESjj9vhx9FPmzaNxYsX19p9Y29rz4QeEzCkgs/OMnx2lnHR62DnZYfHFZ6YrnCn4AInUoxlHC0uJqmoiANFRSQWFpJQWEhqSQlR+flE5efzabXj+hiNDHdx4RxXV851deUcV1ecbGzartFCtILnnnuOK664olVC/t5778XBwQEvLy+++uor9u3bh729PSNHjmTRokUMHDiwatt58+Zx0UUXMXHixE4X9K2ltYI+DXNPSHVe1da1mFmzZvHll1+ecUHWwcGBwMBAlm5ZikO5AyfWnCDrlyyyfs6i+Egxae+nwftpKKOi61h3gq7uhsfV3tj7/+9s/WRZGXEFBezOzyc6L4+ovDz+zssjpaSEbzMz+bay28hWKYaYTFzUpQtjunThfDc3XGw7/O9Q0YkUFBTwzjvv8P3337f4sbXWfPfdd3zyyScsWrSIadOmMXz4cLTWzJs3j0suuYS9e/fStWtXAMLCwggICGD16tU88sgjLV5Pp6S1btQLyAMm17PNfUAO4FBt2RNAMqDOtu/QoUN1Y+Xm5up58+ZpT09PrZTSnp6eet68eTo3N/eMbSsqKnRudK5Oei5JR54bqdepdXodlS+1Tu88f6c+vPSwLjxUWOtnlVdU6Lj8fP1haqq+Pz5eR2zfrg3r1mmqvWzXr9cX7typn01K0jtycnR5RUWj29RQ69ata7Vjt1fNafPevXtbrpA2lJOToysqKvSLL76o+/btq41Go/b19dVz5szRWmtdVFSkZ8yYobt3767t7e31iBEj9MaNG2scY8OGDXrEiBHa2dlZu7q66uHDh+uYmBittdaff/65dnd31xWn/az2799fA7W+5s+frz///HNtNBp1UlJS1T4PPvigDggI0GlpaVprrbdu3aq7du2qS0tLz2hXbm6uNhgM+rvvvquxfMGCBXrkyJHN/8J1MDk5OWddf7afX2CHriNXG3TaqZQyAX0r3xoAP6VUOJCltT6slFoEnKO1Hlu5zSfAfGClUmohEATMARZUFtSiTCYTCxYsYMGCBfX23yqlMIWZMIWZ6P14b0rSS8j8IZOMbzLI+i2Lk3+e5OSfJ0mcmYjrKFe6T+yO542e2Hubz/QNShHs5ESwkxO39ugBQG5ZGZtzcliXnc26EyfYnpvLHydP8sfJk8w9eBBvo5ErPTy4ysODS9zdcZBunnZDWehCtm7CNYYnnniCFStWsHTpUi688EKOHz/O33//DcDs2bP5z3/+w3vvvUdAQABLly7l8ssvJyEhAW9vb8rKyrj66qu56667+PjjjyktLWXnzp3YVP4sbty4kaFDh54xCu3rr78mODiYn376iYiICLTW9O3bl1dffZWbb74ZBwcHnn/+eRYuXMjbb7/Niy++yKeffsqmTZvw8jL/Ef/NN98wfvx4bGv5Kzc3N5eKigrc3d1rLD/nnHNYuHAhhYWFODo6NvprJWpqaP/CMGBdtfcLKl8fAJMBb6CqY09rfVIpdSnwOrADOAG8BCxtfskty+hlxPsub7zv8qYst4zMHzPJ+DKDzB8zydmcQ87mHPY/vJ8uF3fBa5IXntd5YutW88vmYmvLuK5dGVf5p2d2aSlrsrP5JSuLn7OyOFpczNupqbydmorJxobxXbtyg6cnV3h4SN++aJC8vDxefvllXnnlFe68804A+vbty6hRo8jPz2fFihW88847jB8/HoA33niDtWvX8vrrr7Nw4UJycnLIzs7mH//4R1UffEhISNXxDx06hI+Pzxmfm56ejlKKCy64AJPJxP79+ykoKOD888+vCuDnnnuO8ePHExgYyHPPPceaNWvo169f1TG+/fZbnnnmmVrbNWPGDMLDwxk1alSN5T4+PpSWlpKSkmLRC8PWoqHj6NcDdQ4411pPrmVZDHBhUwuzBFsXW7xu8sLrJi/K8srI/CGT46uPk/lTJtlrsslek03CtAQ8rvagx209cB/njsH2zAGYXezsuN7Tk+s9PdFaE5WXx/eZmXybkUFkXh6rjx9n9fHjOBsMXN2tG7d4eTHO3R07uTW/zTXlzNoS4uLiKC4uZuzYsWesS0xMpLS0lPPOO69qmY2NDaNGjWLv3r0AdO3alcmTJ3PZZZcxduxYxo4dyw033FA1QVZhYWHVGXh1UVFRBAQEYDKZANi1axdOTk707du3aptx48YxfPhwnnzySb7//nuGDx9etW7//v0cOHCAyy677Ixjz5w5kz///JM///yz6i+LU079EiksrPf+TNEAkix1sDWZQ3/g1wM5N/1cgt8JpsuYLlQUV3B89XFixsewuedmEmcnkh9X90RDSinCXVx4yt+fHcOGcWDECF4MDGSEiwv5FRV8cuwYV8bE4Lt5Mw8lJLCrk40bFq2relfM+++/z9atW7nwwgv57rvvCA4O5pdffgGgW7dunDhx4oz9o6OjGTRoUNX7Xbt2MXDgwBrzrqxdu5aoqCi01mf8svjmm28YO3Yszs7ONZY//PDDfPrpp6xdu5aAgIAzPjcrKwsAT0/PJrRanE6CvgHsutjhfZc34WvDGZk0kj7P9cEx2JHS9FKOLDnC9tDt7Dx3J6krUykvOPsc530cHXmkVy+2DB3K/hEjWNinD6FOThwvLWVZcjIRkZEM2bGDFcnJnCwra6MWivYuODgYe3t71qxZc8a6wMBAjEYjmzZtqlpWXl7O5s2b6d+/f41tBw8ezGOPPVZ1LeuDDz4AICIioursv7ro6GgGDx5c9X7Xrl013kdFRXHttdfy6quvcs011/D444/X2P/bb7/lmmuuqbFsxowZVSFfvfuout27d+Pj41PrXxmiCeq6SmupV1NG3VTXVqNQKioqdPZf2TpuSpz+w/RH1cidP9z+0Pum79O5MWeO+DnbsbafPKmn79unu27cWDV6x2nDBn1XbKzeeZYr8TLqpnE68qib2bNn6y5duuj33ntP79+/X2/dulUvX75ca631jBkzdI8ePfSPP/6o9+7dq++++27t7OysU1JStNZaHzhwQD/22GN606ZNOikpSa9du1b7+PjoZ555RmutdXR0tDYYDDojI6PqM8vLy7Wzs7P+6quvqpaFhITohQsXaq21TkpK0t7e3nrBggVaa61jYmK0Uqrq+3Ps2DFta2tbNfpGa62nTZumXVxc9Jo1a3RqamrV6/QRcrfffrv+5z//2cJfxfavtUbdWDzYT391lKCvrjS3VKe8m6J3jNjxv6GarNM7L9yp0z9L1+Ul5Q0+VmFZmf4kLU2P/vvvGkM2z42M1J+mpemS8prHkqBvnI4c9OXl5XrRokW6T58+2s7OTvfs2VM/8cQTWuuawyuNRuMZwyvT0tL0tddeq318fLTRaNS9evXSs2bN0iUlJVXbjBw5Ur/22mtV7+Pj4zWg9+/fX7VswoQJumvXrvqjjz7SISEheurUqTXqvPHGG6uGRb777rt61KhRNdZzlqGapxQWFmpXV1f9+++/N/8L18FI0DeQpYMvd1eu3nffvhpn+Zt8NumD/zqoi9OKG3WsuPx8PSM+Xrv98UdV4Pf86y/9wqFDOqvyH6il22sJnTXoW9vPP/+sg4KCdFlZWYsc76qrrtIvvPBCo/d77bXX9KWXXtombW5vWivopY++hZkGmwhaHsSo5FH0e70fTv2dKEkpIWleEpv9NhM7OZa8qIbNkR/s5MQr/fpxdNQoVvTrR4iTE0eLi3nswAH8tmzhkf37Od7K7RGdx+WXX8706dM5evT0aaqa5rzzzuPmm29u9H52dna8+uqrLVKDMFPmXwTtx7Bhw/SOHTuavH97m+RLa0322myOvnqUzO8yzX+oAl0u7kKvR3rR9f+6Nniq5Aqt+SUri6VHj/J75QgJW+C2Hj143M+Pvk61TNFphZo7qVloaGjLFtQG5PmpnUN9bT7bz69SKlJrPay2dXJG38qUUriPdSfsmzBGJIzA90FfbEw2ZK/NJmZ8DNvDtpO6MpWKkop6j2VQiv/z8OC3wYOJHDqUGz09qQDeS0sjeNs2/hkbS1wTnykphLBeEvRtyDHQkX7L+jHyyEgCFgdg9DFSsKeAfXfsY2vgVo68coTy/LMPzzxliIsLqwcMYBVwZ48eGJTio/R0Bmzfzj9jY0koKGjdxgghOgwJeguw62KH3yw/Rh4cScjKEJwGOFF8tJjEhxPZ3HszSc8kUZpd2qBj+QLvhoSQcM453OPtjU1l4Idu28ZdcXEcruPpW0KIzkOC3oIMRgM9bu/B8OjhDPx2IK4jXSnLLCNpXhJbem/h4FMHKc1sWOD7OzryRnAw8eecw52Vk629l5ZG0NatPLJ/P5mlDTuOEML6SNC3A8qg6HZVNyL+imDw2sF0GdOF8pxyDi08xBb/LRyYe6BRgf9uSAix55zDTd27U6w1S48eJWDLFl44fJii8oZ1DQkhrIcEfTuilMJ9jDvha8OJ+DMC98vcKc8r5/Bzh9nSZwsHnjxA6YmGBX4/Jyc+7d+fyKFDGefuTk55OXMOHCBk2zY+TU+nvY22EkK0Hgn6dsrtPDcG/3cwEX9VBn5uOYefPczWgK0cevYQZbkNmwdniIsLvwwezC+DBhHm7Myh4mJuiY3l/L//JlImUBOiU5Cgb+fcRlUG/qYIuozpQll2GQefPMjWwK0c/fdRKGnYccZ17crfw4bxdlAQXnZ2/JWTw/DISKbExXGspIEHEVZtw4YNBAUFUd7A7r2YmBh8fX3JlyG97Z4EfQfhdq4b4WvDGbxmMK4jXSk9Xsr+GfvhNkhblYYur78rxkYppvj4ED9iBLN69cJWKd6tHIO/PDmZcunOsTqlpaU89thjDBo0CGdnZ7y9vbnllls4fPjwGdvOmjWLuXPnnjE3fF3CwsIYOXIkS5e2u+cJidNI0Hcw7he7E/FXBAO/HYjTACdIh7jb49gxdAdZv2Y16BiutrYsDgxkz/DhXN61K9llZUxPSGB4ZCTbcnJauQWiLRUUFLBz507mzp3Lzp07+fbbbzly5AiXX345ZdWmwf7rr7+Ii4vjxhtvbNTx77jjDlasWFHjWKL9kaDvgJQyj9IZHjUcHgP7XvbkR+UTfVk0UZdGkRfdsLl0+jk58VNYGF8NGICfvT1/5+UxcudOHkhIIEf+4bYrV1xxBffeey8zZszA3d0dd3d3Zs2aRUXF2e+odnNz47fffmPixIkEBwdzzjnn8OabbxIbG0tsbGzVdp988gljx46terKT1ppLL72USy65pOrCfV5eHv369WP69OlV+40bN46srCzWW+jZu6JhGvrMWNEOKRsFl8M5888h+dVkDj13iBO/n2BH+A563NmDPs/0qXqoeZ3HUIprPT0Z17UrzyQl8eKRI7yWnMxXx4/zWr9+XGvlT/hZr9Zb5HNH69GN3ufjjz9m8uTJbN68mejoaO6++268vb2ZOXNmo46TU/lXW/UHcm/cuLHG2bxSig8++IBBgwbx4osvMmvWLB588EGMRiMvvvhi1XZGo5Hw8HA2bNjAJZdc0ug2ibYhQW8FbBxt8Jvth/dd3iQ9k0TK6ymkvZvGsc+O0fuJ3vSc2RMbh7P3uzrb2PB8YCC3eHlxT3w8W3JyuG7PHm7w9OTVvn3pYX/2Xxii9Xl7e/Pvf/8bpRQhISHEx8ezdOnSRgV9SUkJjzzyCP/4xz/o2bNn1fLaHg7u4+PDO++8w8SJE8nJyeHjjz9m27ZtVWf91bdLSkpqVttE65KgtyJ2Hnb0e6UfvtN8SZydSOa3mRyce5DUt1MJWBKA5/We9c6UOchkYlNEBMuTk3n84EG+OH6cNSdO8HLfvtzm5dXgmTY7iqacWVvKyJEja3z9R40axVNPPUVOTg6urq717l9WVsatt95KdnY23333XY11hYWFODg4nLHPNddcwy233MLChQtZvHhxjccInuLo6CgP8W7npI/eCjkFORH2TRiDfx+M80BnipKK2DthL1GXRJG/t/6hcAaluL9nz6qLtSfKypgcF8c/YmJIKS5ugxaIllZWVsbNN99MdHQ0a9aswcPDo8b6uh4OXlRUxPbt27GxsWH//v21HjsrK0se4t3OSdBbMfex7gz9eyj9lvfDtqst2Wuz2T5oO/tn7qcsp/6LrX4ODvwUFsbKkBDcbGz4MSuLAdu381FamtxZawFbt26t8XXfsmULPj4+9Z7Nl5aWMnHiRKKjo1m3bh09KudCqq6uh4PPmjWL4uJifvvtN95///0z/hIA84O8hwwZ0oQWibYiQW/lDLYGfO/zZUT8CHzu8wENR18+yrbgbaR/Uv9UCEopbu/Rgz3nnMMVlUMx/xkXx41795IhN1q1qZSUFB566CH27dvHF198wZIlS3j44YfPuk9ZWRkTJkxgy5YtfPrppyilSEtLIy0trUZ3y2WXXcaff/5ZY9+ff/6ZN998k48++ogxY8bw9NNPM2XKFNLS0qq2SUpKIjk5mXHjxrVsY0WLkqDvJOw87AhaHsTQHUNxHeVKSVoJsZNiibo4ivzY+rtzfO3t+SEsjHeDgzHZ2PDF8eOE7djBz5mZbVC9AJg0aRLl5eWMGDGCu+++m7vuuqveoD969CjffvstKSkpDB06FG9v76rX6tWrq7a79dZbiY+PZ8+ePQAcP36cO+64gyeffJIRI0YAMGfOHEJDQ7njjjuqThA+/fRTxo0bR+/evVup1aIlyMXYTsYlwoWIPyNIW5lG4uxEstdns2PwDvwe88PvCT9sHOsenaOU4k5vb8Z06cLtcXFsPHmSK2JieMDXlxcCAnBs4B2VomlsbW157bXXeO211xq8j7+/f4O62dzd3XnwwQdZunQp7777Lp6enjXO3AEMBgMbNmyoel9cXMyKFSv49NNPG94IYRFyRt8JKYPC+05vRuwbgffd3uhSzaGFh9getp2s3+u/u7aPoyPrwsN5ISAAO6V4NTmZ4ZGRROc17EYt0T498cQTBAQENHium0OHDjF37lzOO++8Vq5MNJcEfSdm52FH8FvBRPwZYR6dk1hE9KXRxN4eS0nG2fvfbZRitp8fm4cMIcjRkT0FBZwTGcny5GS5UNuGNm7ciMlkqvPVGK6uro2a6yYoKIh77rmnKWWLNiZdNwK389wYunMoR146wqEFh0hflU7WT1n0XdaX7jd3P+vY+aEuLuwcNoyH9+/n7dRUpick8PuJE7wTHExXO7s2bIV1++mnn3BxcTlj+bBhw9i1a1fbFyQ6FAl6AYDBzkDvOb3xvMGT+HviyV6bTeykWNI/SSdoRRAOvc68meYUZxsb3goO5hJ3d+7et4+vMzKIzM1ldf/+jHRza8NWdD6Ojo707dvX0mWIdk66bkQNTn2dGPz7YILfDca2iy1ZP2axfcB2Ut5MqbdL5sbu3dk1bBjnuLhwuLiYC3bt4uUjR6QrRwgLk6AXZ1DKfLF2+N7hdLu2G+W55cTfG0/UpVEUJp39Vvc+jo5sjIjgoZ49KdOamYmJXLdnDydlNkwhLEaCXtTJ3tueAV8OoP/q/th1syN7TTbbB24n+Y2zX3A1Ggy83LcvXw0YgJuNDd9kZDBMRuUIYTES9OKslFJ0v7E7w/cOx/NGTyryK0i4L4Hoy6MpOlJ01n2v9fQkctgwBjs7s7+wkJE7d/LRaWOzhRCtT4JeNIjR08iA1QPo/5/+2HrYcuLXE2wP225+jOFZzu4DHR3ZPGQIk3v0oLCign/GxTEjIYHSeh6YIYRoORL0olG6T+jO8N3D8bjKg/KT5cTdHseeCXvOOu7e0caG94KDeSMoCDul+HdyMpdERclDyTuhlStXNnp8f23Wr1+PUoqMjIwWqKpxXnzxRfz9/dv8c5tDgl40mn0PewZ+M5Dg94KxcbEh48sMdoTtIPPnuue9UUpxj48PG8LD8TYa+ePkSYZGRhKZm9uGlTdfXl4e8+fPx9PTE4PBgKenJ/Pnzyevla8/3HvvvVx55ZWN2mf06NHcf//9rVRR2/H396/xVCuAc889l9TU1DOmW26vlFJ88cUXFvt8CXrRJEopvO/wZlj0MNwucKMkrYSYK2JIeCCB8sK6b6Ef5eZG5NChnOfqytHiYs7/+28+TU9vw8qbLi8vj5EjR7J48WIyMjLQWpORkcHixYsZOXJkq4e9pZS0w7+8jEYjPXr0sLoH4bQWCXrRLI7+joSvCyfg+QCUnSL5tWQih0WSF1V36Hnb27M2PJwp3t4UVVRwS2wsjx84QHk7H2+/ZMkSEhMTKSqqeRG6qKiIxMRElixZ0ma1TJ48mSuvvJJly5bh6+uLu7s7d9xxBwUFBVXrN2zYwOuvv45SCqVU1eP+9u7dy/jx43FxcaF79+7cfPPNNSYwO3XsF154gZ49e1Y9ctDf35+nn36aW2+9FZPJRI8ePc440z58+DDXXnstLi4uuLi4cN1113H06NE625GYmMjVV19Njx49cHZ2ZsiQIfzwww9V60ePHs2hQ4eYNWtWVTug9q6br776irCwMOzt7enVqxfPPvtsjetH/v7+LFy4kHvuuQdXV1d69uzZoO/Z4sWL6dGjByaTidtuu+2MX+jbt29n3LhxdOvWDVdXV84//3w2b95c43MBJkyYgFKq6n1tbf/555/rracpJOhFsykbhd9jfgzZMgTHYEcK9hYQeU4kR5cdrfNCrdFg4K2gIF7t2xcb4PnDh7lu927y2vF4++XLl58R8qcUFRWxYsWKNq1n48aN7N69m99//53Vq1fz9ddfs2zZMgCWLVvGqFGjuOOOO0hNTSU1NZVevXqRmprKhRdeyMCBA9m2bRu///47eXl5XH311VRUu0C+YcMGoqOj+e9//8uaNWuqli9dupTQ0FB27tzJggULeOKJJ/jqq68AqKio4OqrryY9PZ1169axbt06UlJSuOaaa+r8OcjLy+P//u//+O2334iKiuL666/nuuuuIy4uDjCHd8+ePZk3b15VO2oTGRnJhAkTuO6664iJieH5559n0aJFZ8z0+fLLLxMWFsbOnTt57LHHmD17do1QPt1//vMfnnzySRYsWMDOnTsJDg5m6dKlNbbJzc3ln//8Jxs3bmTbtm2Eh4dzxRVXkFk5hff27dsBePvtt0lNTa16X1vbb7311qq2tyitdbt6DR06VDfHunXrmrV/R9Pe2luWX6b33btPr2OdXsc6HXVFlC5OLz7rPmuysrT7xo2adev0oG3b9KHCwrNu35w27927t8n7KqU0UOfLYDA0+dj1ueWWW/T48eOr3t9+++26Z8+euqysrGrZlClT9NixY6veX3TRRXr69Ok1jvPUU0/piy++uMayrKwsDeitW7dWHbtbt266qKioxna9e/fWl1xySY1ld911lz7vvPO01lr/+uuv2mAw6IMHD1atT0xM1Eop/dtvv2mttX7//fe1s7PzWds6YsQI/cwzz+icnJyqz12yZEmNbdatW6cBffz48aqvz5gxY2psM3/+fO3r61uj/ptuuqnGNn379tXPPPNMnbWMGjVKT5kypcaysWPH6t69e9e5T0VFhe7Ro4f+8MMPq5YB+vPPP69zn1OGDRt21nrO9vML7NB15Kqc0YsWZeNkQ9CKIAZ8NQBbd1uyfspix+AdnFh75vNIT7nY3Z0tlbNgRufnc05kJFtzctqw6oap78JfW18Y7N+/f42ZJn18fDh27NhZ94mMjOSPP/6oMcNlr169AHNXwikDBw7E3t7+jP1HjRp1xvtTjyCMjY3Fx8enxoiUgIAAfHx8an1MIUB+fj6zZ8+mf//+uLu7YzKZ2LFjB4cPHz57408TGxt7xnTJ559/PsnJyeRU+1kaNGhQjW3q+5rFxsbW2ubqjh07xj333ENQUBBubm64uLhw7NixettQW9v//vvvRre9IWRSM9EqPK/1xGW4C7G3xnJyw0miLomi99ze9J7fG4PtmecXQU5ObBkyhAl79rAmO5vRu3bxUWgo17ejh05PmzaNxYsX19p94+DgwH333dem9didNjuoUqpG90ttKioqGD9+/Bl96wBeXl5V/+/s7NwyRVarrTaPPvoo//3vf3nxxRfp168fTk5O3HbbbS16Abj6Zzfla1af22+/nfT0dF5++WX8/f2xt7dn7Nix9bahtrZPmjSpVS5+yxm9aDUOPR0IXxNO7/m9QcGhhYeIujiK4uTiWrd3t7Pj50GDqi7S3rBnD4sPH243k6LNmjWLwMBAHBxqzuTp4OBAYGAgs2bNslBltTMajWc8RGTIkCHs2bOH3r1707dv3xqv2qZBPt2WLVvOeB8aGgpAaGgoKSkpVRd9AQ4cOEBKSgr9+/ev9Xh//vknt912G9dffz2DBg2iZ8+eNf6yqKsdpwsNDWXTpk1nHLtnz54NatfZjltbm0//nAceeIDx48czYMAAXFxczriWYGdnd0Ybamv7wYMHm1zr2UjQi1albBR9nu7D4DWDMfoYObnxJDvCd5D1S+1PsrKrvEj7QkAAAI8dOMC0hATK2sGdtCaTiS1btjB79uwa4+hnz57Nli1bWuRGoJbk7+/Ptm3bSEpKIiMjg4qKCqZPn87JkyeZOHEiW7du5cCBA/z+++9MnTqV3Abc07BlyxYWLVpEQkICb7/9NqtWrap6bu0ll1zCoEGDmDRpEjt27GDHjh1MmjSJIUOGcPHFF9d6vKCgIL7++mt27txJTEwMt9566xl/Mfn7+7Nx40aSk5PrvEHqkUceYcOGDTz99NPEx8fz8ccf89JLLzF79uxGftVqmjFjBh988AFvv/02CQkJLFq0iK1bt57Rho8++oi9e/eyfft2brrpJoxG4xltWLNmDWlpaZw4caLOthcX134S1FwS9KJNuI92Z9jfw3Af505pRinRl0dzYO4BKsrODHBV+fSqz/v3x14p3khJ4do9e8hv4CPuWpPJZGLBggUcO3aM8vJyjh07xoIFC9pdyIO5a8BoNNK/f388PT05fPgwPj4+bNq0CYPBwOWXX86AAQOYPn069vb2tfbJn27mzJlER0cTERHBk08+yb/+9S9uuOEGwPx9+/bbb/H09GTMmDGMGTOGHj168M0339TZdbN06VK6d+/OBRdcwP/93/8xcuRILrjgghrb/Otf/+LIkSMEBgbiWUdX3pAhQ/j888/58ssvGThwIHPmzGHOnDnNvmFs4sSJPP3008ydO5eIiAhiYmKYOXNmjW3ee+898vLyGDp0KDfddBN33nnnGXfOvvTSS6xbt45evXoRERFRZ9tP7/9vMXVdpbXUS0bdNE5Ha29FeYVOWpik1xnMo3L+HvO3Lk6re1TOn9nZumvliJxhO3botOJii426saRTI1AsqbbRL62pPbS5rdXXZhl1IzoEZVD0ntub8LXhGHsYyV6XzY6IHWRvzK51+/Pc3Ng8ZAgBDg7syM3l3J07SW7bkoXo8Boc9EqpaUqpg0qpIqVUpFLqgnq2v0UptUspVaCUSlNKfaSU6tH8koU16HJRF4buHIrbhW6UpJawa8wujiyt/WlUQU5ObB4yhGEuLhwoKuJ+YHs7HH4pRHvVoKBXSk0ElgHPARHAX8DPSim/OrY/D/gQ+AAYAFwD9Ac+bn7JwlrYe9szeM1ges3uBeWQ+Egie2/aS1nemXfHdjcaWTd4MJe5u5MNjNm1i1+yar+gK1pHUlISjz76qKXLEE3Q0DP6mcBKrfXbWutYrfUDQCpQ18DhUcBRrfXLWuuDWustwKvAiOaXLKyJwdZA4AuBDPhyADYuNhz/z3F2jthJwb6CM7Y12dryfVgYlwL5FRX8IyaGzzrIhGhCWFK9Qa+UMgJDgV9PW/UrcG4du20CvJVS/1Bm3YCbgJ+aU6ywXp7XeTJk2xCcQp2q5srJ+P7MoXR2BgNzgEd69qRUa26JjeW1s0yaVZvauoeEaO+a83Or6ttZKeUDJAMXaa3/qLZ8HjBJax1cx37XASsBR8x34P4GXK21PuPp0kqpqcBUAC8vr6GfffZZkxoD5omC2uNQt9Zide0tAF4ATv2kTQb+SY1TklNt/hR4q3LZ7ZWv+iatdXV1xc/Pr0FDCduT8vLyGtMddAbS5pqKi4s5fPhwjSkdqhszZkyk1npYbetaZQoEpVR/zF01zwC/AN7AEuBN4LbTt9dav0Xlv9lhw4bp0aNHN/mz169fT3P272issb36/zSHXzjMwScOwkrolt2NkFUh2LqYf1xPtXk0MCI1lbv37eMDoIuvL0v79sVwljnKc3JySE9Px9fXF0dHxw4zn3lubm6z7vDsiKTNZlprCgsLOXHiBH379sXV1bXRx21I0GcA5YDXacu9gLqe9Pw4sE1rfWqy52ilVD6wUSn1hNa6cX9ri05FKUXvOb0xhZuIvTmWjG8y2DlqJ2HfhuEY6Fhj2zu9velia8tNe/eyLDmZk+XlvB0UhK2h9l7JU/9IUlJSKC0tbfW2tJSioqIzpl6wdtLm/7Gzs8PLy6tJIQ8NCHqtdYlSKhK4FPi82qpLgS/r2M0J8y+H6k69l7H7okE8LvdgyNYh7L56NwV7CogcHsmAzwfAaX/ZXufpyQ9hYVy7ezcr09LILSvjk/79MZ4l7Jv6D8ZS1q9fX3VHZWchbW45DQ3dpcBkpdQUpVSoUmoZ4AO8AaCUWqWUWlVt+++Bq5VS9ymlAiqHW/4b2Km1bvk5OIXVcgpyYsjWIXhc6UHZiTKiLouCr8+8MDWua1d+GzwYNxsbvszI4NrduylsB1MmCNEeNCjotdargYeAJ4FdwPnAFVrrQ5Wb+FW+Tm2/EvOQzPuB3cAXQDxwdcuULToTW1dbBn47EL/H/cx/F/4b4u+Lp6K05jw557q5sS48nG52dvyUlcX4mJh2/cQqIdpKg7tRtNbLtdb+Wmt7rfXQ6iNwtNajtdajT9v+Va31AK21k9baW2s9SfrmRVMpgyLguQBCPwoFO0h9M5Xoy6IpzarZzx7h4sKG8HC8jUbWZWczLjqakxL2opOT/nLRoXhN8oJlVM2Ts3PkTgria95c1d/ZmT/Cw+llb8/mnBwujYriRAe68CpES5OgFx1PKAzZNgTnwc4UJhSyc8TOMx5V2NfJiT/Cw+nj4MD23FzGRkWR0QpP7hGiI5CgFx2SQy8HIv6MwOMqD8qyy4i+LJrU92s+1cff0ZEN4eH0c3Tk77w8Lo6K4riEveiEJOhFh2VrsmXgVwPp9WgvdJlm3537OPDEAXTF/0bk9HJwYEN4OCFOTsTk53NxVBTHJOxFJyNBLzo0ZaMIXBJI0BtBYAOHFx1m7817KS/639BKb3t71g0eTIiTE7vz87l41y4Je9GpSNALq+Bzjw+DfhxUNQNm1CVRlGb+7wJsD3t71oeH09/JiT0FBYyRsBediAS9sBpdL+tKxJ8RGH2N5GzKYeeonRQm/m8OPS+jkXXh4QxwcmJvQQFjpc9edBIS9MKqmAaZGLp16P9G5IzaSc62/832191oZE14OKGV3TiXyGgc0QlI0AurY+9rT8TGCNzHuVN6vJRdY3aR+WNm1Xovo5G1lX320fn5XBodTZaMsxdWTIJeWCVbF1vCfgjD63YvKgoqiLk6hpR3UqrW97C3Z+3gwQQ5OrIrL4/L5A5aYcUk6IXVMtgZCHk/hN5P9oZyiL87nqRnkqomRPO2t2fN4MEEODiwIzeX/4uOJlfCXlghCXph1ZRS9HmmD/1W9AMDJM1LImFaArrcHPY9HRxYGx6OX+V0Cf+IiaFAZr0UVkaCXnQKvvf6MuCLASh7RcobKeyZsKdqrH3vyrD3MRrZcPIk1+3eTXFFRT1HFKLjkKAXnYbntZ4M/m0wtl1syfg6g+jLoyk7ae6qCXR0ZM3gwXja2fHLiRPctHcvpRL2wkpI0ItOpcsFXQj/Ixyjt5GTG06ya/QuStLNwytDnJ35bfBgutja8k1GBrfHxVF+2gNOhOiIJOhFp2MKMxHxVwSOfR3J25XHzvN2UnjQfGPVYJOJ/w4ahMnGhk+PHWNafPwZT7MSoqORoBedkqO/IxGbIjBFmChKLOLv8/8mf08+ACNcXfkhLAwHg4G3UlN57MABCXvRoUnQi07L2N1I+Lpw3C50oySlhL8v/LvqLtqLunThiwEDsFWKJUeOsOiwPOpYdFwS9KJTs3WzZdB/B5kfPp5VRtTYKE6sNz/EZLyHBx+FhqKAuQcPsiI52bLFCtFEEvSi07NxtGHAVwPoPqk75XnlRF8eTcYPGQBM7N6dN4KCAJiekMBn6emWLFWIJpGgFwLzXbShq0LxudcHXazZc+0e0j8zh/pUHx8W9emDBv4ZF8d/MzPPfjAh2hkJeiEqKYOi3/J+9JptfmJV7C2xpL5nfjzhY35+PNKzJ2Vac/2ePWw5edLC1QrRcBL0QlSjlCLg+QD6LOwDGvbdtY/k15NRSrEkMJDJPXpQUFHB+JgYYvPzLV2uEA0iQS/EaZRS9J7bm8ClgQAk3J/A4SWHUUrxdlAQV3p4kFVWxrjoaI4UFVm4WiHqJ0EvRB16PdzLPBkacGD2AZKeScLWYGB1//6c5+rK0eJiLpO57EUHIEEvxFn43utLyMqQqpkvDz51EEeDge/Dwhjg5ERsQQH/iImhUGa8FO2YBL0Q9ehxew9CPwoFGzi08BAH5hygi60tvwweTC97e/7KyeHmvXspk0nQRDslQS9EA3jd7MWA1QNQtooji4+Q+EgiPkYj/x00CHdbW77NzGR6QoJMlSDaJQl6IRrI83pPBnw5AGWnOPryUfY/vJ9QJye+rzYvzsJDhyxdphBnkKAXohG6XdWNAV8NQBkVycuS2f/gfs51deXT0FAMwLykJFamplq6TCFqkKAXopG6XdmNgV8PRNkrkl9LJuGBBK7u1o1/9zOP0Lk7Pp5fsrIsXKUQ/yNBL0QTeFzhwcBvzGGf8noK+x/czzQfHx7r1Ysyrblhzx7+zs21dJlCABL0QjSZx+WVYW80n9nvf3A/z/bpwy3du5NXXs74mBgOyw1Voh2QoBeiGU4P+wMzE3k3OJjRXbqQWlLC+JgYTpaVWbpM0clJ0AvRTB7/52Hus6+8QHv0sYN82b8/oU5O7M7P5/rduymRMfbCgiTohWgBHld4MOCLyqGXLx0le/5Rfhw4EC87O9ZkZ3OvPHtWWJAEvRAtpNs/utF/dX+wgcPPH4bn0/khLAxHg4H309LkcYTCYiTohWhBntd60v/T/mCAQ88cwvO1E3xS7XGEq48ds3SJohOSoBeihXWf0J3QVaGg4ODcgwz9sIiXAs1THt8eG8tf8tAS0cYk6IVoBV6TvAh+JxiAxEcSufF7A9N8fCjWmmt27+ZgYaGFKxSdiQS9EK3E+05v+i033y2bMD2BJzaZuLxrV46XlnKlDLsUbUiCXohW5HufLwFLAgBImBLPitjuDHByYm9BATfu2SNTG4s2IUEvRCvze9SP3vN7QwUcunUf/0nuRXc7O349cYIH9++3dHmiE5CgF6IN+M/3p+cjPdFlmoxbEvg6uzf2SrEiJYXXk5MtXZ6wchL0QrQBpRSBSwLxnupNRVEF5Tcf5MOi3gDMSEjgV5ntUrQiCXoh2ohSiqDlQXS/qTvlueV433KURcXelAM37tlDXH6+pUsUVkqCXog2pGwUIatC8LjSg7KsMi68I5PJRe6cLC/nypgYskpLLV2isEIS9EK0MYOdgf7/6Y/bRW6UpJYw9f5CLixyIrGoiBv37KFURuKIFtbgoFdKTVNKHVRKFSmlIpVSF9SzvVEp9a/KfYqVUoeVUg82v2QhOj4bRxvCvgvDNMREcWIRi2ZDQKEta7KzeVhG4ogW1qCgV0pNBJYBzwERwF/Az0opv7Ps9hlwOTAVCAYmANHNqlYIK2Lrasug/w7CMdiRkpgC3ltgj2sxvJ6SwhsyEke0oIae0c8EVmqt39Zax2qtHwBSgftq21gpNQ4YC1yhtf5Na52ktd6qtV7fIlULYSWMnkYG/zYY+1726K35fPaiMzZl8MD+/WzIzrZ0ecJK1Bv0SikjMBT49bRVvwLn1rHbNcB2YKZS6qhSKkEp9W+llKk5xQphjRx6OTDo10HYetji+Hs+K193orzc/NzZJJkTR7QAVd/DEJRSPkAycJHW+o9qy+cBk7TWwbXs819gNLAG+BfQBXgViNZa31DL9lMxd/Hg5eU19LPPPmticyAvLw+TqfP8Puls7QUrbnMs5r+di+CPCTB/GgQArwHl1trms7Da7/NZNKfNY8aMidRaD6ttnW2zqqqbAdDALVrrkwBKqfuBX5RSXlrr9Ooba63fAt4CGDZsmB49enSTP3j9+vU0Z/+OprO1F6y4zaMhq28WMeNjuPBzzfQetrx+XRnveXoyDayzzWdhtd/ns2itNjekjz4DKAe8TlvuBaTVsU8qkHwq5CvFVv73bBdwhejUul7alZBVIQDc8GoZV/2u+OL4cT62cF2iY6s36LXWJUAkcOlpqy7FPPqmNpsAn9P65IMq/3uosUUK0Zl43eRF4MvmB5U89IJmxFZ4D/ghI8OyhYkOq6GjbpYCk5VSU5RSoUqpZYAP8AaAUmqVUmpVte0/ATKB95VSA5RS52EenvmF1lqepSZEPXo91Itej/VClcGzCxRBcTApNlamSRBN0qCg11qvBh4CngR2AedjHjp56uzcj2pdMlrrPOASwA3z6Jv/ABuAO1uobiGsXsCiALxu88KmUPPS42A6XM41u3fLA0tEozX4YqzWejmwvI51o2tZtg8Y1+TKhOjklFIEvxNMSXoJ/HKCV+Yopv67kNucYvl64EAMSlm6RNFByFw3QrRjBjsDAz4fAP3A86hmyRPw69FMFh6SS12i4STohWjnbF1s4Xlw6ONA3ziY9wwsSEzie7k4KxpIgl6IjqArDPp5ELZdbRm1GWYsg1v37iW+oMDSlYkOQIJeiA7CKdiJsO/CUPaKq76HKz+s4Nrdu8mTi7OiHhL0QnQgbue50f/j/qDg7nfA+7sC7ty3j/qmMhGdmwS9EB2M5/WeBC4131D12GLY9/txXjxyxMJVifZMgl6IDqjnjJ74PuCLXSk88xS8vv4Aa0+csHRZop2SoBeiA1JK0fflvnhc5YFrLjw3B6Zu2sORoiJLlybaIQl6ITooZaPo/0l/TMNc8EmFhx8r4+bI3RTLM2fFaSTohejAbJxtCPt+IHZ+RgbshbFP5vFwfIKlyxLtjAS9EB2cfQ97Bv84CFwNjFkPpc+k8mFaXTOIi85Igl4IK2AaaGLQ5wPRNjDpE/ji5X3E5OVZuizRTkjQC2Eluo7rStDr/QB44CXN46uiZaZLAUjQC2FVfO/xxWuGL7blMHVOCQ/9skduphIS9EJYm5CX+uJwRRdcc+Hie0/wSkySpUsSFiZBL4SVUTaKYasHUj7AgV5HQU8+xMbjcjNVZyZBL4QVsjXZct5P4RR3MzDkb/j17hjSi4stXZawEAl6IayUg58Dw78dRKkRxn5bwStP7aJc+us7JQl6IaxY13O70Ott80icS14q5NVVsRauSFiCBL0QVi7kNl/0o16UVBSy9e4XcHP3wGAw4Onpyfz588mT8fZWT4JeiE5g2FM9ucd0P1+VfkZOdhZaazIyMli8eDEjR46UsLdyEvRCdAIvvvQix8tSKKGkxvKioiISExNZsmSJhSoTbUGCXohOYPny5RTVMYVxUVERK1asaOOKRFuSoBeiE8jMzGzWetGxSdAL0Ql4eHg0a73o2CTohegEpk2bhoODQ63rjBi565a727gi0ZYk6IXoBGbNmkVgYOAZYW+njPjgw6DfRlNeWG6h6kRrk6AXohMwmUxs2bKF2bNn4+npWTWO/tb7H2Re9+V477Vj7R0xMtOllZKgF6KTMJlMLFiwgGPHjlFeXs6xY8d4799LyPugL4UOYLc6mz0vJlm6TNEKJOiF6OSmX9aXX//lDED644fI+E1G4FgbCXohOjmDUjwxYzDf3WrAphx2TtxD0aHax9yLjkmCXgiBp9HIVa+Ese0cMJ6oYPNVu+TirBWRoBdCADDaw53St3qR7AMquoioKXFycdZKSNALIarMGRTAd6+4UOgAOZ8c5+iyo5YuSbQACXohRBUbpXj5qoGseNwGgIRZiWT/kW3ZokSzSdALIWrwsbfnvun9+WwiGMpg14TdFKfIYwg7Mgl6IcQZ/s/DA5enfdkZARwrI+r63VSUVFi6LNFEEvRCiFo9GxTIdy84c8wTCrbksv+h/ZYuSTSRBL0QolZGg4F3zx/IC88YKLGDlBUppH2UZumyRBNI0Ash6hTo6Mgj1wTz6gPm93FT48mLlscOdjQS9EKIs7rFy4tuU7z4+XKgsIKY63ZTml1q6bJEI0jQCyHq9Wq/fvw4x4H4flCcWETc7XHoCrmZqqOQoBdC1Mtka8uHQwaw8F+Q4wKZ32VyZMkRS5clGkiCXgjRIBEuLsw8L5BFj5vfH3jiACfWn7BsUaJBJOiFEA02o2dP3Md35aNJQAXsvWkvxalyM1V7J0EvhGgwpRQrQ0L4caotOyOgNL2UvRP3UlEmN1O1ZxL0QohG6W408v6AUJ55CjK6wcmNJzk496ClyxJnIUEvhGi0yz08uCOsJwvmQbkNHFl8hIzvMixdlqhDg4NeKTVNKXVQKVWklIpUSl3QwP3OV0qVKaV2N71MIUR7syggANtRJt6eYn4fd3schQcLLVuUqFWDgl4pNRFYBjwHRAB/AT8rpfzq2c8dWAWsaWadQoh2xt5g4NPQUL67WbHpXCjLLmPPhD1UFEt/fXvT0DP6mcBKrfXbWutYrfUDQCpwXz37vQt8AGxuRo1CiHYqxNmZl/v14/k5kN4D8iLzSJyVaOmyxGnqDXqllBEYCvx62qpfgXPPst80wAtY2JwChRDt21Rvby7p043586HMFpJfTebYF8csXZaoRtX3TEillA+QDFyktf6j2vJ5wCStdXAt+4QBvwMjtdYHlVJPAzdorQfW8RlTgakAXl5eQz/77LMmNgfy8vIwmUxN3r+j6WztBWlze3QSmAJc+CU88BrgDLwJ+Db9mO29za2hOW0eM2ZMpNZ6WG3rbJtVVS2UUvbAauBRrXWDxlxprd8C3gIYNmyYHj16dJM/f/369TRn/46ms7UXpM3tlXNWFpfqaMKj4IKNYHrJxJDNQzDYN21wX0doc0trrTY35DuQAZRj7oapzguobXJqbyAUeL9ytE0ZMA8YUPl+XHMKFkK0T5d07cojvXrywmw45qvI+zuPxEelv749qDfotdYlQCRw6WmrLsU8+uZ0yUAYEF7t9Qawv/L/a9tHCGEFng0IoK+XM089pSm3g+TXkjn+5XFLl9XpNfRvqqXAZKXUFKVUqFJqGeCDOcBRSq1SSq0C0FqXaq13V38Bx4Diyvfy1AIhrJS9wcAn/ftzONTAiqnmZXF3yfh6S2tQ0GutVwMPAU8Cu4DzgSu01ocqN/GrfAkhOrn+zs68GBjIl9fD1vMV5SfLzfPhyMPFLabBV0m01su11v5aa3ut9dDqI3C01qO11qPPsu/TdY24EUJYn2k+Plzh0ZVnZ2lO9lDkbs/lwBMHLF1WpyVz3QghWpxSivdCQnDwsGPuXI22gaMvHSXzp0xLl9YpSdALIVqFl9HIe8HB7BkI791pXhZ3exzFyTJ/fVuToBdCtJoru3XjXh8fPr4J4kYYKM0oJfbWWHS5PG+2LUnQCyFa1YuBgfR1duSJ2RUUdTOQvT6bQ88dqn9H0WIk6IUQrcrZxoaPQ0PJ6QpPPmYeeZP0dBInN520cGWdhwS9EKLVDXd1Zb6/P5HD4PtJBvPzZm/ZS+mJUkuX1ilI0Ash2sTjfn6McnVl2eQK0sNsKT5czL4p+6hvYkXRfBL0Qog2YWsw8GFoKA5GAw/PKaPCZCDjqwxS30q1dGlWT4JeCNFmAh0deaVvX1J94OWZ5jP5/Q/vJz8238KVWTcJeiFEm7rL25urPDz4YYwm6ko7Kgor2HvzXnkEYSuSoBdCtCmlFG8HB9Pdzo7H7yuluLct+VH5HHhcpkhoLRL0Qog2191o5J3gYAqdYNacMrBVHH35KJn/lSkSWoMEvRDCIv7RrRtTvL2JCYEf7rEDIG5yHCXHSixcmfWRoBdCWMzSwED6ODjw8rUlZI2wpzS9lH13yZDLliZBL4SwGBdbWz4MDQUbmDazGLrYkPlDJinLUyxdmlWRoBdCWNR5bm7M9vMjvTssn2WOpMRHE+GghQuzIhL0QgiLW+Dvz2BnZz4/t5TE6x2pKKqAhciQyxYiQS+EsDhj5V2zRqW4/65CKvoY4QAcfFJO61uCBL0Qol0IM5lY2KcPRY4wf04F2gBHXjrCibUnLF1ahydBL4RoN2b26sUFbm78GVTG+tsADbG3xVKaJbNcNocEvRCi3bBRig9CQjDZ2LDwVigZ6khJcgnx98XLkMtmkKAXQrQrfRwdWRoYSIUNzJhVgnI2cPw/xzn2yTFLl9ZhSdALIdqdKd7ejADivMr58REHAOKnx1N0uMiyhXVQEvRCiHZHKcUsoKutLUtGF5A7zpnyk+XE3R6HrpAunMaSoBdCtEsewIqgIFBw77QCDN1tyV6fzZGlRyxdWocjQS+EaLdu7N6dm7p3J8VN88Hj5onPDs49SN7uPAtX1rFI0Ash2rXX+/XD22jkvfBCjt3sgi7RxN4aS0WJ3DXbUBL0Qoh2raudHe8EBwNw9z9zMfSxJz8qn6SnkyxbWAciQS+EaPeu8PDgbm9vchxh2VwDGODwC4c5+ddJS5fWIUjQCyE6hJcCA/F3cOCrwEIO3+MGFea7ZsvyyixdWrsnQS+E6BBcbG15v7IL595rT2IY6EhRYhEHZsuzZusjQS+E6DBGu7szw9eXQjtYOEej7BQpK1LI+jXL0qW1axL0QogO5bmAAIIcHfnNt4i4Ga4AxN0ZR+kJmfisLhL0QogOxcnGhg9CQjAA919+EnWOMyXJJSQ8kGDp0totCXohRIczsvLxg2U2MPfRMgyOBo59fIzjXx+3dGntkgS9EKJDetrfn4HOzmzyLObvR1wAiL8nnpLjJRaurP2RoBdCdEj2BgMfhIRgqxQPjzkJF5goPV4qc9fXQoJeCNFhDXFxYa6fH9oAjz5cgsFkQ8aXGRz7TOaur06CXgjRoc3t3ZsIk4lI9xK2PmYCIGF6AsWpxRaurP2QoBdCdGh2lV04dkox+7yTVIx1oexEGfH3ShfOKRL0QogOL8xk4ml/f1Dw8APF2LjZkPldJukfp1u6tHZBgl4IYRVm9+rFcBcXot1K+OMxZwD2P7Cf4hTpwpGgF0JYBdvKLhx7pXhyZA7ll7pQll3Gvqn7On0XjgS9EMJqhDo780yfPqDgwelF2HSxIevHLNI/7NxdOBL0QgirMrNXL0a6urLXrZT1p7pwZnTuLhwJeiGEVbFRipUhITgYDMwbkUNZZRdO/D2ddxSOBL0QwuoEOznxXLUuHIObDZk/dN5ROBL0Qgir9GDPnpzn6kqsWynrZ1d24Ty4v1PeSCVBL4SwSjZK8X5ICI4GA/NH5VB6sYmyE2UkTEvodF04DQ56pdQ0pdRBpVSRUipSKXXBWba9Tin1q1LquFIqVym1VSl1VcuULIQQDdPPyYlFAQGg4KHpRRhcbMj4JoPj/+lc0xk3KOiVUhOBZcBzQATwF/CzUsqvjl0uAtYC4yu3/wn4+my/HIQQojU84OvLBW5u7O1axoaZTgAk3J/QqaYzbugZ/Uxgpdb6ba11rNb6ASAVuK+2jbXWM7TWz2utt2mt92utFwCRwDUtUrUQQjSQQSneCw7G0WBg3kW5lF7gTGlGaad6IlW9Qa+UMgJDgV9PW/UrcG4jPssFONGI7YUQokX0dXLi+counEfuL0Y5GTi++jgZ32ZYurQ2oeq7KKGU8gGSgYu01n9UWz4PmKS1Dq73Q5SaDjwPDNRaH6pl/VRgKoCXl9fQzz77rFGNqC4vLw+TydTk/TuaztZekDZ3Fi3d5grgYSAamPMFXPY64AGsBNrJl7Y5bR4zZkyk1npYrSu11md9AT6ABi48bfk8YF8D9r8eKAD+Ud+2WmuGDh2qm2PdunXN2r+j6Wzt1Vra3Fm0RpsT8vO144YN2vD7Ov3bsC16Het07F2xLf45TdWcNgM7dB252pA++gygHPA6bbkXkHa2HZVSNwAfArdprb9vwGcJIUSrOdWFU2EDsx8qRRkVae+mkfV7lqVLa1X1Br3WugTzhdRLT1t1KebRN7VSSt2IOeQna62/aE6RQgjRUu739eVCNzf+9i1j233mUTjxd8dTlldm4cpaT0NH3SwFJiulpiilQpVSyzB36bwBoJRapZRadWpjpdRNwMfAHOAPpVSPylfXFq5fCCEaxaAU71aOwnn8ynzKBzpQlFRE0lNJli6t1TQo6LXWq4GHgCeBXcD5wBX6fxdW/Spfp9wL2AKvYB6Geer1VQvULIQQzdK38kaqclt4cmYZ2MDRZUfJ2Zpj6dJaRYPvjNVaL9da+2ut7bXWQ3W1ETha69Fa69GnvVe1vEbXdmwhhGhrD/j6cr6bG1v6lBE92Qk0xN0VR0VJhaVLa3Ey140QolM6dSOVg8HArAkFVPQxUrCngMOLDlu6tBYnQS+E6LT6VU5nXGIPz800n8kfevYQ+XvzLVxZy5KgF0J0ag/27MkoV1fWDCwjfoIjulSz7+596ArrmeFSgl4I0amdms7YXilm3laI9rIl568cUt5IsXRpLUaCXgjR6QU7ObGwTx/yTbDsQfOyA3MOUHSkyLKFtRAJeiGEAB7u1YsRLi58e24Zhy+xpzy3nITp1vGQEgl6IYTgtC6ce4rRLgYyv8/k+Bcd/yElEvRCCFEp1NmZp/39yewGK+8zx+P+B/dTml1q4cqaR4JeCCGqebRXL4a5uPDhZWUcH2KkJK2EA3MOWLqsZpGgF0KIamwNBt4PDsbWRjHrgRK0nSL1zVSy/8y2dGlNJkEvhBCnGWgy8VTv3hzyh+/+aY7J+KnxVBR3zOkRJOiFEKIWc/z8CDeZeH1iObl9bCmILeDwCx1zegQJeiGEqIVdZReOtlc89aB5rvpDzx6iYF+BhStrPAl6IYSoQ7iLC4/7+REVDhuvtEGXaOLvje9wY+sl6IUQ4iye7N2bgc7OvDilnGJ3A9nrs0n74KxPUW13JOiFEOIsjAYD7wUHk+cGL91nvhib+EgiJcdLLFxZw0nQCyFEPYa7ujKrVy9+uwRizzFQllVG4qOJli6rwSTohRCiAZ729yfYyZGFD1ZQboT0VemcWHfC0mU1iAS9EEI0gIONDe+FhJDqC6smmZfF39sxxtZL0AshRAOd6+bGjJ49+fQmSO+tKIwv7BBj6yXohRCiEZ7t0wc/VwcWPWQeYnno2UMUxLfvsfUS9EII0QhONja8ExxMVDj8ehnoEt3u562XoBdCiEYa7e7OfT4+LL8X8t3gxO8nOPbZMUuXVScJeiGEaIIXAgJw87Jn+d3m94kzE9vtvPUS9EII0QQutra8FRzMz/8HewZCSVoJB588aOmyaiVBL4QQTXRZ165M9unBSw9DuQ2kLE8hZ3uOpcs6gwS9EEI0w0uBgRSFGPn8BkBD/H3x6PL2dWFWgl4IIZrB3c6ON4KC+OB2ONYd8iLzSHkjxdJl1SBBL4QQzXRVt25c17s7r95vfn9g7gGK04otW1Q1EvRCCNEClvXty74xtmwZAeUnyzkwq/08UFyCXgghWkA3o5HXgoL494NQbIT0j9rPpGcS9EII0UImeHoyclA3Pq6c9CxhegIVJZaf9EyCXgghWohSitf79eOXSTYc9YWC2AKOvnLU0mVJ0AshREvytrdn8YB+/PtB8/uDC5IoOlJk0Zok6IUQooXd5uWFx+Vd2XAh6IIKEmda9mlUEvRCCNHClFK8GRTEygcMFDrA8S+Ok/VLlsXqkaAXQohW4OfgwGOjAll1m/l93P2WexqVBL0QQrSSqT4+HLvLlaTeULK/iCMvHbFIHRL0QgjRSgxK8ebAEN58SAFwYGESRYfb/sKsBL0QQrSivk5OTLg+gHWjQRVqYh9KaPMaJOiFEKKVzejZk82znCl0gJNfZ7b5hVkJeiGEaGU2SvHKRf35+Hbz+6jpcW16YVaCXggh2kB/Z2f6P+LP4V6gEks4sPRwm322BL0QQrSRWYF+/PiYAwBJCw9RdLRtLsxK0AshRBuxMxiYc9sA/rgQbAs02x6Ka5PPlaAXQog2FOHiglroQ5E9VHyZTfqazFb/TAl6IYRoY7PODeTXyXYAbL8vjoqy1r0wK0EvhBBtzMHGhhsW9CfFB0wJpWxderBVP6/BQa+UmqaUOqiUKlJKRSqlLqhn+4sqtytSSh1QSt3b/HKFEMI6nOvlzuGnulJIIa89+S+6eXTj4osvxtPTk/nz55OXl9din9WgoFdKTQSWAc8BEcBfwM9KKb86tu8D/FS5XQSwCHhVKXV9SxQthBDWYPJNftzpOI0vSj8jMysTrTUZGRksXryYkSNHtljYN/SMfiawUmv9ttY6Vmv9AJAK3FfH9vcCKVrrByq3fxv4AHi0+SULIYR1WPHyy2SVp1JCSY3lRUVFJCYmsmTJkhb5nHqDXillBIYCv5626lfg3Dp2G1XL9r8Aw5RSdo0tUgghrNHy5cspKSmudV1RURErVqxokc+xbcA23QAbIP205enAJXXs0wP4vZbtbSuPl1p9hVJqKjAVwMvLi/Xr1zegrNrl5eU1a/+OprO1F6TNnUVnaHNm5tmHVmZkZLTI16AhQd/qtNZvAW8BDBs2TI8ePbrJx1q/fj3N2b+j6WztBWlzZ9EZ2uzh4UFGRkad67t169YiX4OG9NFnAOWA12nLvYC0OvZJq2P7ssrjCSFEpzdt2jQcHBxqXefg4MB999V1GbRx6g16rXUJEAlcetqqSzGPqqnN5jq236G1Lm1skUIIYY1mzZpFYGDgGWHv4OBAYGAgs2bNapHPaeiom6XAZKXUFKVUqFJqGeADvAGglFqllFpVbfs3AF+l1CuV208BJgMvtkjVQghhBUwmE1u2bGH27Nl4enqilMLT05PZs2ezZcsWTCZTi3xOg/rotdarlVIewJOAN7AbuEJrfahyE7/Ttj+olLoCeBnzEMwU4EGt9ZctUrUQQlgJk8nEggULWLBgQatdl2jwxVit9XJgeR3rRteybAMwpMmVCSGEaBEy140QQlg5CXohhLByEvRCCGHlJOiFEMLKSdALIYSVk6AXQggrJ0EvhBBWToJeCCGsnNJaW7qGGpRSx4FD9W5Yt250ronTOlt7QdrcWUibG6e31tqzthXtLuibSym1Q2s9zNJ1tJXO1l6QNncW0uaWI103Qghh5STohRDCyllj0L9l6QLaWGdrL0ibOwtpcwuxuj56IYQQNVnjGb0QQohqJOiFEMLKdaigV0pNU0odVEoVKaUilVIX1LP9RZXbFSmlDiil7m2rWltKY9qslLpOKfWrUuq4UipXKbVVKXVVW9bbEhr7fa623/lKqTKl1O7WrrGlNeFn26iU+lflPsVKqcNKqQfbqt6W0IQ236KU2qWUKlBKpSmlPlJK9WireptDKXWhUuo7pVSyUkorpSY3YJ8wpdQGpVRh5X7zlFKqSQVorTvEC5gIlAJ3A6HAq0Ae4FfH9n2A/MrtQiv3KwWut3RbWrHNy4A5wDlAX2A+UA5cYOm2tFabq+3nDhwAfgF2W7odrd1m4CtgG3Ap4A+MAEZbui2t1WbgvMqf5Ycr/22PBHYCayzdlga29wrgOeAGoACYXM/2rkAa8B9gYOV+ucAjTfp8S38BGvGF2gq8fdqyBGBRHdu/ACSctuwdYLOl29Jaba7jGNuAlyzdltZuc2XwzQee7oBB39if7XHASaCbpWtvwzY/Chw6bdkdQJ6l29KEtuc1IOjvA3IAx2rLngSSqRxE05hXh+i6UUoZgaHAr6et+hU4t47dRtWy/S/AMKWUXctW2PKa2ObauAAnWqqu1tTUNiulpgFewMLWq651NLHN1wDbgZlKqaNKqQSl1L+VUqbWq7TlNLHNmwBvpdQ/lFk34Cbgp9ar1KJGARu11oXVlv0C+GD+C65ROkTQY57/wQZIP215OlBXH12POra3rTxee9eUNteglJoO9AQ+bNnSWk2j26yUCsN8Jn+r1rq8dctrFU35PgcA5wODgeuB+4HLgZWtU2KLa3SbtdabMQf7x0AJcBxQwO2tV6ZF1ZVfp9Y1SkcJetFISqnrgSXALVrr5kwS124ppeyB1cCjWuuDlq6nDRkAjfl7u1Vr/QvmsL9eKeVl2dJah1KqP+Z+/Gcw/zVwOebAe9OSdXUUtpYuoIEyMF+IOf2H2AvzBYvapNWxfRkdY0a8prQZAKXUDcAq4Dat9fetU16raGybvTFfyHtfKfV+5TIDoJRSZcAVWuvTuwfam6Z8n1OBZK31yWrLYiv/68eZZ4LtTVPa/DiwTWu9pPJ9tFIqH9iolHpCa320dUq1mLry69S6RukQZ/Ra6xIgEvMIg+ouBf6qY7fNdWy/Q2td2rIVtrwmthml1I2Yu2oma62/aL0KW14T2pwMhAHh1V5vAPsr/7/Or1N70cTv8ybA57Q++aDK/7b7v96a2GYnzL8cqjv1vkPkWCNtBi5QSjlUW3YpkAIkNfpolr4C3Ygr1RMx981NwXwWtwzz1eveletXAauqbX9qeOUrldtPqdy/ow2vbEybb8I8ZG0G5j9rT726WrotrdXmWvZ/mo436qax32cTcAT4HBiAeejhbuBzS7elFds8ufJn+z7M1yjOw3xBOtLSbWlge03872SkAJhX+f9+lesXUW2oKOCG+cz9M8zDK6/DPArHuodXVjZ+GubfZsWYzwgurLZuPbD+tO0vwjzWthg4CNxr6Ta0Zpsr3+taXuvbuu62/D6ftm+HC/qmtBkIxjxKpQDzXzavAy6Wbkcrt/kBYE9lm1MxX5jtael2NLCto+v4t7mycv1KIOm0fcKAP4CiyvbOpwlDK7XWMqmZEEJYO2vs2xJCCFGNBL0QQlg5CXohhLByEvRCCGHlJOiFEMLKSdALIYSVk6AXQggrJ0EvhBBWToJeCCGs3P8DgYmHSZtesj4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdwAAAEiCAYAAABTO2OcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA++0lEQVR4nO3dd3wUdf7H8dcnnRB6IEgNSJMmJaBSJNjOrhx62MWGeirq7+70LGc5Pb3z9GxnO3tDPPWs2EsUBem9915CJ5Ce7++PWWCzJiTE7E7K+/l47CPs7HdmPvvZL/vZmfnOjDnnEBERkfCK8jsAERGR2kAFV0REJAJUcEVERCJABVdERCQCVHBFREQiQAVXREQkAlRwaygzSzUzZ2Z/jsC6VprZK+FeT1URyOs9lbzMkYHlplbmcv1Sk/rEr/lswvm5mtmTZja+AvN1N7MCM+te2THJwangViOB/7jleYz0O9byMrMLzOwmv+OIBDNLNLN7zCzd71ikegsU8FHAA4c6r3NuLvAp8NdKDkvKEON3AHJILg55Pgo4Grg8ZPqEyISzX2egqILzXgB0Bx6rtGiqrkTg7sC/M0Jeex0YC+RGMiApl6r42dwIbAA+r+D8zwLjzKyDc25p5YUlB6OCW404594Ifm5mJwD9Q6cHXkuNYFxV6YsI8LYmnXN7/Y6jvJxzhUCh33HIL1W1z8bMYoGLgFdcxS8V+BWwHRgJ3FlJoUkZtEu5FjCzq8xsmZnlmtkUM+tXQptOZvZfM9tqZjlmNsPMzinn8osdrzOz9MCu7fPN7HYzWxtY5jdm1iGoXQZwGtA2eJd40OtmZjeY2ZzA/JvN7EUzSy5h/Z+b2fFmNsnMcoBbQl47zsymBZaz2MwuKeF9JJvZf8xsY6DdXDO7qhzvv7GZ/dPMZpvZbjPLMrMMMxsc1CYVyAw8vTvo/b4SeL3EY31mNtzMpppZduCzecvM2oS0eSUQb0sz+yCw/kwze9jMossRf4aZLTSzI81svJntNbPVZvbHEtomBt7r6kB/WmJmfzazUr9LzKx+YJlPlPBaEzPLM7N/hORhiJn9K/A+9pjZ+2bWtIT5rwl8TjmBz+05M2tcyvvrYWbfB2JZbmYjAq8PMrOfAzleZGa/CZn/F5+NmQ02s7fNbFUgDxvM7PnQdZdX0DqGmtm/zWxLoC/918xSQpoPApLximbwMu4JLCM0/ofNrNDMBu6b5pzLx9vLMqwi8UoFOef0qKYP4BUgp5TXUgEHzACW4BWgP+F96a8BYoPaHoH3a3cBcCtwHfB1YP6LyhHHSrxf2/uepweteypwE96u1D3ApKB2JwbaZOL9Yr8oeH14u73ygReAq4H7gZ3AHCAhZP1LA+/hIbxd7ScHvbY46LUbgcmB+EYELSMhsNx8vN3bNwDfBNr9OeT9OuCeoOdpwArgH4E4bw2sMxfoGWhTF7gmMO//gt7vMYHXRwZeSw1a7kWBafty+Ddgd+DzSw7pB3mB+F8MrOfdwLzXluPzy8DbPbkReAa4Hu/L3AG3BrUz4MvA9BeB3wPvBZ4/W0afGANsAmJC2v0+MH+3kDzMCOT/euBhoAB4O2TeOwNtvwl8Xo8FPr+ZQHwJ7291YFnXA3MDyzwv8L7vCfSNNcAuoEHQ/CV9Nk/g7c69A7gKeBLIBiYCdrB5S/kM9rWbHVjGaOCRwOc6C4gLantHoG3jkGXEBPrKWqBhYNpgvK3zf5SwzjvwDgU18vu7rLY8fA9Aj1/x4ZWv4G4J/g8FnBmYfnrQtC+BeUCdkGV8GfjPa2XEEfrlmh5Yx4KQL4rRgendg6Z9AqwsYZkDAm0vCZk+KDB9VMj6HXBmKbE54PygaXWAhYEv4KiQ2C4NaheN98MjB2gSND204MbvW07QtEZ4BeaFoGnJofMGvVbsixmIxSsE84M/l6DcPhzSDxxwV8gypwNTy9GPMgLz3xby3r/D+5HUIKTv3B0y/8slfK6hfeLkQJtTQ+adAEwvIQ9fU7xw/QuvQO6LpSneD5qvgegS5r++hPd3cdC0zoFpRcDAoOknBaZfWdpnE5iWWEIeLwi0G3SweUv5DPa1m03xHwuXlxDP68COUpbTFa/wvw4kAcvwfojFl9D2/MCyB5TVR/SonId2Kdd87znntgc933caQXvwdocCJwD/BeoGdqsmB3bbfg60BDpVcN2vOefySlt3GX4HZAGfh8S0EK+QDQ1pv9Y591Epy9oMvL3viXMuG2+ruTXQMzD5NLwt7TeC2hXibTXF4+WoRM65XOdcEYCZJZhZE7yCNQXoW473WpI0IAV4JhDvvnVlANMC8YZ6PuT5eMqXa/AKz1NB6ykMPE/kQK5PC7R7PGTeR4JeL81XeFuZ+wf+mVl74Bi84hDqRReoCgHj8XLaNvD8BCAOeDwQ6z6v4/WP0FiygTf3PXHOLQJ2AIudcz8FtZsU+HvQvLnA+ADz1A/0zX2DFSv6mQM854qPiXgtEOfpQdOa4O2xKSmu+Xhbrhfh/RhpjfcjsqRxFvuWkVzCaxIGKrg13+rgJ0HFt1Hgbwe8XYX34BWc4Me+L9JmlbFuDvwHbxTasASd8H6hbyohrpQSYlp+kGUt21cQgywO/E0N/G0LLA358gZvKz243S+YWVTgOOZyvC/2LYE4TwMaHCSug9lXWBaV8NqCEuLJd85tCJm2nfLlGmCTc25XyLSScrTJObcjpN0ivEIcGtN+gby+AZxpZvUCky/C2935VgmzlNV3SsxPYD1LSohlXQl9YCfeLuTg+XeGrKdEZtbazMYGlrET7/NeEXi5op85eLEHx1MQWG5qaAgHWcZjeD/2jgIedM5NL6XdvmW4Ul6XSqZRyjVfaaMr9/1n2/ej61G8c/NKMjdM6z6YKGAr3jG2koT+ws8usVVk3IZ3fPlVvOOKW/He+23A4RGKoaKnZUXSa3jjCH6Ll6sLga+ccxtLaPtr+k5JSlveIa8nMBDtS7zd2g/i/QDag9dnPyf8GzJbOPgPgrZ44zIAehyk3b5lbKmMoKRsKriyb8uwwDn3tQ/rL+3X9TK8QVU/O+eyfuU6DjezqJAtnH27yVcG/q4CeptZdMhWbpeQdiU5F8hwzo0Mnmhm94a0O5QtiVWBv53xvtyDdSkjnopIMbP6IVu5JeXoRDNrELQluK9dVFkxOefmmtl04GIzWxCYLzRH5RWcn31b4gRGS3fEG3QVLj3wPoORzrlXg9bdsRKW3ZGgz9vMYoB2wPdBbRbg5bCJc25r8MxmZnjH1HPw9lDdbWYXuRJOHQws1+EdppEI0C7lWs45txlvcMxVZtYy9PWSTsWoZHuAhoEvimBv4/XPu0qIKdrMyrurFLzdzyOC5q8DXIk3IGx2YPIneFssFwS1i8IbubpvcE5pCgnZIjKzAXjHJ4PtOy+4PLFPxdudfrWZJQQtdzDe8d1PyrGMQxGFNzp933r2Pc/G6x8E1hmFN8As2P8F/o4rx3pexTsmfAveiOv3KxjvV3gjeEeHnJJ0Id4hh8rOT7B9P8hC++wvTqOqgKvNLD7o+SVAQ4rndt8x55KOFd8EDMEb/X0vXr99wsxalNC2L7AwZIyHhJG2cAXgWrz/xLPN7Hm8rctmeMeAuuId5w2XqXjF8DEzmwQUOefGOud+MLOngD+ZWU/gC7zC1wE4B68Qv1LOdSwBnjaz3nhF9iK8LaMLg7Z6n8c7nejFQLvlwNnA8Xijd7f+YqkHfATcY2av4Q3u6RhY1ny849CAN1jLzOYB55nZYrxdzyucc5NCF+icyzezP+Hthh1vZm/g/SAYDazDOwWpMm0EbjTvHN+5eO89Hbg9aGv2E7xCd6+ZtcUbBX0cMBxvsE95Dj28hXdqznC8UcwVOhTgnNtiZvcB9wFfmtkHeAOdrsc7jeaFiiy3nBbi9alHzKwVsA04BWhVScv/zszewjtuewPe5/Fq0OsT8Y4Zn0jxreHOeKeOjXXOvROYdjneKOUXgFOD2sbiFebnKilmKQ+/h0nrUfEH5Tst6M8lvPaLU1MC7V8G1uNtOazDO6Z7fjniWEnJpwWdV0pMI4OmJeJ9mWzFOw7pQua5HO+82b1450fOxfvCbhOy/s8PEtvneIVhGt6utiXBMQS1TcYrvJvwivs84Kqy8oc3WvYfeMU8G2/AysmBz2dlyLxH4Y2EzQks55XA9JGUcPoIXmGaGmi/De8Sg23L0w/wBsK5kvIS0i4Dr4gcifeDIRtvMNEtJbStG8j/2kA/WYp3rDo6pF2xPhHy2oeB93pcCa/ty8PRIdP39an0kOnXBD6n3MDn9h+CTuEKfn+l9Y1SPt9nS4gpNWha50C/2hn4XMbgbVmH9o0SP9eDvO+heKPDt+CN0n8XOKyE9o8E9y28EdyT8EaCh56fe2lg2VcFTTslMK1TWf1Dj8p7WCD5IjWSma3E+7I92e9YqirzrvjV3DnXpay2lbS+d/CuAd7W/XLkcK1k3g1HXsa7EMrP5WjfFu/Y9dnOuc8qsL6P8PYmnX2o80rF6RiuiESMmTXDu4DG6yq2FeecW4W3NX/7oc5r3m35TkXXUI44HcMVkbAzs3bAQLxDBEXA0/5GVP05526o4Hxz0Xe/L5R0EYmEIXi7TNfgHT9f63M8IhGnY7giIiIRENYt3OTkZJeamlppy9uzZw9169attOVVd8pHccrHAcpFccpHccrHAeHIxbRp07Y4535xDYOwFtzU1FSmTp1aacvLyMggPT290pZX3SkfxSkfBygXxSkfxSkfB4QjF2a2qqTpGqUsIiISASq4IiIiEaCCKyIiEgEquCIiIhGggisiIhIBKrgiIiIRoIIrIiISAbq0o4gcVGGRI3N3Lluyctm+N49te/LYmpXH7pwC9uQVeH9zC8gtKCS/0JFfWEReQRHOAebdpT3KjJhoIyE2mviYKBJio6kbF039OrHUT4ilfp0YGtSJo2m9OJomJZBcL47EOH09Sc2iHi1Syznn2Lgrh5Vb9rJq6x5WbdvL6m17Wb8jm407c9i8O5fCoshfArZefAwtGtahRcMEWjaqQ6tGibRLrkv75Lq0aZJIfEx0xGMS+TXKVXDNrB5wHzAMaAbMAG50zk0JY2wiUsl2Zuczb/1O5q/fxeJNu1m8KYulm7PIyi046HzJSXEkJ8XTuG7c/kf9hFiSEmKoGx9DUnw0CTHRxMVEERsdRUy0EWWGc+BwOAf5hUXk5BeRW1BIbn4RWbne1vGunHx2ZuezY28emVl5bNmdS2ZWLrtzC1i0aTeLNu3+RTxRBq0aJdK5eT26NK8X+Fuf9sl1iYqycKVP5Fcp7xbuC0BP4FJgLXAR8LWZdXXOrQtXcCJScTn5hcxdt5Ppq7czc80O5q7bxepte0ts27huHO2S69K2cSJtm9SlbZNEWjSsw2ENEmhWPz7iW5POOXZm57N2ezbrd2Szbkc2q7buZeXWPazYsoc1ga3w1dv28tX8TfvnS4qPoXvL+vRs1ZAjWzWkb9tGNG+QENHYRUpTZsE1szrAcGC4cy4jMPkeMzsDuBbdxFikSsjKLWDKym38vGwrX83KZs1XX5BfWHxXcHxMFEccVp9uLerTuXk9OjarR6eUJJokxfsUdcnMjIaJcTRMjKN7ywa/eD2voIiVW/ewcONuFm3cxaKNu5m3fhcbdubw8/Jt/Lx82/62rRrVoXWdPNbVWcWAw5NJbZKImbaCJfLKvD1fYHfyLuBk59wXQdN/BAqcc+kh7UcBowBSUlL6jh07ttKCzcrKIikpqdKWV90pH8XVtnwUFjmW7yxizpZC5m0pZMWuIoIPtRrQMsk4vGE0HRpGkdogmhZ1jegavMt1R24RK3cWsWJnEct2FrFsRyHZIXvLGycYXZtE07VJNN2To6kfV3PzUZra9n/lYMKRi6FDh05zzqWFTi/X/XDNbAJQCJwHbATOB14FljrnOpc2X1pamtPdgsJH+SiuNuRj+548vl24mW8Xbmb8kkx25RyoJtFRRveWDTimfRMSs9Yy8owh1E+I9TFa/xUWORZt3M2YryaxPaYJE5ZtYfve/P2vm8GRrRpyXJdmDO3cjO4t69eKrd/a8H+lvMJ0t6ASC255j+FeDLyEd/y2EJgOvAX0rbQIRaRE63Zk8/ncjXw1fyNTVm4vNmI4tUkiQzo15dhOTenfrjH1AgU2I2NjrS+24P0I6dqiPie0jSU9vQ9FRY4FG3cxYelWfliSyaTl25i5Zgcz1+zgX18tpkWDBE7q1pyTuzenX2rjGr03QCKvXAXXObcMGGJmdYH6zrkNZvY2sDys0YnUUht2ZvPpnI18Mns9M1bv2D89JsoY3DGZE45IIb1zU9o20U3ED0VUlNGtRQO6tWjAVce2Z29eAT8t3cp3izbzzYJNrN+ZwysTVvLKhJU0rhvHKd2bc+aRLeiX2lijn+VXO6TzcJ1ze4A9ZtYI+A1wS1iiEqmFdufk89mcjbw3fS2TVhwY9FMnNprjujTjpG4ppHduRoM62nKtLIlxMZzYNYUTu6ZQdFZ3Zq3dwefzNvLF3I2s3LqXNyet5s1JqzmsQQKn9zyMYb1b0bVFfb/DlmqqvOfh/gbvMpALgQ7APwP/fjl8oYnUfEVFjh+XbuHdaWv5Yt5GcguKAG808dDOzTj9yMM4rkszXXUpAqKijN5tGtG7TSP+fHIXFmzYzUez1vPxrPWs25HN8+NX8Pz4FXRrUZ9z+7bizF4taVw3zu+wpRop7//iBsCDQCtgG/AecIdzLv+gc4lIiTbtyuGdqWsYO2UNa7dn759+VLvGDO/TilN6NN9/PFYiz8w79tu1RX1u+U1nZqzZzvsz1vHRzPXMW7+Leevn87dPF3BSt+ZceFQbjmnfpFYMtpJfp7zHcP8L/DfMsYjUaM45JizbyqsTVvLNws37Bz+1bFiH36W15rd9WtK6caLPUUqoqCijb9vG9G3bmDtP68rXCzbx7rS1/LA4k3GzNzBu9gbaJ9flgqPacE7fVjRM1FavlEz7qUTCbE9uAf+bsY7XJqxkyeYswBv8dEr35pzXvw2DOyRrQE41kRAbzek9W3B6zxZs2JnN2MlrGDtlNcu37OH+cQt4+MtF/LZPKy4bkErHlHp+hytVjAquSJhs2pXDyz+t5M1Jq9gdOF+2Wb14LjyqLecf1Zpm9XTJwerssAZ1uPnETtxwXAe+XbiZ139exfglWxgzaTVjJq1mcMdkLh/UjvROTbW7WQAVXJFKt2jjbp4fv5wPZ67bf2nFvm0bcemAVE7u1py4GN2GuiaJiY7ipG7NOalbc5Zu3s3LP63kf9PXMX7JFsYv2UKX5vW4ZsjhnNbzMGKj9dnXZiq4IpVk1pod/Pu7pfsvph9lcGqP5lw5uD192jTyOTqJhA7N6vG3YT34028689bkNbz80woWbtzNTW/P5J9fLGLUse0Z0a81CbG6tWBtpIIr8itNWbmNJ75ZwvglWwDvlJ4R/VpzxaB2ujBFLdUwMY5r0w/n8kGpfDBjHc/9sJzlmXu4+6N5PPXdUkYd254Lj2pLnTgV3tpEBVekgqat2s6jXy3mx6Veoa0bF81Fx7TlykHtaVqvat19R/wRHxPNiH5tOLdva76cv4knv13CvPW7uH/cAp79fhmjjm3PxUenqvDWEiq4Iodo9lrvursZizIBqBcfw2WD2nH5wFSdEiIliooyTu7enN90S+HbhZt54pslzFq7kwc+XcgL41dww3EdGNGvjY7v13AquCLltGLLHh7+YhHj5mwAvC3aywa246rB7WmQqItUSNnMjOOPSOG4Ls3IWJzJI18uYu66Xfzlw3k898NybjqhE8N6t9RNE2ooFVyRMmzencMT3yxh7OQ1FBQ54mOiGDkglauHHK5L+0mFmBlDOzcjvVNTPp+7kUe+WszSzVn88Z1ZvPjjCm47pQvHdmrqd5hSyVRwRUqRk1/I8z8s55nvl7E3r5AogxFprbnpxI4c1qCO3+FJDWBmnNLjME7q1pwPZqzjkS8XsWDDLi55aTKDOyZz2ylH6GYJNYgKrkgI5xwfzVrPPz5byPqdOQCccEQKt57cWVcPkrCIjjKG923FaT0P45UJK3nqu6WMX7KFH5eO57x+bfjDSZ1ITtJAvOpOBVckyOy1O7j7o3n770Hb9bD63Hn6EQw4PNnfwKRWSIiN5pohhzMirTVPfLuE1yeu4q3Jq/lk1npuPKEjlxyTqoFV1ZgKrgiwbU8e//xiIWOnrME5SE6K50+/6cQ5fVtrAItEXKO6cdx9RjcuPKot94+bT8aiTO4ft4Axk1Zz71ndGNxRx3erIxVcqdUKixxvTV7Nw18uYsfefGKijCsGt+OG4zuSFK//HuKvDs2SeOWy/ny3cDP3fTKf5Vv2cPGLkzmt52H85bSuNG+g63FXJ/pGkVpr/vpd3P7+HGau2QHAoA7J3HNmNzo0S/I3MJEQQ7s0Y2CHZF74cTlPfrOUcbM3kLFwMzef2ImRA1KJ0TWaqwUVXKl19uYV8NjXS3jxxxUUFjlS6sdz1+ndOLVHc93VRaqsuJgofp/egTOPbMFfP57Pl/M3cf+4Bbw/Yx3/GN6T7i0b+B2ilEEFV2qV7xdncvv/5rBuRzZmMHJAKn84qRP1EnThCqkeWjVK5D+XpPHtwk385YN5zFu/izP//SNXDGrHzSd2IjFOX+tVlT4ZqRV27s3nvnHzeXfaWgC6tajPA8N6cGTrhv4GJlJBx3VJ4aibm/Cvrxbz8k8reH78Cj6bu5F/DO/JwA4aVV8VqeBKjff53I385cO5ZO7OJS4miv87sRNXDmqn415S7dWNj+Evp3flrF4t+PN7c5i/YRcXvjCJC45qw22ndNGemypG3zhSY+3Ym8eNY2dwzRvTyNydS1rbRnx242CuGXK4iq3UKD1bNeTD6wfyx5M6ERttjJm0mpMfG8/4JZl+hyZB9K0jNdJ3izZz0qM/8OHM9dSJjeaeM7ry36uP4fCmGoEsNVNsdBTXH9eRT24YTI+WDVi3I5uLX5zMnR/MYW9egd/hCeXYpWxm0cA9wEXAYcAG4E3gHuecPkXxV24WTHgCprzAkL3bcJMa8029Mxi9+lj2kkBa20Y8fO6RpCbrRvBSO3RuXo/3fz+A535YzmNfL+aNn1fz09KtPPK7I+nTppHf4dVq5TmGeytwHXApMAfoCbwK5AL3hS80kTLkZsELJ8D2FVCQgwFkb2XQ3jf5IO5rxg8Zy8ih3XWlKKl1YqKjuG5oB9I7N+X/3p7Fok27OeeZCVw3tAOjj+/od3i1Vnl2KQ8APnbOfeycW+mc+wj4CDgqvKGJlGHCE/uLbbAEy6dDbCZXRH2sYiu1WrcWDfjw+oGMOrY9Dnjy26Wc8+xENu8t8ju0Wqk8BfdHYKiZdQEws67AccCn4QxMpExTXvhFsd0nqjAXprwY4YBEqp6E2GhuP/UI3rrqaFo0SGDWmh3c9VM2789Y63dotY455w7ewLv0zv3AbUAh3m7ovznn7iyl/ShgFEBKSkrfsWPHVlqwWVlZJCVp0Ms+tT0fQzLOxii9/zqM79M/iFxAVUht7xuhlA/PnnzHy3NzmbqpEIBjWkRzSdd46sTU3j1B4egbQ4cOneacSwudXp6Cex7wT+BPwDygF/A48Cfn3EE3IdLS0tzUqVMrGvMvZGRkkJ6eXmnLq+5qaz5y8gt58NMFjJ5+Mk1sd+kNE5PhlmWRC6wKqa19ozTKxwHOOf765teMXVRIdn4h7ZLr8tQFfWrtje7D0TfMrMSCW55dyv8EHnbOjXXOzXHOvQ78C2+LVySiVmzZw2+fnsCrE1cxpuhECqJKuSl3TAL0uyKywYlUA2bGkFaxfHzDILo0r8eKLXs4++mfeHPSKsraAJNfpzwFNxFvV3KwwnLOK1JpPpuzgTOe/JH5G3bRtkkiQy+/n5gm7b3iGiwmARq1gwGj/QlUpBro0CyJD64byPn925BXUMQd78/lhrdmkJWrsz3DpTxF82Pgz2Z2mpmlmtkw4P+A98Mbmognv7CI+z6Zz7VvTicrt4BTezTnkxsG0b1dS7jyaxh4IyQm4zBvN/LAG73p8TpmJ3IwCbHRPPjbHjx+Xi/qxkXzyewNnPnvH1my6SCHaqTCynMe7g1459s+DTTDu/DF88BfwxiXCACbduVw3ZvTmbpqOzFRxu2nHsFlA1MP3EYvPgmG3g5Db+d7HacTqZCzerWke8sG/P6N6SzatJuznvqJfwzvyRlHtvA7tBqlzC1c59xu59xNzrm2zrk6zrn2zrnbnXMln48hUkmmrNzGaU/8yNRV22leP4G3rz6aywe10z1rRcLg8KZJvH/dAM7q1YK9eYXc8NYM7v14HvmFOme3sug4rFQ5zjlem7iS8//zM1uychlweBPGjR5E37aN/Q5NpEZLjIvhsRG9uPfMbsRGGy//tJILX5jElqxcv0OrEVRwpUrJyS/klndnc9eH8ygoclw1uB2vXd6fJkmljEYWkUplZlw6IJWxo46hWb14Jq/YxplP/sjcdTv9Dq3aU8GVKmPTrhxG/Odn3pm2loTYKB4/rxd3nNZVt9IT8UHfto345IZB9G7TkPU7cxj+zAQ+mLHO77CqNX2TSZUwa80Ozvz3j8xas4NWjerwv2sHclavln6HJVKrNaufwNhRRzMirTW5BUXc9PZMHvx0AYVFOl+3IlRwxXcfzlzH756byKZdufRv15gPrxtYa696I1LVxMdE8/fhPbjvrG7ERBnP/bCcq1+fpvN1K0AFV3xTVOR45MtF3Dh2JrkFRZzfvzVvXHGUjteKVDFmxsXHpPLa5f1pUCeWrxds4pxnJrB2+16/Q6tWVHDFFzn5hYweO4Mnv11KdJRxzxldeWBYD+Ji1CVFqqoBHZL54LqBtE+uy8KNuzn7qZ+Ytmq732FVG/p2k4jbkpXL+c//zCezN5AUH8NLI/sxcqDOrxWpDtol1+X93w9kUIdktmTlccHzPzNu9ga/w6oWVHAlopZs8n4Vz1i9g5YN6/DutccwpFNTv8MSkUPQIDGWVy7rx4VHtSG3oIjrxkzn2e+X6eYHZVDBlYiZuGwrv31mAmu3Z3Nkqwa8f90AujTX4CiR6igmOor7z+7Obad0AeDvny3kjg/mUqArU5VKBVci4qNZ67n0pcnszingpK4pgZPqE8qeUUSqLDPj6iGH8/SFfYiPiWLMpNVc+dpU9mgEc4lUcCWsnHP854dljH5rBnmFRYwckMozF/WlTly036GJSCU5tcdhjLnqaBrXjSNjUSbnP/+zLgdZAhVcCZuiIse9H8/ngU8XAnDHqUdw9xldiY7S4CiRmqZv20a8e80xtG5ch9lrdzL8mQms3LLH77CqFBVcCYu8giJufHsmr0xYSVx0FE+c35urjm2vkcgiNVj7pkn879qBdG9Zn1Vb9zL8mQnMWrPD77CqDBVcqXRZuQVc8eoUPp61nqT4GF65rB9n6r6aIrVC03rxjB11DIM7JrN1Tx7nP/8zPy7Z4ndYVYIKrlSqrVm5XPD8z4xfsoXkpDjGjjqaAR2S/Q5LRCJo3/n1ZwfurXv5K1P4dI7O1VXBlUqzfkc25z43kdlrd9KmcSLvXjOA7i0b+B2WiPggNjqKf/2uFyMHpJJX6J2rO2bSar/D8pUKrlSK5ZlZnPvsRJZn7qFL83q8e+0xpCbX9TssEfFRVJRx9xld+b8TO+Ec3P7+HJ7OWOp3WL5RwZVfbcGGXfzuuYms25FNnzYNeVvn2IpIgJkx+viO3Hd2d8zgoc8X8dDnC2vlValUcOVXmb56OyOem8iWrDwGd0zmjSuPokFirN9hiUgVc/HRbXlsRC+io4ynM5Zx78fzKapl99VVwZUKm7hsKxe9MIldOQX8plsKL1yaRmJcjN9hiUgVdVavljxzYR/ioqN4ZcJKbnlvdq26mX2ZBdfMVpqZK+ExLhIBStX0w+JMRr48mb15hQzr3ZKnLuhDfIyuHiUiB3dSt+a8ODKNOrHRvDttLaPHziC/llx/uTxbuP2Aw4IefQAH/DeMcUkV9s2CTVz56lRyC4o4r19rHjn3SGKitbNERMpncMemvHZFf+rFxzBu9gauHzOdvIKaX3TL/JZ0zmU65zbuewCnArtQwa2VPp+7gWvemEZeYRGXHtOWB4b1IEqXahSRQ9QvtTGvX3kU9RNi+GLeJn7/5jRyCwr9DiusDmmzxLzr8l0BvOGcyw5PSFJVjZu9gevGzCC/0HHV4Hbcc2Y3FVsRqbBerRsy5qqjaZgYy9cLNjPqtWnk5NfcomuHMjTbzE4CvgB6OedmldJmFDAKICUlpe/YsWMrI04AsrKySEpKqrTlVXeRzMfkDQU8OzuXIgent49leMfYKnddZPWPA5SL4pSP4qpaPtbsLuKhydnszoduTaK4sU8CcdGR+X4JRy6GDh06zTmXFjr9UAvuO0Bb51z/8rRPS0tzU6dOLX+UZcjIyCA9Pb3SllfdRSofH89az01vz6SwyHHDcR34vxM7VbliC+ofwZSL4pSP4qpiPhZv2s0Fz09iS1Yugzsm8/wlaSTEhn8gZjhyYWYlFtxy71I2s2bAWcDzlRmYVG3BxXZ0FS62IlK9dUqpx1tXHUVyUhzjl2xh1Os1b/fyoRzDHQnkAm+FJxSpaj6ds+FAsT2+Izer2IpIGHVMqceYq46mSd04flicydU1rOiWq+AGBktdCYx1zmWFNySpCr6ct5HRb83Yv2V78wkdVWxFJOw6BYpu47pxfL84k2veqDmjl8u7hZsOdES7k2uF7xZu5rox0ykoclwz5HBt2YpIRHVuXo8xVx1F47pxZCzK5IYxNePiGOUquM6575xz5pybHO6AxF/jl2Ry9RvTyC90XD6wHbee3FnFVkQirkvz+rx+RX/qJ8Tw5fxN/N9/Z1X7y0Dq8kCy36TlW7nqtankFRRx8dFt+cvpR6jYiohvurVowGtXHEVSfAwfz1rPLe/OrtY3PFDBFQBmrdnBFa9OJSffu1zjvWd2U7EVEd/1at2Qly/rR53YaN6bvpY7P5xbbW/tp4IrLNy4i0temkxWbgFnHtmCv+lyjSJShfRLbcyLl6YRHxPFmEmr+ftn1fN+uiq4tdyKLXu46IXJ7MzO54QjUnjkd0cSrWIrIlXMgA7JPHtRX2KijOd+WM7TGcv8DumQqeDWYut3ZHPh8z+zJSuXgR2a8O8LehOru/6ISBU1tEszHh3RCzP45xeLeHXCSr9DOiT6dq2ltmblctGLk1i/M4e+bRtF7DJqIiK/xhlHtuCBYT0AuPujebw3ba3PEZWfCm4tlJVbwMiXp7A8cw9dmtfjpZH9SIyL8TssEZFyOb9/G+449QgAbnlvNl/P3+RzROWjglvL5OQXctWrU5mzbidtGify2uX9aVAn1u+wREQOyVXHtuf6oR0oLHJcN2Y6k1ds8zukMqng1iIFhUXcOHYGE5dvpWm9eN644iia1U/wOywRkQr5w0mdOL9/G3ILirji1Sks2LDL75AOSgW3lnDO8ZcP5/LFvE3UT4jhtcv706ZJot9hiYhUmJlx/9ndOblbc3bnFHDJS5NZvXWv32GVSgW3lnjs6yW8NXkN8TFRvDSyH0ccVt/vkEREfrXoKOOx83pxTPsmZO7O5eKXvHvqVkUquLXAGz+v4vFvlhBl8OT5vUlLbex3SCIilSYhNpr/XNKXbi3qs2rrXi5/ZQp7cgv8DusXVHBruM/nbuSuD+cC8LdhPTipW3OfIxIRqXz1EmJ5+bJ+tG5ch9lrd3LdmOlV7g5DKrg12NSV2xg9dgZFDm4+wRtcICJSUzWrl8Brlx+4rd/t/5tTpS4BqYJbQy3LzOLKwJ1/LjiqDaOP7+B3SCIiYdcuuS4vXppGQmwU70xby7++Wux3SPup4NZAmbtzGfnyZHbszef4Ls34q+78IyK1SO82jXjqgj5ERxlPfruUtyav9jskQAW3xtmbV8AVr05hzbZserZqwJMX9CZG10cWkVrm+CNSuP/s7gDc+cFcvl+c6XNEKrg1SkFhETeMmcHstTtp3bgOL16qSzaKSO11fv82/D79cAqLHL9/Yxrz1/t7YQwV3BrCOcdfP5nPNws30zAxllcu60/TevF+hyUi4qs/ntSZM49swZ68Qi5/ZQobdmb7FosKbg3x8k8reW3iKuKio/jPxWkc3jTJ75BERHwXFWX889ye9E9tzMZdOVz28hR25+T7E4sva5VK9fX8Tdw3bj4AD53Tk/7tdGELEZF94mO8C2O0b1qXhRt3c8NbMyjw4RzdchVcMzvMzF41s0wzyzGz+WY2JNzBSdnmrtvJ6LEzcA5uOqEjZ/du6XdIIiJVTsPEOF4e2Y9GibFMXrSGH5//AzzUniEZZ8ND7eG7ByA3K6wxlDmixswaAj8BPwKnAZlAe2BzWCOTMm3PKeLWV6ewN6+QYb1bcuPxHf0OSUSkymrbpC4vnH8E9V4/mTYbNoHlYwB7t8JPj8P8j+DKryE+PIfkyjOE9RZgg3PukqBpK8ISjZTb3rwCHpuey6ZdRfRPbczfh/fQubYiImXou/Z1CmIyiSkKOY5bkAPbV8CEJ2Do7WFZd3l2KZ8NTDKzt81ss5nNNLPrTd/uvikqcvzxnVms2lVEm8aJPHtxX+Jjov0OS0Sk6pvyAjFFpdxNqCAHprwYtlVbWdeZNLOcwD8fBf4L9AKeBP7snPt3Ce1HAaMAUlJS+o4dO7bSgs3KyiIpSaNv/7ckj4+W5ZMQ7bjrmERaJGnsG6h/BFMuilM+iqvN+RiScTZG6XXPYXyf/sGvWsfQoUOnOefSQqeXp+DmAVOdcwOCpj0ADHPOHXGwedPS0tzUqVMrGPIvZWRkkJ6eXmnLq44+nLmOG8fOJMrgpj7xjD73BL9DqjLUPw5QLopTPoqr1fl4qL13zLY0iclwy7JftQozK7HglmfTaAMwP2TaAkC3nomwmWt28Kd3ZwPwl9O70rOpriIlInJI+l0JMQklvxaTAP2uCNuqy1NwfwI6h0zrBKyq/HCkNJt25TAqcPef8/u3YeSAVL9DEhGpfgaMhkbtfll0YxK86QNGh23V5Sm4jwJHm9kdZtbBzM4FRgNPhS0qKSYnv5BRr09j8+5c+rdrzF/P0t1/REQqJD7JO/Vn4I2QmIzDvN3IA28M6ylBUI7TgpxzU8zsbOAB4C/A6sDfp8MWleznnOOO9+cya80OWjaswzMX9iFWd/8REam4+CTv1J+ht/N9BI9nl+sgoHNuHDAuzLFICV78cQXvTV9Lndhonr8kjSZJuiGBiEh1pE2lKuyHxZk88OkCAB753ZF0bVHf54hERKSiVHCrqNVb93LDWzMocnDDcR04tcdhfockIiK/ggpuFbQ3r4BRr09lZ3Y+x3dpxs0ndPI7JBER+ZVUcKsY5xy3vjeHhRt30z65Lo+e14uoKI1IFhGp7lRwq5gXxq/g41nrqRsXzXMX96V+QqzfIYmISCVQwa1Cflq6hQc/OzBIqmNKPZ8jEhGRyqKCW0Ws35G9f5DU9UM7cHJ3DZISEalJVHCrgNyCQq59czrb9uRxbKem3HyiBkmJiNQ0KrhVwF8/nr//SlKPj+hFtAZJiYjUOCq4Pntv2lrenLSauJgonrmoD43qxvkdkoiIhIEKro/mr9/F7e/PAeCvZ3ajZ6uG/gYkIiJho4Lrk105+Vz75jRyC4oYkdaa8/rr9sIiIjWZCq4PnHP86Z1ZrNq6l24t6nPvWd38DklERMJMBdcHL/64gi/mbaJeQgxPX9iHhNhov0MSEZEwU8GNsGmrtvH3zxYC8M9zjqRtk7o+RyQiIpGgghtBW7NyuX7MDAqKHFcOasfJ3Zv7HZKIiESICm6EFBU5bv7vLDbszKFPm4bcekoXv0MSEZEIUsGNkGe+X8YPizNplBjLvy/oQ2y0Ui8iUpvoWz8Cpqzcxr++WgzAv0b0okXDOj5HJCIikaaCG2bb9uRxw5gZFBY5rh7SnqGdm/kdkoiI+EAFN4yKihx/fGcWG3fl0LdtI/54Ume/QxIREZ+o4IbRCz8u59uFm2lQJ5Ynzu+t47YiIrVYmRXAzO4xMxfy2BiJ4KqzmWt28NDniwB45NwjaanjtiIitVpMOdstAtKDnhdWfig1x66cfG54azoFRY7LB7bjhK4pfockIiI+K2/BLXDOaau2HJxz3PH+XNZsy6Zbi/rceoqO24qICJhz7uANzO4BbgF2ALnAJOB259zyUtqPAkYBpKSk9B07dmylBZuVlUVSUlKlLS8cflibz0tz84iPhnsH1KF53fAdt60O+Ygk5eMA5aI45aM45eOAcORi6NCh05xzaaHTy1NwTwHqAQuBZsCdQBegm3Nu68HmTUtLc1OnTq1w0KEyMjJIT0+vtOVVtqWbd3PGkz+RnV/II+ceyfC+rcK6vqqej0hTPg5QLopTPopTPg4IRy7MrMSCW+YuZefcZyEL+hlYDlwK/KvSIqzmcgsKueGtmWTnFzKsd8uwF1sREaleDnl/p3MuC5gHdKz8cKqvhz5fxIINu2jbJJH7zu7udzgiIlLFHHLBNbMEvF3KGyo/nOrph8WZvPjjCmKijMfP601SfHnHoomISG1RnvNwHzazIWbWzsyOAt4F6gKvhj26amBrVi5/eGcWADef2IlerRv6G5CIiFRJ5dkUawW8BSQDmcDPwNHOuVXhDKw6cM5xy7uzydydS/92jblmyOF+hyQiIlVUeQZNnReJQKqjN35exTcLN1M/IYZHR/QiOsr8DklERKooXdy3gpZu3s394xYA8MBve+jSjSIiclAquBWQV1DETW/PJLegiOF9WnF6zxZ+hyQiIlWcCm4FPPb1Yuau20WrRnW458yufocjIiLVgAruIZq8YhvPfL+MKINHR/SiXkKs3yGJiEg1oIJ7CHbn5HPz2zNxDq5NP5x+qY39DklERKoJFdxDcM9H81m3I5seLRtw4/Gd/A5HRESqERXccvp87gbem76WhNgoHh3Ri7gYpU5ERMpPVaMcMnfncvv7cwG47ZQj6NBMt7USEZFDo4JbBucct/1vNtv25DGoQzIXH93W75BERKQaUsEtwzvT1vL1gs3US4jhoXN6EqWrSYmISAWo4B7Emm17+evH8wG498xutNDVpEREpIJUcEtRVOT44zuzyMot4ORuzRnWu6XfIYmISDWmgluKVyasZNKKbSQnxfG3Yd0x065kERGpOBXcEizPzOKhLxYC8MCwHjRJivc5IhERqe5UcEMUBnYl5+QXMax3S07q1tzvkEREpAZQwQ3x4o/Lmb56B83qxXPPGd38DkdERGoIFdwgSzfv5uEvFwPwj+E9aZCoGxOIiEjlUMENKCgs4g/vzCavoIjfpbViaJdmfockIiI1iApuwPPjVzBrzQ4Oa5DAnafrHrciIlK5VHDxdiU/+rW3K/nvw3tSX/e4FRGRSlbrC25hkeNP7x7YlTykU1O/QxIRkRrokAuumd1mZs7M/h2OgCLtpR9XMGP1DprXT+CO07QrWUREwuOQCq6ZHQ2MAmaHJ5zIWp6ZxcNfLgLgwd/2oEEd7UoWEZHwKHfBNbMGwJvA5cD2sEUUIUVFjlvfm01uQRHD+2hUsoiIhNehbOH+B3jXOfdduIKJpNcmrmTKyu00qxfPXRqVLCIiYWbOubIbmV0FXAMc7ZzLN7MMYK5z7voS2o7C2+1MSkpK37Fjx1ZasFlZWSQlJf3q5WTuLeLOn7LJLYQbesfTNyWmEqKLvMrKR02hfBygXBSnfBSnfBwQjlwMHTp0mnMuLXR6mZXGzDoDDwCDnHP5ZbV3zv0Hb2uYtLQ0l56efujRliIjI4NfuzznHJe8NJncwmxO63kYfxjRp3KC80Fl5KMmUT4OUC6KUz6KUz4OiGQuyrNpdwyQDMwLukVdNHCsmV0D1HXO5YYpvkr33vR1jF+yhYaJsbpWsoiIREx5Cu4HwNSQaS8DS/C2fPMqOaaw2bw7h/s+mQ/AXad3pWk93XZPREQio8yC65zbAewInmZme4Btzrm54QkrPO75aB47s/MZ0qkpw3q39DscERGpRWrNlaY+n7uRT+dsJDEumr8N607Q7nEREZGwq9DwXOdceiXHEVa7cvK560NvY/yW33SmVaNEnyMSEZHaplZs4f7js4Vs3p1L7zYNufiYVL/DERGRWqjGF9wpK7fx5qTVxEYbf/9tT6KjtCtZREQir0YX3NyCQv78nnfZ52uGHE7n5vV8jkhERGqrGl1wn/5uGcsy99C+aV2uG9rB73BERKQWq7EFd+nm3TydsRSAB4f1ICE22ueIRESkNquRBbeoyHHb/+aQX+g4v39rjmrfxO+QRESklquRBfedaWuYsnI7yUlx/PnkI/wOR0REpOYV3C1ZuTzw6UIA/nJ6Vxok6qbyIiLivxpXcP82bgE7s/MZ3DGZM49s4Xc4IiIiQA0ruOOXZPL+jHXEx0Txt7N76PKNIiJSZdSYgpuTX8idH3iXb7zxhI60aaLLN4qISNVRYwruU98tZdXWvXROqcdVg9v7HY6IiEgxNaLgLt2cxbPfLwPgb8O6ExtdI96WiIjUINW+Mjnn+MsHc8kvdIxIa01aamO/QxIREfmFal9wP5y5nonLt9IoMZY/n9LF73BERERKVK0L7s69+dw/bj4At516BI3qxvkckYiISMmqdcH955cL2ZKVR7/URpzTp5Xf4YiIiJSq2hbcmWt28Oak1cREGfef3YMo3edWRESqsGpZcAuLvIFSzsEVg9rpPrciIlLlVcuCO2byauas28lhDRIYfXxHv8MREREpU7UruFuycvnn597NCe46vSt142N8jkhERKRsZVYrM7sOuBpIDUyaB9zvnBsXxrgOyM2CCU/AlBcYsncbe8bX5/L845l7+KWc3L15REIQERH5tcqzebgWuBVYgrdFfCnwgZn1dc7NDmdw5GbBCyfA9hVQkIMBSYU7uSbmY6Ky52F530J8UlhDEBERqQxl7lJ2zn3onPvMObfUObfYOXcHsBs4JuzRTXhif7ENlmD5xO1a5b0uIiJSDRzSMVwzizaz84AkYEJ4Qgoy5YVfFNv9CnJgyothD0FERKQymHOu7EZmPYCJQAKQBVxY2jFcMxsFjAJISUnpO3bs2AoHNyTjbIzS43MY36d/UOHlV3dZWVkkJWmX+j7KxwHKRXHKR3HKxwHhyMXQoUOnOefSQqeXt+DGAW2ABsA5wFVAunNu7sHmS0tLc1OnTq1YxAAPtYe9W0t/PTEZbllW8eVXcxkZGaSnp/sdRpWhfBygXBSnfBSnfBwQjlyYWYkFt1y7lJ1zeYFjuNOcc7cBM4GbKzXCkvS7EmISSn4tJgH6XRH2EERERCpDRc/DjQLiKzOQEg0YDY3a/bLoxiR40weMDnsIIiIilaHMgmtmfzezwWaWamY9zOxBIB14M+zRxSfBlV/DwBshMRmHebuRB97oTdcpQSIiUk2U5zzc5sAbgb87gdnAKc65L8IZ2H7xSTD0dhh6O9/ruIOIiFRTZRZc59zICMQhIiJSo1W7aymLiIhURyq4IiIiEaCCKyIiEgHluvBFhRdulgmsqsRFJgNbKnF51Z3yUZzycYByUZzyUZzycUA4ctHWOdc0dGJYC25lM7OpJV29o7ZSPopTPg5QLopTPopTPg6IZC60S1lERCQCVHBFREQioLoV3P/4HUAVo3wUp3wcoFwUp3wUp3wcELFcVKtjuCIiItVVddvCFRERqZZUcEVERCJABVdERCQCqlTBNbPfm9kKM8sxs2lmNriM9kMC7XLMbLmZXROpWCPhUPJhZulm5kp4dIlkzOFgZsea2Udmti7wnkaWY54eZva9mWUH5rvLzCwC4YbdoeYjcGvNkvrGyREKOWzM7DYzm2Jmu8ws08w+NrPu5ZivRvaPiuSjpvYPM7vOzGYHcrHLzCaa2WllzBPWflFlCq6ZjQAeBx4AegMTgM/MrE0p7dsBnwba9QYeBJ40s+GRiTi8DjUfQboBhwU9loQzzghJAuYCNwLZZTU2s/rAV8AmoF9gvj8B/xfGGCPpkPIR5GSK941vKz+0iEsHngYGAMcBBcDXZta4tBlqeP9I5xDzEaSm9Y+1wK1AHyAN7/18YGY9S2ockX7hnKsSD2AS8HzItCXAg6W0/wewJGTaC8BEv9+LT/lIBxyQ7HfsYc5LFjCyjDbXAruAOkHT7gTWERiZX1Me5cxHaqBvpPkdbwTykQQUAmeof5Q7H7Wpf2wDrvarX1SJLVwziwP6Al+GvPQl3i+1khxTQvsvgDQzi63cCCOrgvnYZ6qZbTCzb8xsaFgCrPqOAcY754K3/r4AWuB9udRW/zOzzWb2k5md43cwYVIPb8/d9oO0qU39ozz52KfG9g8zizaz8/B+gEwopVnY+0WVKLh4F4+OxtuUD7YJaF7KPM1LaR8TWF51VpF8bMD7hTYc+C2wCPimrOPgNVRpfWPfa7VNFvBH4HfAqcA3wNtmdpGvUYXH48BMYOJB2tSm/lGefNTY/hE4JpsF5ALPAsOcc3NKaR72fhFTGQsR/znnFuEV2X0mmlkq3jGI8b4EJVWCc24L8EjQpKlmlgzcArzhT1SVz8z+BQwCBjnnCv2Ox2/lzUcN7x+LgF5AA+Ac4FUzS3fOzfUjmKqyhbsF7zhDSsj0FGBjKfNsLKV9AdX/tlMVyUdJJgEdKyuoaqS0vrHvNalhfcPMHgXOB45zzi0vo3mN7x+HmI+S1Ij+4ZzLc84tdc5Nc87dhre1f3MpzcPeL6pEwXXO5QHTgBNDXjqR0ve3Tyyl/VTnXH7lRhhZFcxHSXrh7WqubSYCg80sIWjaicB6YKUvEVU9vaghfcPMHudAcVlYjllqdP+oQD5K0osa0j9CRAHxpbwW/n7h96ixoNFgI4A84ErgCLxjD1l4N/IFeA14Lah9O2AP8Fig/ZWB+Yf7/V58ysdNwNl4v0q74Z0m5YDf+v1eKiEXSXhfAL2AvcBdgX+3Cbz+IPBNUPsGeL9IxwLd8Y5p7wL+4Pd78SkflwIXBPpRZ7zjdXnAzX6/l0rIxVOBz/Y4vONs+x5JQW1qTf+oYD5qZP8A/g4Mxhvw1CPwvouAU/zqF74nJSRBv8f7JZGLt4V3bNBrGUBGSPshwPRA+xXANX6/B7/ygXe8ZQneeZnb8I7bnur3e6ikPKTj/XgIfbwSeP0VYGXIPD2AH4AcvF/qd1NDTvk41HwEvlDn4/1A3QVMBS7y+31UUi5KyoMD7glqU2v6R0XyUVP7R+B9rgp8f24GvgZ+42e/0N2CREREIqBKHMMVERGp6VRwRUREIkAFV0REJAJUcEVERCJABVdERCQCVHBFREQiQAVXREQkAlRwRUREIkAFV0REJAJUcEVqCDNramYbzOzuoGk9zSzHzM71MzYRQZd2FKlJzOw3wMd41xmfiXdd3MnOucv8jEtEVHBFahwzeww4E/ge724pvZxzWb4GJSIquCI1jZnFA7PwbtU4wDk3yeeQRAQdwxWpiVKB1ni3ZWvvbygiso+2cEVqEDOLBX4GFgOT8O7neaRzbrWvgYmICq5ITWJmfwcuAHoCO4HPgATgOOdckZ+xidR22qUsUkOY2RDgD8Alzrkdzvs1PRLoCtzqZ2wioi1cERGRiNAWroiISASo4IqIiESACq6IiEgEqOCKiIhEgAquiIhIBKjgioiIRIAKroiISASo4IqIiETA/wO3NOHHeNlR7gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%run LagrangeInterpolation.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Theory\n",
    "Given some function $f\\in C[a,b]$. Choose $n+1$ distinct nodes in\n",
    "$[a,b]$ and let $p_n(x) \\in \\mathbb{P}_n$ satisfy the interpolation\n",
    "condition\n",
    "\n",
    "$$\n",
    "p_n(x_i) = f(x_i), \\qquad i=0,\\dots,n.\n",
    "$$\n",
    "\n",
    "What can be said about the error $e(x)=f(x)-p_n(x)$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The goal of this section is to cover a few theoretical aspects, and to\n",
    "give the answer to the natural question:\n",
    "* If the polynomial is used to approximate a function, can we find an\n",
    "  expression for the error?\n",
    "\n",
    "* How can the error be made as small as possible? \n",
    "\n",
    "Let us start with an numerical experiment, to have a certain feeling\n",
    "of what to expect."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example 1: Interpolation of $\\sin x$\n",
    "\n",
    "<div id=\"interpolation:ex:example_sin_interpolation\"></div>\n",
    "Let $f(x)=\\sin(x)$, $x\\in [0,2\\pi]$. Choose $n+1$ equidistributed\n",
    "nodes, that is $x_i=ih$, $i=0,\\dots,n$, and $h=2\\pi/n$. Calculate the\n",
    "interpolation polynomial by use of the functions `cardinal` and\n",
    "`lagrange`. Plot the error $e_n(x)=f(x)-p_n(x)$ for different values\n",
    "of $n$. Choose $n=4,8,16$ and $32$.  Notice how the error is\n",
    "distributed over the interval, and find the maximum error\n",
    "$\\max_{x\\in[a,b]}|e_n(x)|$ for each $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Define the function\n",
    "\n",
    "\n",
    "def f(x):\n",
    "    return np.sin(x)\n",
    "\n",
    "\n",
    "# Set the interval\n",
    "a, b = 0, 2 * pi  # The interpolation interval\n",
    "x = np.linspace(a, b, 101)  # The 'x-axis'\n",
    "\n",
    "# Set the interpolation points\n",
    "n = 8  # Interpolation points\n",
    "xdata = np.linspace(a, b, n + 1)  # Equidistributed nodes (can be changed)\n",
    "ydata = f(xdata)\n",
    "\n",
    "# Evaluate the interpolation polynomial in the x-values\n",
    "l = cardinal(xdata, x)\n",
    "p = lagrange(ydata, l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Max error is 1.20e-03\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAAEPCAYAAADvfdaQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABXk0lEQVR4nO3dd3hUZfbA8e+ZSS+kkoTQAoQmSA0dpair4q6ruyqKZWHtvayufXV3Xdsqdn+uuotdcO0FFVBB6b2E3nt6b5PMzPv7YwaWkkACmblJ5nyeZx7InTvvPe9NmTNvFWMMSimllFKNzWZ1AEoppZRqmTTJUEoppZRPaJKhlFJKKZ/QJEMppZRSPqFJhlJKKaV8QpMMpZRSSvlEkNUBtDSJiYkmLS2t0corLy8nMjKy0cprbgK9/qD3APQeBHr9Qe9BU67/smXL8owxrWt7TpOMWojI6cDdwEAgFZhkjHmrPq9NS0tj6dKljRbL7NmzGT16dKOV19w0lfpX1FQwJXMK0zZOo8hRRGxoLOO7j2dS70lEBEf49NpN5R5YKdDvQaDXH/QeNOX6i8jOup7T7pLaRQGZwO1ApcWxKItV1FQwYfoEpqydQqGjEIOh0FHIlLVTmDB9AhU1FVaHqJRSTZImGbUwxkw3xjxgjPkYcFsdj7JOTU01ry15hd0lu3G4HIc953A52F2ym38texXj1h8TpZQ6knaXqIDkdrnJ2beN3O2ZVORsw120F3vZPsKqcoiszifSXUqUKSdKKvm8Q1uq7fZay6l2V/PZuv9w27ePUSqRlEskFfZWlIe0xhGRgoluQ1Bce6JTu5HSpQ8xcbV2WyqlVIskunfJsYlIGXDLscZkiMh1wHUAycnJA6dOndpo1y8rKyMqKqrRymtuGqP+jsoyKnO2YS/aQmz5NpKrd9PO7CdCHMd9rdsI/Tq1w4jUeY4Yw+odu+sVS4FpxV57G3JC0yiP7ozEpxOR2AF7UEidrwn0nwHQexDo9Qe9B025/mPGjFlmjMmo7TltyWgExpjXgdcBMjIyTGMOzmnKg338oaH1N243u7euJWvNbMyuBbQpXkkHs/foEwUKaEV2cHvKIjvgikrFFptKWHx7IhLaEtEqkai4RKKi44j97xgKHYV1XjM2LJ6aB5dTWphHWXEeFUW5VOTvprpgL6ZkLyFle4mt3EUb517ipYR4dwlUbvSM9smBahPEtpCuFCVmEN5lBGn9xxKTkHzC96AlCvR7EOj1B70HzbX+mmSoZq8gdz9bF32De8ss0ooW0YECOhzyvMMEszO4E0UxPaFNX2LS+pHSpQ/xca2Jr0f547uPZ8raKUeNyQAItYcyvvt4goNDiE9KJT4ptc5y3C4XWfu2k7NtNeU7VxCcu4bkso20Zx89atbD/vWw/12YC1vtnclNHklUr1/hcoY2/KYopWrldrvZs2cP5eXlVofSIDExMaxfv97v1w0ODiYpKYlWrVqd0Os1yVDNjjGGnRtXsW/hf0nYM5OuNZsYJP/r9iukFTsi+1CdOoSEU0aT1msI3UJO/I16Uu9JzNw1kz2lew5LNELtobSLbsek3pPqVY7NbielfTop7dOB3x08XlyYx46VP1GxeS4xucvoUr2BLq5tdNm3Dfa9QycTyooVGbi6nUe3kRfTKiHphOuiVKDLy8tDROjevTs2W/OZ+1BaWkp0dLRfr2mMobKykr17Pa3BJ5JoaJJRCxGJAtK9X9qADiLSDygwxuyyLLAAt23dUrLmvUfqvpmkmT2keY9XE8Sm0FMpb3c6yf3Po+Mpg+jfiH88IoIj+GDcBz5bJyMmLpG+Yy6GMRcDUFVZzurFMyhf9z3JuQvo7N5B//J5sGIezuUPkRnWl4r08+g+9kpiElIao4pKBYyioiLS0tKaVYJhFREhIiKCtm3bsm/fPk0yGlEG8NMhX//V+3gbmGhFQIGqoiiHBe/8haQdX9LFvZ3O3uPFRLI5ZiRBvX5D9+G/pXfUiTXl1VdEcAQ397+Zm/vf7NPrAISFR9Jn1IUw6kIAvv3iI+KrthO5bTo9qlbT27EC1q6gOvMJVkYOxtXr9/QacylhEf79lKNUc+RyuQgODrY6jGYlPDycmpqaE3qtJhm1MMbMBuqeTqB8yuGoJPPHaQStepdzKpdh83aFlBDJhvixRPa/mO5DziHjJLpAmpPwmCSG/PYS4F6K8rLZ9Mt/Cdv4Gb0ql9GvYgEsWUDJkr+wKvEcEk67hvS+I6wOWakmTY4xW0wd7WTulyYZqsnYu20du2a8TPesrxhICeDpCsmMGoGt7yWcctrvGRwWbnGU1opNTGbwhbcAt5CbtYutP75L/LbP6ebcxJC8T+GzT9nyZRcKelxG73OvIyIqxuqQlVIBTJMMZSm3y82aXz7DLHqdPhWLaOtttdhh60B21/EUtTqVs8/7rcVRNk2tUzrQesKDwINsy1xEzpw36ZE7nXTXVlj7GCVrJ7Mw+Xw6nH07qZ1PsTpcpVQA0iRD+VRdG4td1uUSNs94l5S1/6av2QN4Wi1Wx5xBq9NvoGv/0aTZbMyePdvS+JuLzr2H0Ln3EKoqy1k66z2iVk2hh3M9Q7On4n57GisjhxI66k56DDqLSleVZZu9KaVOjNvt5vrrr+fjjz+moKCAn376ib59+9KjRw/mz59Ply5djluGw+Gga9eufPrpp2Rk1Lp2VqPTJEP5zIGNxQ6d+lnoKOTfq17nu6WvMG3fPiKMIYd4tqVdSrdzbyYjuZ3FUTdvYeGRZPzmevjN9Wxe+QuFP71Mv6JZnrEb3y5gxawe3J8WQZ4pPex7MmXtFGbumskH4z7QREOpJuj7779nypQpzJ49m86dOxMfH8+DDz7IuHHj6pVgAISGhnLPPfdw77338sMPP/g4Yg+dw6N8ZkrmlKPWlgCoETf7g4TJsR1YOvBp4h/YwNCJTxCvCUaj6trvNAbfOY2Sm1axoN3VFBPJ/Mh95FXn1LrZ257SPUzJnGJRtEqpY9m2bRtt2rRh+PDhpKSk4HQ6efPNN7n66qsbVM7ll1/O3LlzWbt2rY8iPZwmGcpnpm2cVusqmQAOm40ZybFk/OZ6ggJklohVEpPbMeyayQT9aR3vxSXjqGN9AIfLwbSN0/wcnVLqeCZOnMj999/Prl27EBHS0tKYPn06IsKIEf+bTfb3v/+dlJQUcnJyDh677LLLGDBgANXV1QDEx8czYsQIPvzwQ7/Ert0lyicKsvdQWFV4zInARY4iv8WjIDI6lnKqj3mOfk9UoEm77xtLrrvjyfPqfe4LL7xASkoK77//PkuWLMFut/PYY48xcODAw6aXPvDAA8yYMYM//vGPfP3117zzzjt88cUXLF++nJCQ/23COHjwYObMmdOo9amLtmSoRlVSlMfCN+8g7NUBxLldxzw3NjTWP0Gpg453z1u53Cz//l2M2+2fgJRSxxUTE0N0dDR2u52UlBRat27Nzp07SU09fK8ku93Oe++9x9y5c/nzn//MLbfcwrPPPkuPHj0OOy81NZUdO3b4JXZtyVCNwlFZxsqPn6bn1jcZSjkIjKqI4ZtWDmrM0SvFHdhYTPnXsTZ7C3EbLispZsCuW9i05EVqxjxKrxH1/7SlVHPUkBaFpqSyspLk5OSjjnfs2JEXXniBiRMnct5553HjjTcedU54eDiVlZX+CFNbMtTJcbtcLPvy/yh6ui9Dtr5AK8rJDOnL+nGfcP/1P9AhpgOh9sPHXDR0YzHVeCb1nkS76Ha1f09iu9A75TryiKWbcxO9Zk5g5VNns3P9MouiVUrVJTExkcLCwlqf+/nnn7Hb7ezevRuH4+gPFAUFBbRu3drXIQKaZKiTsH7eV2x/PIOBy+8j2eSxzZbGqtH/odf9c+g5+MyDG4tN6jWJuNA4BCEuNI5JvSbpVEmLHOt78uGvP2TUpQ8TcfdqFna8gXITRr/KhbSbegaLX7yCvKzdVoevlPLq378/69atO+r4p59+yvvvv8+PP/5IcXEx999//1HnZGZmMmDAAH+Eqd0lquH2b1tL9if30K98HgDZJLCj711k/OYG7EGH/0j5c2MxVT/H+55ERMUwdNJT5GffQuZHDzMw7wsGF3xF2f/NYmH6dfS/+D5CwzRBVMpKZ599Nvfeey/5+fkkJCQAsHfvXq699loef/xxTj/9dN59913Gjh3LuHHjOPPMMw++9pdffuHvf/+7X+LUlgxVbxWlBSx+/WYS3j6dfuXzKDehzO94I63uWc2QC285KsFQzVtCcnuG3PoW+y7/iVXhQ4iSSoZufYHcp/qzcuYHOjhUKQudeuqpDB48mKlTpwJgjGHixIn079+fO++8E4DTTjuN++67jz/84Q/k5+cDsGDBAoqLi7nooov8EqcmGeq4jNvNyq//RcWz/Rm87z1CxMniVmdTdt0ihk96kvDIKKtDVD7UoVs/+t47g9Wj/81OWzvamSz6zbuRNU//it1b1lgdnlIB4bbbbjtqRsgjjzzCiy++iMvlQkSYOXMms2bNOmxa69///nf27t17sLVj8uTJ3HPPPYSH+2ezSU0y1DHtWreYjU+eRr+lfyaRIjYE9WD9r79g8F0fkdy2k9XhKT/qM/oiUu9bzqLuf6aECPpULSH53dEseOMOKstLrQ5PqYBzzjnncPPNN7Nnz556ne9wOOjTp8/Blg5/0CRD1aqytIglr11P6rSz6VGdSQGtWHjq3+h6/3x6Zoy2OjxlkeCQUIZc9iA1Ny5mSey5hIiTYXunUPxPTxeKUsq/brvtNjp27Fivc0NDQ3n44Yf91ooBmmSoIxlD5sx3KH12AIOypiIY5if8Htutyxn6+9ux2+1WR6iagITk9gy6Yyobxn3CVntnUsil37wbWfH0OLJ2b7E6PKVUE6Ej9dRBubs3kf3hLfSuWATABns3+PVkhvc/zeLIVFPVY/CZuAYsZtF/n6L3hpfoXzGPijeHs7DbzWRccj9BwSHHL0Qp1WJpS4bC7axh6Yd/I/LNkfSuWESJiWBu9wdIv28+PTTBUMdhDwpmyGUPUX7dQpZHnU6EOBi6eTI7nhzKllXzrA5PKWUhTTIC3K61C9j+5DAyNj5LhDhYFDmGiusWMvKyewkKDrY6PNWMJLXtxIC7v2L16a+TRWvSXVtJ+/TXLPjXzTowVKkApUlGgKquLGfpG7eS+tE4ujg3s59Elg7/F4Pv/oyUtvUbRKRUbfqMHU/0n5ayMOkSbBiG7X+PgmcGkjn3S6tDU0r5mSYZAWjr0hnk/DODjL3vYMMwN/FiIu5YSsavLj1sfrVSJyoyOpahN73BlvM/Y7stjbYmm96zrmTRi1dSUpRvdXhKKT/RJCOAVJUVsez/rqbL1xfTzr2P7dKezHM/ZuQtbxITG2d1eKoF6jZwDO3uW8zCjjdQbewMKfiSquczWPXjVKtDU0r5gSYZAWLzwq8penYQA7M/psbY+bnNJJLvWUSfoWce/8VKnYTgkFCGTnqKfZfOYFNQN5IooO/P17PkuUsoLsi1OjylWozNmzeTnJxMcXFxvc7PycmhdevW9V7M60RoktHCVZUVsfzViXT97nJSTA6bbF3YfMFXnH7980RERFodngogaT0z6HLfAham30mVCWZQ8fc4XhysrRpKNZIHHniAm266iZiYmHqdn5SUxFVXXcUjjzzis5g0yWjBtiz6hqJnBzEg5zOqjZ057a6jw5/nc0r/EVaHpgKUPSiIoVc8Su4VP7Ah+JRDWjXGU1yYZ3V4SjVbu3fv5vPPP2fSpEkNet2kSZN4//33KSgo8ElcmmS0QI6KEpb93zWkfzvhYOvF1gu/ZtQ1/yQsLMzq8JSifde+dL33FxZ2vcvbqvEdjhcGsXr2J1aHplSTNG7cOG644QZuv/124uLiiIuL45577sHt3Q152rRp9O7dmw4dOhx8zdVXX02vXr2orKwEwOVycdppp/HrX//64Dm9e/cmNTWVTz/91Cdx64qfLcy25T8Q+vXNDHTvp8bYWdDuagZf+XdNLlSTYw8KYujlj7B702+p+Oh6ujs3kDT7jyxa/Tm9Jr5IVCsdjKz84NH6dS00/nXrN27iUO+//z4TJ05kwYIFrF69mmuvvZY2bdpw11138csvv5CRkXHY+S+++CL9+/fn7rvv5pVXXuEf//gHmzdvZvXq1YedN3jwYObMmcM111xzUlWqjSYZLYTTUcmqd++h/+73sIlhq3Sk6jevcPpAXbFTNW3tu/XDdd88FnzwVwZufZUhBV+y77mF7DzrBXoNH2d1eEo1GW3atOHFF19EROjRowebNm1i8uTJ3HXXXezcuZN+/foddn5kZCTvv/8+I0aMICEhgSeeeIIvv/ySpKSkw85LTU1lyZIlPolZk4wWYPe6hbg+uY6Brp24EH5OvpKMPzylAztVs2EPCmLYVX9n+9pf4/r0RtJdW0n5fgILV11Kv4nPWh2easlOoEXBKkOHDj1sLaNhw4bx8MMPU1JSQmVlZa0t1oMGDeLBBx/k0Ucf5aabbuLcc8896pzw8PCDXSqNTcdkNGNup5Nl7z1E8rRxpLl2skvakHnOfzn9xpc1wVDNUqdeQ+h47wIWtrsaN8LQ7A/J+ucQSvdvsjo0pZq0xMRECgsLjzpujGHu3LnY7Xa2bt2KMeaocwoKCmjdurVP4tIko5nK3rGeLU+fxsAtLxEiLubGXUDMHQvpO+wsq0NT6qQEh4Qy9JrJbPvtZ+yWVNLcuzlnw30smHIvzppqq8NTyjKLFi06LElYuHAhqamptGrViv79+7Nu3bqjXjN58mSWL1/Ozz//zMKFC3nppZeOOiczM5MBAwb4JGZNMpoZ43az4vMXiH5rFN2q15FDPEtHvsHI298mJibW6vCUajTdBowm8e7FLGx9McHiYtjO19j61Eh2b15ldWhKWWLfvn3ccccdbNy4kY8//ph//vOf3HnnnQCcffbZLFy4EKfTefD8VatW8eCDD/LGG28wfPhwXn31Ve69917Wrl178JyKigqWLVvGOeec45OYNcloRopz97Lm2fPov/IvROBgceRobDctIOPMS6wOTSmfCI+MZujNb/JF2iNkk0B350YS3juLRR89jfFO3VMqUFx++eW4XC6GDBnCtddey9VXX30wyRg3bhzh4eF8//33AFRVVXH55ZczYcIEfve73wEwYcIELrroIiZMmIDD4QDgiy++oEOHDpx2mm8mCejAz2MQkZuAe4A2wFrgDmPML76+bkVNBVMypzBt4zQKHYXETY1jTGhvrl39NX1MMSUmgrX9/8LQ869HbJonqpYvJm0AYb9dzNIpN5BRMpMh6/7B6qdnEDfhZb7M/4FpG6dR5CgiNjSW8d3HM6n3JCKCI6wOW6lGFRQUxMsvv8zLL7981HN2u50HH3yQyZMnc9555xEWFkZmZubB511uF/lV+Tz4woO43C62l28n3hXP5Ocm85e//MV3Mfus5GZORMYDLwA3AXO9/34rIqcYY3b56roVNRVMmD6BPaV7cLg8mWaho5BvKuewKjWch/I60uay/zAsrauvQlCqSYqJSyTjro9Z/u1/6LzoYdIdS7n0+wvZExJKDS7A87syZe0UZu6ayQfjPtBEQwWUa6+9loKCAoqLiw9bWtzldrG9ZDvVruqDYzpcbhcbd21kzHljuGS871rD9WNw3e4C3jLGvGGMWW+MuRXYD9zoy4tOyZxyWIJxgMNmY1dIGAvPHk9bTTBUABtw7h+puW4eTyZ0Z1+Q7WCCcYDD5WBP6R6mZE6xKEKlrGG323nggQeO2rskvyr/sATjgPjEeCbeMpECh2+WFAdNMmolIiHAQGDGEU/NAIb78trTNk47KsE4oAY3H236yJeXV6pZaJ2axk+JITjq6C50uBxM2zjNz1Ep5TvTp0+vtZukPgqqCmqdugqeKa4FVb5LMrS7pHaJgB3IPuJ4NnDU3ugich1wHUBycjKzZ88+4QsXOo6e53yoIkfRSZXf3JSVlQVUfWuj96D2e1DkKDrma1rS74r+DDTePYiJiaG0tPTkA/Izl8t1wnG73K7jPn+8squqqk7o/muS0QiMMa8DrwNkZGSY0aNHn3BZcVPjjploxIbGcjLlNzezZ88OqPrWRu9B7ffgeL8r0UGRLea+6c9A492D9evXExUVddjKmc1BaWkp0dHRDXqN2+2mPH8vdjG4qLu+dpv9mGUbYwgLC6N///4Nuj5od0ld8gAXkHzE8WQgy5cXHt99PKH20FqfC7WHMr77eF9eXqlm45i/K24343P2seCdh3Edsm6AUna7nZqaGqvD8LmqygqqszYQXZNHvMtdZ4ohIsSHxR+zrMrKSoKDg08oDk0yamGMqQaWAUcun3kWMN+X157UexLtotsd9ccz1B5Ku+h2TOo9yZeXV6rZqOt3JcQWQoKJ4JqSYoZte5FNT53O3m3rLYpSNTWxsbFkZ2cf3CK9pTHGUJq/n5CCTYThoIYgoiLbExIUelTrjYgQYg8hISyhzrIqKirYu3fvUZuq1Zd2l9RtMvCuiCwG5gE3AKnAa768aERwBB+M++DgOhk691+p2h3vd2XL3Om0mXM3PWvWUv72aBb3uY9BF96ua8sEuMTERPbs2cPGjRutDqVBqqqqat0A7VBOpxN3eT4hxjN5oNoWTlBUArbiLNzGTXlNOeU15biNG5vYiAyOJDI4kk3H2BsoODiY5ORkWrVqdUJxa5JRB2PMNBFJAB7CsxhXJjDOGLPT19eOCI7g5v43c3P/m7UvVqljOPR35Uh9xlxE0akjWf7WdQwom8PgNY+yasu3tP3DmySmdLAgWtUU2Gw2OnRoft//2bNn1zkmwhjDgi/foNfyR4mRcoqIYs+Ix+l71h/8HOXRNKU/BmPMq8aYNGNMqDFmoDHmZ6tjUkrVX2xiCv3v+pylA5+mhEj6Vi7C/tpwVnyra2ioliE3Zz+LnrmQ4SvuIUbKWRsxGHPjAno3gQQDNMlQSrVwYrOR8ZvrqbzmF9aEDiCOUvovuoOlk39PcUGu1eEpdcIWz/wI8+owhpb/RCWhrOz7CKfc/T1xyU2npUaTDKVUQEhu14Xe9/7Aop4PUmFCySiZhePFwaz66b9Wh6ZUgxQWFjD3uSsYPO9akihkU8gplE6aTb8L72pyY46aVjRKKeVDYrMxZPyfyb/yR9YHn0ISBfSdcw2LX7ic0mLfrXqoVGNZNucrKl4Ywsjir6g2QSzrdgdd7/2FpI6nWB1arTTJUEoFnPbpvel27y8s6HIH1SaIwYVfU/bcYDJ/+cLq0JSqVVVVJb+8+EcG/nQFbclhW1AX8i6fwcAJf0XsTXcOhyYZSqmAZA8KYtiVf2Xfpd+z2Z5OG3Lp/cNVLHrpD5SVHHt5f6X8afkv33DKgts5reATnMbGsrTr6HjvQlK7DbQ6tOPSJEMpFdDSemaQdu98FnS8kWpjZ0j+55Q+N4i12qqhLFZcUswvL15Nv1mX00Gy2WFPY/8l0xk48Z/Yg0OsDq9eNMlQSgW84JBQhk16kj0Xf8sWexfamFx6/XAVi1+8gpKifKvDUwFo2ewvKJ08iNMKPsaNMCvm97S/bzHtew2zOrQG0SRDKaW8OvceQsd7F3hbNYIYXPAVVc9nsPpH3TZe+UdhQR7znr+CgbOvoh0HWi++Iaj/VdiDa9+rpynTJEMppQ5xoFVj36Uz2BjUnSQK6PPzdSydfBGFufutDk+1UMYYFn33PtUvDmJE0VdUGzvLOt/kbb0YbnV4J0yTDKWUqkVaz4Gk3zefhel3UWlCyCiZCa8MZumXr2Fa6OZayhr79+5kyT9/y5CFN5FMAVuCu5N3+SwGXvVEs2y9OJQmGUopVQd7UBBDr3iE/Ct/IjOkH3GUkLH8XtY8fRb7djSvDbZU0+NyuZn30WQiXh/G4Io5VBDKsh730Pne+aR2G2B1eI1CkwyllDqOdum96XXfTyzu83eKiaRP1VJip5zGwnf/Qk21w+rwVDO0ee1S1j15OiPW/fXgniMVV89l4KUPYQtquuteNJQmGUopVQ9iszH4d7dRfcNClkaPIUIcDN36AnueHMyGxTOtDk81E2VlJcx97TY6fvQrTq1ZQwGtWD3kGXrdM4PE9t2sDq/RaZKhlFIN0DqlAxl/+pxVo/7NXkmmk3sHPaZfxKIXr6QoL8vq8FQTZYxh6cypFD+TwcistwkRF8sSf0vw7cvoc+61IGJ1iD6hSYZSSp2AvmMuIv7u5SxoO8mziFfBl5iXM1j88WTcLpfV4akmZNeWtax4+lwy5l1PW++01K3nf8rAW94hOi7J6vB8SpMMpZQ6QeGRUQy79nn2XzqTtSF9iaOUwZl/ZcsTQ9m0fLbV4SmLVZSXMu/NP5H87igGVC6gjHCWdP8T7e9bQpcBZ1gdnl+0nNElSillkY49B2K6z2bZt/+h/ZJ/0M25CfcXF7Dk53PpdOlTJKZ0sDpE5UfG7Wbpt1Not+QJRpALAstifkWny55lUID9LGhLhlJKNQKx2Rh43jVE/GkFC9pcgRMbg4qmE/5/g1jwzsM4qiqsDlH5waYVc1n3xEgGLbmLNuSyzZ7GxnM/YuCd/yU+wBIM0CRDKaUaVVR0LMOuf4XsK+ewMmIYkVLFsG0vkvdUP5ZNn6ILebVQWXu2sej5y0j//Nf0qllLAa1YeuojdLx/Gd2HnG11eJbRJEMppXygffqp9Pvzd6we8xY7bO1pa7IZuPgONj0+lA2Lvrc6PNVISooLmP/67cS8MYQhRdNxYmNRygSC71hJxu/vwt6C1rw4EYFde6WU8rE+oy7EOfw8Fn3+Il3WvkR350b49hJW/DyC+N/8nY49BlodojoBVZUVrPjseXpsepXhlILAiqhRJF/4BEO69LI6vCZDkwyllPKxoOAQhlx8N2VnX82C/z5G313v0r98Hq4Pz2BJ7Nm0veCvxLbrwJTMKUzbOI0iRxGxobGM7z6eSb0nWR1+wKqoqTjqe3Jx14vou91Bt8zXGEYeABuCT8F2zj/oP3CsxRE3PZpkKKWUn0S1imPY1c+Su+9W1nzyCAPyvmJQ8XcUvT2LCzt0Ji/YRbW7BoBCRyFT1k5h5q6Z3Bh9o8WRB56KmgomTJ/AntI9OFyepeMLHYW8tfp12jtr+EDy2S4dKRl2L33OuAyx6eiD2uhdUUopP2udmsaQW98m5w/zWBJzNu/FRJIvFQcTjAMcLgd7Svcwq2SWRZEGrimZUw5LMA6otgm7g4L5R4/f0fGB5fQ963JNMI5B74xSSlmkbeeeDLrzIz5Maoujjjcqh8vB3NK5fo5MTd0w9agE44Bqm/CL2dKiNjLzFb1DSillsVJn+TGfL3eX+SkSVVKUx7ovX6DYXXjM/USKHEX+C6oZ05YMpZSyWGxo7LGfd7lY+/hIVsx4D5fT6Z+gAsy+7RtY9NoN2J7rzdBtLxJ7nPVMjvc9Ux6aZCillMXGdx9PqD201ueCjY0Li6voVb2G/vNvJvsfPVn43iMU52f7OcqWx7jdZP7yOSufPoeUt4YyJOtDoqSSzJB+jI4bU+f3JNQeyvju4/0cbfOk3SVKKWWxSb0nMXPXzKMGGobaQ2kX3Y6uSVexsHAN7be8T1uTTeqW56l88VUWx51Jq+F/pHvGGTr4sAEKc/exacabtNk6jd7uPQBUE8SqmLHEjr6V3v1Pp3NNBauPmF0C//ue6NTi+tEkQymlLBYRHMEH4z6oc52MxfMWM3TcI7icD7JyzsfYlrxOn6plDC6aDtOns/O79uzvcjFdz7yahOR2VlenSXI5a1g/7wuql7xD79K5DBEXADnEs6XDJXQ97xYGJbc/eP7xvicRwRFWVaVZ0SRDKaWagIjgCG7ufzM397+5znPsQUH0O+NSOONSdm9Zw54fXqfr/i/p6N5Nx82TcW56nlXhA6nu+Tt6jp1AVHSs/yrQBBm3my0r51Cw8APSc2bQmyIAXAgrw4fi7n8lvUddzPDQ2rtF6vM9UcemSYZSSjVD7dNPpX36S9RUP8OK2f9FVrxLr4ol9K1aAiuWULn8EZZFDcX0+DVdT7uImNgEq0P2C7fLxZYVsylc9intsmbR1WQdfG63pLK7/fl0/tX19GvX2cIoA4cmGUop1YwFh4TS/1dXwK+uoCh3Pxt/epdWmz+nZ81aBpb/DMt+pnrp/awO70dFx7NIzTiP9l16t6gxHOUlBWxZ9C2ODTPolD+HbhQefC6PWDYnnU380Cvo1m8k7VtQvZsDTTKUUqqFiG3dhiGX/Bn4M9m7NrNj3kdEbf+OHo419KlaBhuXwcYn2SdJ7I4bRlD6aDr0O4PWqR2tDr1BHFXlbF89l6J1s4nZ9zPpjvX09Y6xANhPIjtbjyWq/wX0GPQrhgUHWxhtYNMk4wgich1wGdAfiAE6GWN2WBqUUko1UHKHriR3eBB4kKLc/Wye9zG2rT/QpXQJqSaH1IIvYPEXsBj2SAr7WvXDtB1EXHoGHXpkEBYRdULXrW1TsZMZLGncbvbv3EjJ5l9YuOkTYnOX0rl6Ez3kf+uFuBDWB/ekMOU0Wg88n/Q+I2ijLRZNgiYZR4sAZgBfAM9ZHItSSp202NZtGHTBrcCtuJxONq2eS8Gqb4nMWkznqnW0I4t2xd9B8XewDpxf2Nhub09+VFdqYrsQktyN2PankJzWk6hWcXVep65NxQ5s9PbBuA/qTDSc1Q5y920nb0cmFfs3IvlbiCzZQvvqraRSzvmHnOtG2GbrSE5cf4LTx5I+ZBw941s34h1TjUWTjCMYY54HEJEMi0NRSqlGZw8KotuA0TBgNADOmmq2rFtM/rrZ2PevJLFsI+1du+nk3kmnkp1QMgt2AUs8ry814eTbEykJScIRmogrJAYTFouEx/ClbGZX9Q5qcB12TYfLwa6iHTz20bWc7+wIVUXYHCUEOwqJdOQQ58wlwRTRRgxtaok5nxi22ztRndyP8M7D6NR/LJ0TktChm02fJhlKKRXAgoJDSO87kvS+Iw8eq6ooZcu6xZTsysSVu5mw4q3EVe0i2ZVNtFQS7d4NVbuh6vCy/tyhLTV2e63XqcHFXMcKHt/1da3PuxFyiCc3pB1lUWmY+HTC2nSnTY8hJLXpyJqff2b06NGNVW3lJ2KMsTqGJsnbkrGEeozJ8I7juA4gOTl54NSpUxstjrKyMqKiTqxvtCUI9PqD3gPQe9BU6m/cbqory6guy8NdnotUFWGvKcfuLCPYWc4dyeug7j3FEAPP5PanJigSV3AkJqQVRCQQFJVIWFQ89qC6B2g2lXtglaZc/zFjxiwzxtTa+h8QLRki8hjw4HFOG2OMmX0i5RtjXgdeB8jIyDCNmW3Pnj07oLP3QK8/6D0AvQfNpf5xU0+n0FFY5/OxYXH86p53Tqjs5nIPfKW51j8gkgzgeeC945yzyw9xKKVUizW++3imrJ1y2F4fB+imYoEpIJIMY0wekGd1HEop1ZIdb6M33VQs8OiYjCOISAqQApwCvA+cB+wDdhljCurx+lxgZyOGlEhgJ0iBXn/QewB6D5pP/QVbUFxQSlBUUGtsBOHG6Sxz5joLnVkY3CdRcvO5B77RlOvf0RhT6xxiTTKOICKPAo/U8tQkY8xb/o0GRGRpXQNqAkGg1x/0HoDeg0CvP+g9aK711yXRjmCMedQYI7U83rI6NqWUUqo50SRDKaWUUj6hSUbT97rVAVgs0OsPeg9A70Gg1x/0HjTL+uuYDKWUUkr5hLZkKKWUUsonNMlQSimllE9okqGUUkopn9Ako4kSkZtEZLuIVInIMhE5zeqY/ElETheRL0Vkr4gYEZlodUz+JCL3i8gSESkRkVwR+UpEelsdl7+IyM0istpb/xIRWSAi51kdl1W8Pw9GRF62OhZ/EZFHvXU+9JFldVz+JiJtRORt79+BKhFZJyKjrI6rvjTJaIJEZDzwAvA40B+YD3wrIh0sDcy/ooBM4Hag0uJYrDAaeBUYDowFnMAsEYm3Mig/2gPcCwwAMoAfgc9FpI+lUVlARIbi2eV5tdWxWGAj0OaQx6nWhuNfIhILzMOzt+15QE/gViDHwrAaRGeXNEEisghYbYy59pBjm4GPjTH3WxeZNUSkDLglkBdEE5EooBi4wBjzldXxWEFECoD7jTH/sjoWfxGRGGA5cA2elYgzjTG3WBuVf3hXX77IGBMwLXhHEpHHgVHGmBFWx3KitCWjiRGREGAgMOOIp2bg+VSrAlM0nt/XuvfRbqFExC4il+Jp3ZpvdTx+9jqeDxc/WR2IRTqLyD5v1/FUEelsdUB+dgGwSESmiUiOiKwUkVtERKwOrL40yWh6EgE7kH3E8Ww8G7epwPQCsBJYYHEcfiMip3pbsRzAa8CFxpg1FoflNyJyLZAOPGR1LBZZBEwEzgGuxfP3b76IJFgZlJ91Bm4CtgFn4/k78CRws5VBNURAbPWuVHMmIpOBkcBIY4zL6nj8aCPQD4gBLgLeFpHRxphMS6PyAxHpjmdM1khjTI3V8VjBGPPtoV+LyEI8b7Z/ACZbEpT/2YClh3STrxCRrniSjGYxCFhbMpqePMAFJB9xPBkIuJHVgU5EngMuA8YaY7ZZHY8/GWOqjTFbjDHLvH9kVwJ3WhyWvwzD06q5VkScIuIERgE3eb8OtTY8/zPGlAFrga5Wx+JH+4F1RxxbDzSbSQCaZDQxxphqYBlw1hFPnUXg9UcHNBF5gf8lGBusjqcJsAGB8ub6OZ6ZFP0OeSwFpnr/X21JVBYSkTCgB5433kAxD+h+xLFuwE4LYjkh2l3SNE0G3hWRxXh+yG4AUvH0SwcE72yKdO+XNqCDiPQDCowxuywLzE9E5BXgSjwDvwpF5MB4nDLvJ7oWTUSeBL4BduMZ9DoBz7TegFgrwxhTBBQdekxEyvH8/Lf47iIAEXkG+ArYBSQBDwORwNtWxuVnz+EZh/IgMA3Pkga3AQ9YGlUD6BTWJkpEbgL+jGdueCZwpzHmZ2uj8h8RGQ3UNqL+bWPMRL8GYwERqesX86/GmEf9GYsVROQtYAyewX7FeNaI+Kcx5nsr47KSiMwmsKawTgVOx9NtlAssBB42xhzZfdCieRehexxPi8YuPGMxXjLN5M1bkwyllFJK+YR2lzSyxMREk5aW1mjllZeXExkZ2WjlNTeBXn/QewB6DwK9/qD3oCnXf9myZXnGmNa1PadJRiNLS0tj6dKljVbe7NmzGT16dKOV19wEev1B7wHoPQj0+oPeg6ZcfxGpcyDqcZMM79rpv8MzfSoNCMfTP7Yc+NYYozMelFJKKXWUOqewikiqiLyJZ7rQg0AInilUM/BMnxkFzPTuCDfeH8EqpZRSgWpTdin/mrOVeVvyqHG5rQ6nXo7VkrESz1ShDGPM2tpOEJFwPFPs7hKR9saYZxo9QqWUUiqAFVVU8+46B7Nn/ILL7ZmsERsRzBk9khl3agpjeyTRVLczOVaS0csYk3usFxtjKoEPgQ9FpNZBH0oppZRqOKfLzYeLd/HszE0UVTixCYw7NYWNWaVszS3nk+V7+GT5Hm4Zk87dZx+5ZlfTUGeScbwE4wARCTbG1NT3fKWUUkodm9ttuGPaSr5e7VngtGe8jeeuGkGPlFYAbMkp5ZvVWTz/wyZem7OV3/RNpXtKtJUh16pey4qLyJe17XwnIj2BxY0elVJKKRXA/jljI1+v3k90aBCvXTGAPw8KO5hgAKQnRXP7mV25YkhHnG7DA5+twe1ueute1XfvkjhgjYj86sABEbkFzx4bq30R2CHXuUlEtotIlYgsE5HTjnP+KO95VSKyTURuaGiZInKdiPwkIkUiYkQkrZGrpZRSStXqw8W7+L/ZW7HbhFevGMA5vdvUOebinnO60zo6lGU7C5m2dLefIz2++iYZo/Dsm/GViLwkItOBvwNXG2P+4KvgvLNWXsCzpGp/PBuEfSsite5AJyKdgOne8/oDTwAvicjvG1hmBJ5ZNI82cpWUUkqpOv2yOZeHPvdsT/OPC3pzWtdjD3dsFRbMX359CgBPfruBvDKHz2NsiHolGcYYtzHmb8CTePaxPxM41xjzoS+DA+4C3jLGvGGMWW+MuRXPlNob6zj/BmCfMeZW7/lv4Jkhc3dDyjTGPG+MeQKY64tKKaWUUkfamFXKTe8tx+U23Di6C5cOrt+O7r/u04bTuiZSXFnDP75Z7+MoG6a+YzJCReQl4F48n+7n4GnV+K2vAhOREGAgnhaFQ80AhtfxsmG1nP89kCEiwSdYplJKKeVTVTUubvlgOaUOJ+f1acM9v6r/bBER4bELehMaZOOzFXuZvyXPh5E2TH2XFV/mPXekMWYpgIj8CZgqIu8ZY671QWyJgB3IPuJ4Np6WlNqkALNqOT/IW56cQJnHJSLXAdcBJCcnM3v27BMt6ihlZWWNWl5zE+j1B70HoPcg0OsPLf8efLSxms05NaRECOcnF/Pzz3MOe74+9R+XZuezLW6e/Wop1f3DfBht/dU3yVgI3GaMqThwwBjzrIjMAt7zSWTNiDHmdeB1gIyMDNOY68s35fXq/SHQ6w96D0DvQaDXH1r2PVixq5Dvvp+PTeCVPwxjYMe4o86pT/3T+1bw2VM/sbFIGD7ydEKC6jvs0nfqOybjmkMTjEOOrwIyGj0qjzzABSQfcTwZyKrjNVl1nO/0lnciZSqllFI+UVXj4u7/rsJt4NrTOteaYNRXu7gIuidHU+ZwsmRHQSNGeeKOtXdJvVb1MMY4GnJ+fRljqvF005x1xFNn4ZkRUpsFdZy/1Ltg2ImUqZRSSvnEczM3sTW3nC6tI7nzrG4nXd7YnkkA/LA+56TLagzHasnYLCIPiUi7uk4QEZuInCsiM/HMOmlsk4GJInKNiPQUkReAVDzTaRGRd0TknUPOfw1oKyLPe8+/BpgIPFPfMr3lpohIP+DAd/wUEeknIvE+qKNSSqkAtGxnIW/8sg2bwDMX9yUs2H7SZZ7Rw5Nk/LjhyKGH1jjWmIzTgH8A20RkDZ4dWPcBVXgW5zoFGApU4llz4o3GDs4YM8270uhDQBsgExhnjDmwd32HI87fLiLjgOfwTEndh2csyScNKBM8U2EfOeTrb7z/TgLeaqTqKaWUClDVTjf3fbIat4HrR3Wmf4cT7yY5VP8OccRGBLMjv4JtuWV0bh3VKOWeqGPtXbIZuERE2gOX4Ek6BgPheMY2rMAz2HG6McZne84aY14FXq3judG1HJsDDDjRMr3PP4ouxKWUUspH/jNvO5tzykhLiODOM0++m+QAu00Y3a01n6/cx48bcixPMo478NMYs9sY86wx5gJjTH9jTA9jzEjvgldf+zLBUEoppVqavUWVvDBrMwB/+23vRukmOdSYg10m1o/LaPD8FhGJEhFrUyOllFKqmfrrl2uprHFx3qltOL3bsZcNPxGjurXGbhMWby+gpKqm0ctviHonGSJyh4jsAoqBYhHZLSJ3Sl27tiillFLqMD+sz2bGumwiQ+w87N1zpLHFRoQwsEMcTrdh7mZrV/+s77LiT+MZo/AvPNM9z8IzG+MvwFO+Ck4ppZRqKSqrXTzy5VoA7jyrGykxvluVs6lMZa1vS8Y1wDXGmH8YY370Pv4BXAtc7bvwlFJKqZbhlZ+2sKewkh4p0UwcnubTax2Yyjp7Yw5ut/HptY6lIWMyVtdxzPp1S5VSSqkmbGtuGf/6eSsA/7iwN0F23751pidF0S4unPzyalbtKfLptY6lvrV8h9oX27oReLfxwlFKKaVaFmMMf/kikxqXYXxGewZ29P26jiJyyMJc1nWZ1HeDtFBggoicjWezNIAheFbKfF9EXjxwojHmtsYNUSmllGq+vly1j3lb8omNCObec3v47brD0xN5e8FOVu0p9ts1j1TfJKMHsNz7/47ef7O8j56HnGddx49SSinVxJRU1fDYN+sBuO+cHsRHhvjt2j1SPFuKbcwq8ds1j1SvJMMYM8bXgSillFItzeQZm8gtdTCgQyyXZLT367Xbx0UQHmwnu8RBUUU1sRH+S3AOOJHFuC4TkUhfBKOUUkq1FJl7i3lnwQ5sAo9dcCo2m3+XlbLZhG7JnrUzN2aV+vXaB2M4gdf8C0hu7ECUUkqplsLlNjz4eSZuAxOHd+KU1FaWxNHd22WyKbv5JBm6wqdSSil1DFPmbWfV7iJSWoVx51ldLYujW7InydjQjFoylFJKKVWHnfnlPDNjI+BZEyM6LNiyWHqkeFpQmlNLxrnAvsYORCmllGrujDHc98kaqmrc/LZfKmf0tHZ0QbcUz5iMDVmlGOP/CaANTjKMMXONMVW+CEYppZRqzj5cvJsF2/KJjwzhkd/0sjocWkeFEh8ZQmmVk6wS/791HzfJEJGeIvI3EZkjIjtFJEdE1orIuyIyQURC/RGoUkop1ZTtL67kiemeNTEePb+XX9fEqIvI/2aYWDEuo84kQ0QGiMgsYAUwApgPPAM8ALyNZ+GtfwD7ROReTTaUUkoFKmMMD36WSanDyZk9k/lNnzZWh3TQwXEZFiQZx1qM6zPgaeBiY0xhXSeJyDDgTuBuPEmHUkopFVDenr+DHzfkEB0WxGMX9Eak6UzEPDDDxIq1Mo6VZHQ1xlQfrwBjzAJggYhY3y6klFJK+dnafcU8Pn0DAE/9vg8pMWEWR3S4A2tlbLRghkmd3SWHJhgiUmcyIiJxR56vlFJKBYKKaie3friCapebywZ3YNypTaeb5IADYzI255ThdLn9eu36bpC2UEQuM8ZsPvSgiJwBvAX4d0F2FRAqq11k5jnZvWAHWSVVZBU7yCmtwuU2hAfbCQuxExFsJ6lVKJ0To+iSFEXn1pG0snBOulIqsDzyxVq25ZbTLTmKv/z6FKvDqVV0WDBtY8PZW1TJzoIKurSO8tu165tkbAZWiMhdxpjXRSQYeBK4BXjKZ9GpgLM1t4yfNuQwZ1Mui7YXUO10w9K1DSqjXVw4g9PiGdQpnsGd4umcGNmk+keVUi3DFyv38t9lewgNsvHSZQMID7FbHVKduqdEs7eoko1ZpU0vyTDGXCYiVwEviciv8bRcxAJnGGPm+jA+FSCW7SzgpR+3MHtj7sFjItCplY2hPdqS3CqMlFZhJMeEEWyzUVnj8jyqnewtqmJbbhlbc8vZllvGnsJK9hTu5dMVewFIbhXK2B5JjO2RzIj0BCJC6ptbK6VU7dbuK+aBT9cA8PCvTzk47qGp6p4SzY8bctiYVerXLp16/7U1xrwjIr2AewAnMFYTDHWy5m/J48UfN7NwWwEA4cF2zu2dwqjurRmZnsiapQsYPbpPvctzuQ0bs0pZvD2fJTsKWbS9gOwSBx8u3s2Hi3cTEmRjRJcEzj21DWf1TCauCcxjV0o1L/uLK/njW0sor3ZxQb9ULh/SweqQjqu7RTNM6pVkiEgi8B9gJDARGAXMEJH7jTEv+C481VKVVNXwyBdr+czb2hAdGsTEEWlMGtHppBawsduEU1JbcUpqKyaO6IQxhnX7S/hxfQ6zNuSwancRP23M5aeNudhtwvAuCZzbuw1n90omIUqXelFKHVuZw8kf31pKdomDwWnxPHVRn2bRHWvVbqz1bclYA2wC+hljdgHviMg3wOsiMs4Yc7bPIlQtzqJt+dz10Sr2FlUSFmzj5tHp/GFEmk8GbIoIvVJj6JUaw61ndCWntIpZ63L4NnM/87fm88vmPH7ZnMfDX2QytHM8552aqgmHahFKq2rIK6smv8zh+bfcQWmVk3KHkzKHkwqHixqXGwO4jcEYCLILESF2IkOCCA+xExMeTOvoUJKiw0iKDiUlJoyw4KY77sDXnC43t3ywnPX7S+iUGMm/rhxIaFDzuB+dW0ditwk78supqnH57ftY3yTjFeBxY8zBuS/GmE9FZBEwxSeReYnITXi6aNoAa4E7jDG/HOP8UcBkoBeejdyeNsa81pAyvauXPgNcBoQDPwA3GWP2NGLVAk6Ny82zMzbxr5+3Ygz0aRfDc+P7+XUQUlJ0GBOGdGDCkA4Ullczc1020zP3M3dzHvO25DNvSz4Pfb6GIZ0SOKd3Cr/qlUybmHC/xdfUVVQ7yS+rpqiihtKqGkodTkqrnFRUO3HUuKl2uXE43UdNkxOBELudkCAboUE2QoNtRIUG/e8RFkRcRAhxESFNevBcU1JcUcOuggp2FVSws6Cc3QUV7C2qIqu4kv1FVZQ6nI1+TRFIjQmnc+tIOidGkp4cTZ+2MfRoE91s3mxPlDGGR75cy+yNucRHhjBl4qBm1d0aGmSnc2Ikm3PK2JJTRu+2MX65bn0Hfj5Wx/G9wK8aNaJDiMh44AXgJmCu999vReQUb4vKked3Aqbj6dq5Ak/3zqsikmuM+aQBZT4P/BZPkpGPJ2n5WkQGGmNcvqpvS1ZV4+Lm95fzw4YcbAK3jE3ntjO6Emw/kY2AG0dcZAiXDGrPJYPaU1RRzYx12Uxfs595W/JYsC2fBdvyeeTLtfRtH8sZPZIY2yOJXqmtmkXTaEMZY8gvr2ZPYSW7CyrYV1RJVkkV2SVV7C+uIrfUQX5ZNZU1vv/xDw2yER8ZQmJUKK2jQ0mMCqGioJpdoTsODgBOiQkjMSoUu63lfS8OVVJVw468chbsc7Jy1iZ25lewPa+cHfnlFFXUHPO1YcE2kqLDSIgKISEylITIEGIigokMCSIy1E5kaBDBdhs28SQPglDjclNR7aKi2jOourCihtxSz9TxnFIHWcVV7C2qZG9RJb9szjt4rWC70D0lmr7tYhnaOYGhnRNoHd1yWgNdbsMDn65h2lLPuK7XrxxIWmKk1WE1WLeUaDbnlLEhq9T6JENEOhljttenEPH81W1njNndaJF53AW8ZYx5w/v1rSJyDnAjcH8t598A7DPG3Or9er2IDMGz5Pkn9SlTRGKAq4FJxpiZ3vpdCewEzgS+b9Qa1mFvUSVv/LyN7P3VjB7tjyv6TmlVDde8vZRF2wuIjQjm9SszGNwp3uqwDhMbEcIlGe25JKM9xZU1/Lghm+8ys5izKZdVu4tYtbuIyTM3kRQdyujurRneJZFhXRJIbtW0VvY7FpfbsK+okh355ezIr2CX99+d+eXsLqisVwIREmQjMTKE2IgQosOCiA4LJjosiIgQO6FBdkKDbYTYbQTb5bBkzO02VLvcVDs9LR1VNS7Kq12UVdVQ5m0NKayoprCiBofTzf5iT3JzqK+3HT6VOcgmJLcKo02MJ+loExNGm5hwz7+x4aS0CqN1dNNORNxuQ26Zg90FFewurGBn/oFHOTvzK8gvP3SNw8OWKSIixE6H+Ag6xEfQPj6C9nHhtIuLoE1sGKkx4cRGBDd6QlzjcrO7wJPobMstZ31WCWv2FLMlt4zMvSVk7i3h/UWez2rpSVEM75LAmO5JDOuS0Gy7WRxOF3dMXcm3mVmEBdv4vysGkpHWtP5+1VeP5Gi+Yb9fx2VIXfvLi0gW8A3wpnfp8NrOiQMuBW4DXjHGvNxogXmWKa8ALjPG/PeQ468AvY0xo2p5zc/AGmPMzYccuxj4AIgA5HhlishYPN0jScaY3EPOWQt8bIx5pJbrXgdcB5CcnDxw6tSpJ1d5ILfCzT0/V9IqxPDiWP91JzS2kmrD5KVV7ChxExsq3DMojLZR9W+9KCsrIyrKuvo7nIa1+S5W5bpYneui0HH470tKhNA93k7nGBudYmykRtkIauQ3tYbcA7cxFFQZcioM2eVusivcZFcYsivc5FQYnMdY7C8iCFpH2EgMF+LDhPgwG3Fhnv/HhgrRIUKYHZ+25BhjqHZ5fm5Kqw3F1YZihyGn1EGFO5hCh6GwylDocFNajzWGBYgJFeK8dYgJEWJCPY/oECEqWIgKEaKDISJYCLY1Tv1cbkO5E0q99SipNpQcjN1QWOWmoMqQX2lw1v4nGIAQGyRFCAmhbtq2CiE5UkiJsJEc4alDU2lVq3Iadpa42VzoYkOBm01FLqoPyVmDbdAzwU6/1nb6J9mJC2t4C6YVfwscTsNLKxxk5rsID4I7B4bRLc6aZKkx6r8s28lLKxz0TrRzd0bjfUAaM2bMMmNMRm3PHau7pCeeHVe/ERE3sAzPGIcqIA44xXvOYjxjGhr7E34iYAeyjziejadFoTYpwKxazg/ylif1KDMFcAF5tZyTUttFjTGvA68DZGRkmNGN0PTgdhseWfg9JdUu+g4a3qz6/g7ILqliwhsL2VHipmNCBO9dPYT28RENKmP27Nk0xv08GQdGNRtjWL+/lJ8357JwWz5LtheQVeEiq8LJHO9ondAgGz3btCI9KYpOiZF0aR1Jp8QoUlqF0So86ITeFA7cA2MMpQ7PmIj93n73/cWV7CuuOtjNsaewkhpX3e9araND6ZQYSVpCBB0TIklLiKRjgueTcEx4010ptbafg6oa18HunP3FlewvriKruIp93vuSXVJFXlk1RQ5DkeMY7+SHCLbLwTEiEcFBhATZPA+7jaAjWmeMMTicntaZaqebKqeL0ionZVXOBnUrJUSG0M7bCtEhIeLg96ZjQgTJ0WHYbNIkfg8aotrpZvWeIn7enMdPG3JYs7eY1d5E/Z11MKBDLOf0TuGcXm3okFC/vwn+vge7Cyq45cMVZOZXkBgVwtt/HEyvVP90MdSmMep/SmkV8e3206ddLAM7xjVOYMdRZ5Lh3Xn1HhH5C3AenvENHfEMhMzDs93798aYTH8EGmhsNqFbcjQrdxexMbuUoZ0TrA6pQcodTv741hK25pbTIyWad/44mKRm1LVQG5H/TY+9YVQXnC43mftKWLK9gNV7i1mzp4gd+RWs3F3Eyt1FR70+JMhGa+84g+iwIMKD7USE2L0DHQXjHeHvNoaKGhcVDiflDhdZ+ZXULPiB/LJqquux70Dr6FA6xkeQ5k0mPP9GkpYYSVRoy1mILCzY7n0zrrtvvNrpJrfMM5Ygt9RBbpnD829pFQXl1RSW11BYUU1BeTUlVTXUuAyFFTUUHme8w/HYxLOUc0JkCAlRIcRHhpAQFeoZT9Lqf907bePCW+TicCFBNjLS4slIi+eus7qRU1rF7I25zFqXzZxNuSzfVcTyXUU8Pn0Dp7aNYdypbTjv1PonHL5kjOHT5Xt55Mu1lDmctI0N592rB9PZjwPUfSUpOoxJIzr59ZrHGpNxOjDfGFMJfOx9+FMenhaF5COOJwNZdbwmq47znd7ypB5lZuFp7UgEco84p85ZLb7Q/UCSkdW8kgyny82tH65g7b4S0hIi+ODaoSe19kVTFWS30a99LP3axx48VlxRw9r9xQf7rA8M0sspcVDmcB4cNNdwnvEJESF24iNDSGnlGXeQ6n2zau/tl2/XQt+0TlRIkI22seG0ja3fDCGH00VZlSe5K692elopvGNJao6aMSOE2G0HZ8yEBduIDgsmKtQzRqWpdGU0BUnRYQfHPJU7nMzZlMu3mVn8uD6bNXuLWbO3mKe+20Dvtq0Yd2obzu3dhk4WDKwsqqjmwc8y+WbNfgDO6ZXCE787tVm2JDcVx/pr9BOeKZ45IrINGGSMyfdPWJ5dXUVkGXAW8N9DnjqL/w3iPNIC4MIjjp0FLDXG1ADUo8xlQI332Afe17TD0zU0/0TrcyKs3J73RBlj+NvX6/hxQw5xEcFMmTS4RSYYdYmJCGZ4l0SGd0k86rmKaic5JQ7yyhyUe0fvV1S7DjatC4KI51NwWLCdqNAgIkOD2LBmJWeePoyEyFCd3uljoUF2QqPsJDT/D61NVmRoEONObcO4U9tQVeNizqZcpq/Zz6x12QcHjz793UZ6pERzbu82nHVKMj3b+HbJ7nKHk/cX7eT1n7eTV+YgMsTOo+f34qKB7TRZPEnHSjIKgU5ADpDGMbaF96HJwLsishiYh2f2SCrwGoCIvANgjLnKe/5rwC0i8jzwL2AEnhVKL6tvmcaYYhH5N/C0iOTwvymsqzl6vIdPHUwy/LwM7Mn4z7wdvLNgJyF2G69flWHJp5GmKiIkiLTEoAZPfavaZaddnPXNyEo1trBgO2f3SuHsXikHE47vMrOYtS6bDVmlbMgq5blZm2gTE0b3Vk6cSdkM6RxPdCMt3FdcUcNb83cwZf72g1OCMzrGMfmSfk2i66YlOFaS8QkwR0T2AwZYKiK1jmYyxnT2RXDGmGkikgA8hKdVJRMYZ4zZ6T2lwxHnbxeRccBzeKak7gNuO7BGRj3LBLgDTxfLNP63GNdV/l4j4+AysFmlGGOafEb944ZsHvtmHQDPXNKXQc10mpdSyv8OTTgcThfzt+TzXWYWP27M8Q7uhdnvLMUmcEpqKwanJTC4Uxw927SiXVxEvaYqG2PYmlvOnE25np2et+Xj8E65GtgxjlvGpDO6e+sm/7e2OTlWknED8CXQFc8n+SmA3z9SG2NeBV6t47nRtRybAww40TK9zzuAW70PyyRGhRIdAqUOJ/uKq+rdr2yFrOIq/vTRKoyBP53VjfP7plodklKqmQoNsjOmRxJjeiThdhvW7ivh398tYmd1JGv2FB/sVvnPPM9STiF2G2mJEXROjCI2Ipgw76DqkCAbRRU1ZBVXkVXiWUgst9Rx2LVGpidy85h0hnaO1+TCB441u8TgWScDEekLPGuMaT7t9i1Euygb6wvcbMoqbbJJhsttuH3qCgorahjVrTU3j0m3OiSlVAthswmntovhgvQQRo8eQUW1kxW7ili8vYDluwrZklPG/uIqNmWXsSm77LjlxUeGcHrXRO9Oz61b1MqkTVF9lxWf5OtAVO3aRXuSjA1ZpYzpkWR1OLV65actLNpeQOvoUJ69pC+2JrzColKqeYsICWJEeiIj0v83uLrM4WR7bjnb88sPrlNS5X3EhAd7lqOP8Uwfbhsbrn+j/EjnujVxB1bH9Pf2vPW1ZEcBz8/ahAg8d0k/EnX3UqWUn0WFBnFquxhObWfdYlmqdtbtTqXqpV2051u0oQnOMCmqqOb2D1fgNnDDqC6M7Hr0tE2llFKBS5OMJu5AS8bWnLKjts+22l++WMu+4ir6d4jlrrO6WR2OUkqpJkaTjCYuPEhoFxdOtcvNjvxyq8M5aNa6bL5ctY/wYDvPj+9n6ZbtSimlmiZ9Z2gGuicfWJTr+COn/aGkqoaHPvdsWXP32d2PuXeEUkqpwFXvJENEtovIjCOOzRKRrY0fljrU/1b+LLE4Eo8npm8gq8TTTTJxeJrV4SillGqiGjK7ZA6eFTQPtQTY3XjhqNo0pT1M5m/N48PFuwix23j6933qtcqeUkqpwFSvJENEIoCrj1xW2xhzv0+iUodpKnuYVFQ7ue+TNQDcOjadrsm+3bRIKaVU83bc7hIRsQPFQHffh6Nq0zkxiiCbsLOggspqv26fcpjnZ21mV0EFPVKiuWF0F8viUEop1TwcN8nwtl7sBAJnv+4mJiTIRufWkRgDm3Osac3YkFXCv+duRwSevqiPziZRSil1XPV9p/g78KSI6GpLFumWbF2XidtteOizTFxuw5VDO9KnXazfY1BKKdX81Hfg591AJ2CviOwBDluwwRjTp7EDU4frnhzN1+y3JMn4ePkelu4sJDEqlD/9SnvNlFJK1U99k4yPfRqFOi6rZpgUllfzxPT1ADx0Xk9iwoP9en2llFLNV313Yf2rrwNRx3YgyVi/vwRjDCL+mTr69PcbKayoYWjneH7bL9Uv11RKKdUyNGgXVhEZC5wCGGCtMWa2L4JSR+sQH0FKqzCySqpYu6+E3m19v9vg8l2FTF2yiyCb8NgFvf2W2CillGoZ6jXwU0TaishiYCZwL3Af8IOILBIR/XjrByLCmB5JAPy4Icfn13O63Dz0WSbGwLWndyY9SdfEUEop1TD1nV3yIuAC0o0x7Y0x7YGu3mMv+io4dbixfkwy3l24k3X7S2gbG86tY9N9fj2llFItT327S84CRhtjth84YIzZJiK3AT/4JDJ1lBHpCYQE2Vi1p4i8MgeJUaE+uU52SRXPztgEwKPn9yIipEG9akoppRTQsF1YTT2PKR+JCAliWOcEjIHZG3N9dp3HvllPmcPJmT2TOOuUZJ9dRymlVMtW3yTjB+AlEWl/4ICIdACeR1sy/OqMnge6TLJ9Uv7czXl8tWofYcE2HvlNL59cQymlVGCob5JxGxAJbBORnSKyE9jqPXabr4JTRxvT3ZNk/LIpj2qnu1HLdjhd/OWLTABuHduV9vERjVq+UkqpwFLfzvZ8YDAwGujhPbbeGDPLF0GpurWPj6BbchSbsstYuqOA4emNt9L7v+ZsY1teOV1aR3LtaZ0brVyllFKB6bhJxiG7sPY1xszEM41VWWhsj2Q2ZZfxw4acRksyNmSV8NKPmwH4+wW9CQnSDdCUUkqdHN2FtRk6MJX1p0aaylrjcnP3f1dR4zJcNrgDw7voPnhKKaVOXpPdhVVEQkXkJRHJE5FyEflSRNrV43U3ich2EakSkWUiclpDyxWRF0RkqbeMHY1ctZM2oEMsMeHBbMsrZ3te+fFfcByvzd5K5l7PmhgPntezESJUSiml6p9k3A2MxLML61YRWX3ow0exPQ/8HrgMOA1oBXzt7b6plYiMB14AHgf6A/OBb70zYRpSrg14G3inkerSqILsNkZ1aw2c/MJc6/eX8KK3m+Tpi/oQFaprYiillGocTXIXVhGJAa4GJnnHgSAiV+LptjkT+L6Ol94FvGWMecP79a0icg5wI3B/fcs1xtzqfe5u4FeNX8OTd0bPJL5ctY8fN2Rz9chOJ1TGod0klw/pwIhGHESqlFJK1WfgZzCeqaqvGGN2+j4kAAYCwcCMAweMMbtFZD0wnFqSDBEJ8b7umSOemuF9zQmV21SN6tYau01YsDWftfuK6ZXa8A3TXvpxC2v3ldAuLpz7x2k3iVJKqcZ13CTDGFMjIjcCr/ohngNS8OyLknfE8Wzvc7VJBOzec458zZknUe5xich1wHUAycnJzJ49+0SLOkpZWVmd5Y1tb2fmTie3vj2fh4aGYWvALqmL9zt5dZUDgAnphqUL5jZGuI3uWPUPFHoP9B4Eev1B70FzrX99u0tmAGOB/5zMxUTkMeDB45w25mSuYQVjzOvA6wAZGRlm9OjRjVb27Nmzqau8jGFOznx2DtuKq9gb1okrh6XVq8xlOwt5c9ZCAO4/twfXj+rSSNE2vmPVP1DoPdB7EOj1B70HzbX+9U0yfgAeF5E+wDLgsCkNxphP61nO88B7xzlnFzAUT6tEInDoJh3JwC91vC4PTyvFkZttJANZ3v9nnUC5TVZUaBCPnn8KN7y3nKe/28jZvVJIahV2zNfszC/n2neWUu10M2FIB647XRfdUkop5Rv1TTJe9v5b2xLiBs8b93EZY/I4uqviKCKyDKjBs/vrB95j7YCeeGaM1FZ2tfd1ZwH/PeSps4BPvP9vcLlN3dm9UjijRxI/bMjhb1+v4+UJA+o8t6iimklvLaGgvJpR3Vrzt/N7IQ3oYlFKKaUaol5TWI0xtmM86pVgNIQxphj4N/C0iJwpIv2Bd4HVwMGlzEVkg4jccshLJwMTReQaEekpIi8AqcBrDSw3XUT6eV8bIiL9vI8mtyCZiPDo+b0IC7bx9er9zNlU++6sC7flM/5fC9mWW06PlGhentCfILuu6qmUUsp3mvKiCHcATmAaEI6ny+Yq7wqkB3TH0/UBgDFmmogkAA8BbYBMYNwRs2LqU+6bwKhDvl7h/bcTsOMk69Xo2sdHcMeZ3Xjy2w3cPnUFv+2bytm9UxicFk9WSRVPTN/AN2v2A9AhPoL/TBxEdFiwxVErpZRq6Y6ZZIjIfDxv0kXer58A/mmMKfB+nQgsN8Z0qLuUE2OMcQC3eh91nXNUW78x5lWOMROmnuWObkisTcHVIzvx04YcFm0v4O0FO3l7wU7iIoKpqHbhcLoJC7Zx46h0rju9M+Ehjd74pJRSSh3leC0ZQzl8z5KbgTeAAu/XdqCtD+JSDRRst/HhtUNZvbeY7zKz+H5t1sElx8/vm8p95/YgNTbc4iiVUkoFkoZ2l+gowSbMZhP6tY+lX/tY7j2nO1tzyxARurSOsjo0pZRSAagpj8lQJ0FESE+KtjoMpZRSAex40wuM93HkMaWUUkqpYzpeS4YA74mIw/t1GPCGiFR4vw71WWRKKaWUatbEmLobJkRkSn0KMcZMarSImjkRycWzq2tjSaQeC5i1YIFef9B7AHoPAr3+oPegKde/ozGmdW1PHDPJUNYTkaXGmAyr47BKoNcf9B6A3oNArz/oPWiu9dclH5VSSinlE5pkKKWUUsonNMlo+l63OgCLBXr9Qe8B6D0I9PqD3oNmWX8dk6GUUkopn9CWDKWUUkr5hCYZSimllPIJTTKUUkop5ROaZDRRInKTiGwXkSoRWSYip1kdkz+JyOki8qWI7BURIyITrY7Jn0TkfhFZIiIlIpIrIl+JSG+r4/IXEblZRFZ7618iIgtE5Dyr47KK9+fBiMjLVsfiLyLyqLfOhz6yrI7L30SkjYi87f07UCUi60RklNVx1ZcmGU2QiIwHXgAeB/oD84FvRaSDpYH5VxSQCdwOVFocixVGA68Cw4GxgBOYJSLxVgblR3uAe4EBQAbwI/C5iPSxNCoLiMhQ4DpgtdWxWGAj0OaQx6nWhuNfIhILzMOzxcd5QE/gViDHwrAaRGeXNEEisghYbYy59pBjm4GPjTH3WxeZNUSkDLjFGPOW1bFYRUSigGLgAmPMV1bHYwURKQDuN8b8y+pY/EVEYoDlwDXAI0CmMeYWa6PyDxF5FLjIGBMwLXhHEpHHgVHGmBFWx3KitCWjiRGREGAgMOOIp2bg+VSrAlM0nt/XQqsD8TcRsYvIpXhat+ZbHY+fvY7nw8VPVgdikc4iss/bdTxVRDpbHZCfXQAsEpFpIpIjIitF5BYREasDqy9NMpqeRMAOZB9xPBtI8X84qol4AVgJLLA4Dr8RkVO9rVgO4DXgQmPMGovD8hsRuRZIBx6yOhaLLAImAucA1+L5+zdfRBKsDMrPOgM3AduAs/H8HXgSuNnKoBrieFu9K6UsJiKTgZHASGOMy+p4/Ggj0A+IAS4C3haR0caYTEuj8gMR6Y5nTNZIY0yN1fFYwRjz7aFfi8hCPG+2fwAmWxKU/9mApYd0k68Qka54koxmMQhYWzKanjzABSQfcTwZCLiR1YFORJ4DLgPGGmO2WR2PPxljqo0xW4wxy7x/ZFcCd1oclr8Mw9OquVZEnCLiBEYBN3m/DrU2PP8zxpQBa4GuVsfiR/uBdUccWw80m0kAmmQ0McaYamAZcNYRT51F4PVHBzQReYH/JRgbrI6nCbABgfLm+jmemRT9DnksBaZ6/19tSVQWEpEwoAeeN95AMQ/ofsSxbsBOC2I5Idpd0jRNBt4VkcV4fshuAFLx9EsHBO9sinTvlzagg4j0AwqMMbssC8xPROQV4Eo8A78KReTAeJwy7ye6Fk1EngS+AXbjGfQ6Ac+03oBYK8MYUwQUHXpMRMrx/Py3+O4iABF5BvgK2AUkAQ8DkcDbVsblZ8/hGYfyIDANz5IGtwEPWBpVA+gU1iZKRG4C/oxnbngmcKcx5mdro/IfERkN1Dai/m1jzES/BmMBEanrF/OvxphH/RmLFUTkLWAMnsF+xXjWiPinMeZ7K+OykojMJrCmsE4FTsfTbZQLLAQeNsYc2X3QonkXoXscT4vGLjxjMV4yzeTNW5MMpZRSSvmEjslQSimllE9okqGUUkopn9AkQymllFI+oUmGUkoppXxCkwyllFJK+YQmGUoppZTyCU0ylFJKKeUTmmQopZRSyic0yVBKKaWUT2iSoZRqlkSktYjsF5FHDjnWR0SqRORiK2NTSnnosuJKqWZLRM7Gs4nWKDxbwS8FFhtjJlkZl1LKQ5MMpVSzJiLPA+cDc4DTgH6BsFOtUs2BJhlKqWZNREKBVUBXYLgxZpHFISmlvHRMhlKquUsD2gMG6GxtKEqpQ2lLhlKq2RKRYGAhsAlYBDwC9DXG7LI0MKUUoEmGUqoZE5EngQlAH6AY+BYIA8YaY9xWxqaU0u4SpVQzJSKjgD8BVxljioznE9NE4BTgXitjU0p5aEuGUkoppXxCWzKUUkop5ROaZCillFLKJzTJUEoppZRPaJKhlFJKKZ/QJEMppZRSPqFJhlJKKaV8QpMMpZRSSvmEJhlKKaWU8on/BxmpzITPLS6kAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot f(x) og p(x) and the interpolation points\n",
    "plt.subplot(2, 1, 1)\n",
    "plt.plot(x, f(x), x, p, xdata, ydata, \"o\")\n",
    "plt.legend([\"f(x)\", \"p(x)\"])\n",
    "plt.grid(True)\n",
    "\n",
    "# Plot the interpolation error\n",
    "plt.subplot(2, 1, 2)\n",
    "plt.plot(x, (f(x) - p))\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"Error: f(x)-p(x)\")\n",
    "plt.grid(True)\n",
    "print(f\"Max error is {max(abs(p - f(x))):.2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 1: Interpolation of $\\tfrac{1}{1+x^2}$\n",
    "\n",
    "Repeat the previous experiment with Runge's function\n",
    "\n",
    "$$\n",
    "f(x) = \\frac{1}{1+x^2}, \\qquad x\\in [-5,5].\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Insert your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "## Error Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Taylor polynomials once more.**\n",
    "Before we turn to the analysis of the interpolation error\n",
    "$e(x) = f(x) - p_n(x)$, we quickly recall (once more)\n",
    "Taylor polynomials and their error representation.\n",
    "For $f \\in C^{n+1}[a,b]$ and $x_0 \\in (a,b)$,\n",
    "we defined the $n$-th order Taylor polynomial $T^n_{x_0}f(x)$\n",
    "of $f$ around $x_0$ by\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "T^n_{x_0}f(x) &:= \\sum_{k=0}^{n} \\frac{f^{(k)}(x_0)}{k!}(x-x_0)^k\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note that the Taylor polynomial is in fact a polynomial of order $n$\n",
    "which not only interpolates $f$ in $x_0$, but also\n",
    "its first, second etc. and $n$-th derivative $f', f'', \\ldots f^{(n)}$ in $x_0$!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "So the Taylor polynomial the unique polynomial of order $n$ which\n",
    "interpolates the *first $n$ derivatives*\n",
    "of $f$ in a *single point $x_0$*. In contrast,\n",
    "the interpolation polynomial $p_n$ is the unique polynomial of order $n$\n",
    "which *interpolates only the $0$-order* (that is, $f$\n",
    "itself), but in *$n$ distinctive points* $x_0, x_1,\\ldots x_n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For the Taylor polynomial $T^n_{x_0}f(x)$ we have the error\n",
    "representation\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "f(x) - T^n_{x_0}f(x) = R_{n+1}(x_0) \\qquad\n",
    "\\text{where }\n",
    "R_{n+1}(x_0) = \\frac{f^{(n+1)}(\\xi)}{(n+1)!} (x-x_0)^{n+1},\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "with $\\xi$ between $x$ and $x_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Of course, we usually don't know the exact location of $\\xi$\n",
    "and thus not the exact error,\n",
    "but we can at least estimate\n",
    "it and bound it from above:\n",
    "\n",
    "$$\n",
    "|f(x) - T^n_{x_0}f(x)| \\leqslant\n",
    " \\frac{M}{(n+1)!} h^{n+1}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "M=\\max_{x\\in[a,b]}|f^{(n+1)}(x)| \\qquad \\text{and} \\qquad h = |x-x_0|.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The next theorem gives us an  expression for the interpolation\n",
    "error $e(x)=f(x)-p_n(x)$ which is similar to what we have just\n",
    "seen for the error between the Taylor polynomial and the original function\n",
    "$f$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theorem 1: Interpolation error\n",
    "\n",
    "Given $f \\in C^{(n+1)}[a,b]$. Let $p_{n} \\in \\mathbb{P}_n$ interpolate $f$ in\n",
    "$n+1$ distinct nodes $x_i \\in [a,b]$. For each $x\\in [a,b]$ there is at least\n",
    "one $\\xi(x) \\in (a,b)$ such that\n",
    "\n",
    "$$\n",
    "f(x) - p_n(x) = \\frac{f^{(n+1)}(\\xi(x))}{(n+1)!}\\prod_{i=0}^n(x-x_i).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Proof.**\n",
    "We start fromt the Newton polynomial $\\omega_{n+1} =: \\omega(x)$\n",
    "\n",
    "$$\n",
    "\\omega(x) = \\prod_{i=0}^{n}(x-x_i) = x^{n+1} + \\dotsm.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, the error in the nodes, $e(x_i)=0$. \n",
    "Choose an *arbitrary* $x\\in [a,b]$, $x\\in [a,b]$, where $x\\not=x_i$,\n",
    "$i=0,1,\\dotsc,n$. For this fixed $x$, define a function in $t$ as:\n",
    "\n",
    "$$\n",
    "\\varphi(t) = e(t)\\omega(x) - e(x)\\omega(t).\n",
    "$$\n",
    "\n",
    "where $e(t) = f(t)-p_n(t)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Notice that $\\varphi(t)$ is as differentiable with respect to $t$ as $f(t)$. The\n",
    "function $\\varphi(t)$ has $n+2$ distinct zeros (the nodes and the fixed x). As a\n",
    "consequence of [Rolle's theorem](https://en.wikipedia.org/wiki/Rolle's_theorem),  the derivative\n",
    "$\\varphi'(t)$ has at least $n+1$ distinct zeros, one between each of the zeros\n",
    "of $\\varphi(t)$. So $\\varphi''(t)$ has $n$ distinct\n",
    "zeros, etc. By repeating this argument, we can see that $\\varphi^{n+1}(t)$\n",
    "has at least one zero in $[a,b]$, let us call this $\\xi(x)$, as it does depend on the fixed $x$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Since\n",
    "$\\omega^{(n+1)}(t)=(n+1)!$ and $e^{(n+1)}(t)=f^{(n+1)}(t)$ we obtain\n",
    "\n",
    "$$\n",
    "\\varphi^{(n+1)}(\\xi)= 0 = f^{(n+1)}(\\xi)\\omega(x) - e(x)(n+1)!\n",
    "$$\n",
    "\n",
    "which concludes the proof."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Observation.** The interpolation error consists of three elements: The derivative of the\n",
    "function $f$, the number of interpolation points $n+1$ and the distribution of\n",
    "the nodes $x_i$. We cannot do much with the first of these, but we can choose\n",
    "the two others. Let us first look at the most obvious choice of nodes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Equidistributed nodes\n",
    "\n",
    "The nodes are *equidistributed* over the interval $[a,b]$ if $x_i=a+ih$, $h=(b-a)/n$. In this case it can\n",
    "be proved that:\n",
    "\n",
    "$$\n",
    "|\\omega(x)| \\leq \\frac{h^{n+1}}{4}n!\n",
    "$$\n",
    "\n",
    "such that\n",
    "\n",
    "$$\n",
    "|e(x)| \\leq \\frac{h^{n+1}}{4(n+1)}M, \\qquad M=\\max_{x\\in[a,b]}|f^{(n+1)}(x)|.\n",
    "$$\n",
    "\n",
    "for all $x\\in [a,b]$. \n",
    "\n",
    "Let us now see how good this error bound is by an example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 2: Interpolation error for $\\sin(x)$ revisited\n",
    "\n",
    "Let again $f(x)=\\sin(x)$ and $p_n(x)$ the polynomial interpolating $f(x)$ in\n",
    "$n+1$ equidistributed points on $[0,1]$.\n",
    "An upper bound for the error for different values of $n$\n",
    "can be found easily. Clearly,\n",
    "$\\max_{x\\in[0,2\\pi]}|f^{(n+1)}(x)|=M=1$ for all $n$, so\n",
    "\n",
    "$$\n",
    "|e_n(x)| = |f(x)-p_n(x)| \\leq\n",
    "\\frac{1}{4(n+1)}\\left(\\frac{2\\pi}{n}\\right)^{n+1}, \\qquad x\\in[a,b].\n",
    "$$\n",
    "\n",
    "Use the code in the first Example of this lecture to verify the result\n",
    "for $n = 2, 4, 8, 16$. How close is the bound to the real error?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Insert your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "## Optimal choice of interpolation points\n",
    "So how can the error be reduced? For a given $n$ there is only one choice: to\n",
    "distribute the nodes in order to make\n",
    "$|\\omega(x)|= \\prod_{j=0}^{n}|x-x_i|$ as small as possible. We will first do this\n",
    "on a standard interval $[-1,1]$, and then transfer the results to some arbitrary\n",
    "interval $[a,b]$.\n",
    "\n",
    "Let us start taking a look at $\\omega(x)$ for equidistributed nodes on the\n",
    "interval $[-1,1]$, for\n",
    "different values of $n$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "newparams = {\"figure.figsize\": (6, 3)}\n",
    "plt.rcParams.update(newparams)\n",
    "\n",
    "\n",
    "def omega(xdata, x):\n",
    "    # compute omega(x) for the nodes in xdata\n",
    "    n1 = len(xdata)\n",
    "    omega_value = np.ones(len(x))\n",
    "    for j in range(n1):\n",
    "        omega_value = omega_value * (x - xdata[j])  # (x-x_0)(x-x_1)...(x-x_n)\n",
    "    return omega_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Plot omega(x)\n",
    "n = 8  # Number of interpolation points is n+1\n",
    "a, b = -1, 1  # The interval\n",
    "x = np.linspace(a, b, 501)\n",
    "xdata = np.linspace(a, b, n)\n",
    "plt.plot(x, omega(xdata, x))\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"omega(x)\")\n",
    "print(f\"n = {n:2d}, max|omega(x)| = {max(abs(omega(xdata, x))):.2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Run the code for different values of $n$. Notice the following: \n",
    "* $\\max_{x\\in[-1,1]} |\\omega(x)|$ becomes smaller with increasing $n$. \n",
    "\n",
    "* $|\\omega(x)|$ has its maximum values near the boundaries of $[-1, 1]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "A a consequence of the latter, it seems reasonable to move the nodes towards the boundaries. \n",
    "It can be proved that the optimal choice of nodes are the *Chebyshev-nodes*, given by\n",
    "\n",
    "$$\n",
    "\\tilde{x}_i = \\cos \\left( \\frac{(2i+1)\\pi}{2(n+1)} \\right), \\qquad i=0,\\dotsc,n\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let $\\omega_{Cheb}(x) = \\prod_{j=1}^n(x-\\tilde{x}_i)$. It is then possible to prove that\n",
    "\n",
    "$$\n",
    "\\frac{1}{2^{n}} = \\max_{x\\in [-1, 1]} |\\omega_{Cheb}(x)| \\leq \\max_{x \\in [-1, 1]} |q(x)|\n",
    "$$\n",
    "\n",
    "for all polynomials $q\\in \\mathbb{P}_n$ such that $q(x)=x^n + c_{n-1}x^{n-1}+\\dotsm+c_1x + c_0$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The distribution of nodes can be transferred to an interval $[a,b]$ by the linear transformation\n",
    "\n",
    "$$\n",
    "x = \\frac{b-a}{2}\\tilde{x} + \\frac{b+a}{2}\n",
    "$$\n",
    "\n",
    "where $x\\in[a,b]$ and $\\tilde{x} \\in [-1,1]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "By doing so we get\n",
    "\n",
    "$$\n",
    "\\omega(x) = \\prod_{j=0}^n (x-x_i) =\n",
    "  \\left(\\frac{b-a}{2}\\right)^{n+1} \\prod_{j=0}^n (\\tilde{x}-\\tilde{x}_i)\n",
    "  = \\left(\\frac{b-a}{2}\\right)^{n+1} \\omega_{Cheb}(\\tilde{x}).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "From the theorem on interpolation errors we can conclude:\n",
    "\n",
    "**Theorem (interpolation error for Chebyshev interpolation).**\n",
    "\n",
    "Given $f \\in C^{(n+1)}[a,b]$, and let $M_{n+1} = \\max_{x\\in [a,b]}|f^{(n+1)}(x)|$.  Let $p_{n} \\in \\mathbb{P}_n$ interpolate $f$ i $n+1$ Chebyshev-nodes $x_i \\in [a,b]$. Then\n",
    "\n",
    "$$\n",
    "\\max_{x\\in[a,b]}|f(x) - p_n(x)| \\leq \\frac{(b-a)^{n+1}}{2^{2n+1}(n+1)!} M_{n+1}.\n",
    "$$\n",
    "\n",
    "The Chebyshev nodes over an interval $[a,b]$ are evaluated in the following function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def chebyshev_nodes(a, b, n):\n",
    "    # n Chebyshev nodes in the interval [a, b]\n",
    "    i = np.array(range(n))  # i = [0,1,2,3, ....n-1]\n",
    "    x = np.cos((2 * i + 1) * pi / (2 * (n)))  # nodes over the interval [-1,1]\n",
    "    return 0.5 * (b - a) * x + 0.5 * (b + a)  # nodes over the interval [a,b]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 3: Chebyshev interpolation\n",
    "\n",
    "\n",
    "**a)**\n",
    "Plot $\\omega_{Cheb}(x)$ for $3, 5, 9, 17$ interpolation points on the interval $[-1,1]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Insert your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**b)**\n",
    "Repeat Example 3 using Chebyshev interpolation on the functions below. Compare with the results you got from equidistributed nodes.\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "  f(x) &= \\sin(x), && x\\in[0,2\\pi] \\\\ \n",
    "  f(x) &= \\frac{1}{1+x^2}, && x\\in[-5,5]. \n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Insert your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**For information**: \n",
    "[Chebfun](http://www.chebfun.org/) is software package which makes it possible to manipulate functions and to solve equations with accuracy close to machine accuracy. The algorithms are based on polynomial interpolation in Chebyshev nodes.\n",
    "<!-- --- end exercise --- -->"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}