10_Trigonometric_Interpolation_And_The_FFT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ★ Trigonometric Interpolation And The FFT ★"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import modules\n",
    "import math\n",
    "import cmath\n",
    "import numpy as np\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 10.1 The Fourier Transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lemma : Primitive roots pf unity\n",
    "\n",
    "Let $\\omega$ be a primitive bth root of unity and k be an integer. Then\n",
    "\n",
    "### $\\sum_{j=0}^{n-1}\\omega^{jk} = \n",
    "\\left\\{\\begin{matrix}\n",
    "n & \\text{if k/n is an integer } \\\\ \n",
    "0 & \\text{otherwise} \n",
    "\\end{matrix}\\right.$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definition\n",
    "\n",
    "The Discrete Fourier Transform (DFT) of $x = [x_0,\\cdots,x_{n-1}]^T$ is the n-dimensional vector $y = [y_0,\\cdots,y_{n-1}]$, where $\\omega = e^{-i2\\pi/n}$ and\n",
    "\n",
    "### $y_k = \\frac{1}{\\sqrt{n}}\\sum_{j=0}^{n-1}x_j\\omega^{jk}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definition\n",
    "\n",
    "The magnitude of a complex vector $v$ is the real number $||v|| = \\sqrt{\\bar{v}^Tv}$. A square complex matrix F is unitary if $\\bar{F}^TF=I$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "\n",
    "Find the DFT of the vector $x = [1, 0, -1, 0]^T$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.+0.j 1.+0.j 0.-0.j 1.+0.j]\n"
     ]
    }
   ],
   "source": [
    "np.set_printoptions(precision=3)\n",
    "x = np.array([1, 0, -1, 0]).T\n",
    "w = complex(math.cos(math.pi * 2 / x.size), -math.sin(math.pi * 2 / x.size))\n",
    "F = np.empty((x.size, x.size), dtype=complex) # Fourier matrix\n",
    "for i in range(x.size):\n",
    "    for j in range(x.size):\n",
    "        F[i, j] = pow(w,(i * j))\n",
    "y = (1 / math.sqrt(x.size)) * np.matmul(F, x)\n",
    "print(np.round(y)) # [0, 1, 0, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.+0.j 1.+0.j 0.+0.j 1.+0.j]\n"
     ]
    }
   ],
   "source": [
    "# Or use numpy fft (fast fourier transform)\n",
    "print(np.fft.fft(x) / np.sqrt(x.size))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lemma \n",
    "\n",
    "Let $\\{y_k\\}$ be the DFT of $\\{x_j\\}$, where the $x_j$ are real numbers. Then \n",
    "\n",
    "1. $y_0$ is real\n",
    "2. $y_{n-k} = \\bar{y_k} \\text{ for k} = 1,\\cdots,n - 1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 10.2 Trigonometric Interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lemma\n",
    "\n",
    "Let $t = j/n$, where j and n are integers. Let k be an integer. Then $\\cos2(n-k)\\pi t = \\cos2k\\pi t$ and $\\sin2(n - k)\\pi t = - sin{2k\\pi t}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Corollary\n",
    "\n",
    "For an even integer n, let $t_j = c + j(d - c)/n$ for $j = 0,\\cdots,n - 1$, and let $x = (x_0,\\cdots,x_{n-1})$ denote a vector of n real numbers. Define $\\overrightarrow{a} + \\overrightarrow{b}i = F_nx$, where $F_n$ is the Discrete Fourier Transform. Then the function \n",
    "\n",
    "### $P_n(t) = \\frac{a_0}{\\sqrt{n}} + \\frac{2}{\\sqrt{n}}\\sum_{k=1}^{n/2-1}(a_k\\cos{\\frac{2k\\pi (t-c)}{d-c}}-b_k\\sin{\\frac{2k\\pi (t-c)}{d-c}})+\\frac{a_{n/2}}{\\sqrt{n}}cos{\\frac{n\\pi (t-c)}{d-c}}$ \n",
    "\n",
    "satisfies $P_n(t_j) = x_j$ for $j = 0,\\cdots,n-1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def even_trigonometric_interpolant(t, a, b, c, d):\n",
    "    return a[0] / np.sqrt(a.size) + 2 / np.sqrt(a.size) * \\\n",
    "        np.sum([a[k] * np.cos(2 * k * np.pi * (t - c) / (d - c)) - \\\n",
    "        b[k] * np.sin(2 * k * np.pi * (t - c) / (d - c)) for k in range(1, int(a.size / 2))]) + \\\n",
    "        a[int(a.size / 2)] / np.sqrt(a.size) * np.cos(a.size * np.pi * (t - c) / (d - c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example \n",
    "\n",
    "Find the trigonometric interpolant for $x = [1, 0, -1, 0]^T$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xd4VGX6//H3nU4ooYViKpDQQUqoFkBZirpgXRVUrNjYXUVdddldXZUVy9eOuljWhgWx4SoiVRRFCdIJkFASCJCEXkJLcv/+yOAvYQNJmMmcKffruubKzDnPmfkcEnLnOc85zxFVxRhjjDkuxOkAxhhjfIsVBmOMMeVYYTDGGFOOFQZjjDHlWGEwxhhTjhUGY4wx5VhhMMYYU44VBmOMMeVYYTDGGFNOmNMBTkfjxo01OTnZ6RjGGONXFi9evENVYytr55eFITk5mfT0dKdjGGOMXxGR7Kq0s0NJxhhjyrHCYIwxphwrDMYYY8qxwmCMMaYcKwzGGGPK8UhhEJE3RSRfRFaeZL2IyAsikiUiy0WkW5l1o0Qk0/UY5Yk8FZk8eTLJycmEhISQnJzM5MmTa+qjjDHGr3mqx/AWMOQU64cCqa7HaOAVABFpCDwE9AJ6Ag+JSAMPZfrN5MmTGT16NNnZ2agq2dnZjB492oqDMcZUwCOFQVXnA7tO0WQ48I6WWgjUF5HmwGBgpqruUtXdwExOXWBOy7hx4ygsLCy3rLCwkL/+9a+e/ihjjPGog0eKWL5lD58s3sKE6Ws4cKSoxj/TWxe4xQGby7ze4lp2suX/Q0RGU9rbIDExsVofnpOTc5Llmxk+cQGpTeqUPprWIS25IfWiwqv1/sYY4y5VZWXuPjK27SMzfz+Z+QfIzDtA7p5Dv7UJDxWGdzmDds3r1WgWv7nyWVUnAZMA0tLStDrbJiYmkp39vxf8xcQ2o05kKPPXFTB18RYA6kaFcUPfZG48uwX1oyM8kNwYY05OVZm7Np8X52SxJGcPAJFhIbSKrUNacgOubpJASpO6pDatQ1LDaMJCa/6cIW8VhlwgoczreNeyXKD/CcvnefrDx48fz+jRo8sdToqOjmbis08xcmRvAPYWHmP1tn28/eMmXpiTxRs/bOTaPsncfE4LGteJ9HQkY0yQKylRvl2dx0tzM1mZu4+4+rV4dHgH+rVuQlyDWoSGiGPZvFUYpgFjRORDSgea96rqNhGZAfyrzIDzIOBBT3/4yJEjgdKxhpycHBITExk/fvxvywFiosPp06oRfVo1Yu32/bw0N4t/z1/PWz9uZETPJG7t15Km9aI8Hc0YE2SKS5SvV2zjpTlZrM3bT1KjaJ68rDOXdIsj3Au9gaoQ1Wodlan4TUQ+oPQv/8ZAHqVnGoUDqOqrIiLAS5QOLBcCN6hqumvbG4Hjo8DjVfU/lX1eWlqaemMSvfUFB5g4N4svlm4lNES4Mi2Bu3/Xmoa17RCTMab6/rt8K8/OXMf6goO0iq3NmPNS+H3nM7xyeAhARBaralql7TxRGLzNW4XhuJydhbzyXRZTF2+hSd0o/n1tdzrGxXjt840x/u1IUTEPT1vNB7/k0LZZXcacl8LQjs29frjICkMNWL5lD7e+u5hdB48y4bJOXNI13usZjDH+JW/fYW5/bzG/5uzh9v6tuHdQG8fGD6paGHzjgJaf6Bxfny//eDZnJtTn7o+W8ciXqykqLnE6ljHGRy3O3sVFL/7Amu37mTiiG/cPaevooHJVWWGopsZ1Ipl8cy+u75vMmws2cs0bP7PzwBGnYxljfMz7P+dw1aSFREeE8tkdZ3Fh5+ZOR6oyKwynITw0hIeHdeD/rjiTX3P2MOylBazM3et0LGOMDzhSVMyDny7nr5+toG+rxky782zaNKvrdKxqscLghsu6x/PJbX1RVS575Uc+/XWL05GMMQ7K23eYqyct5INfNnPngFa8eX0PYqL9byYFKwxu6hQfw5d/PJuuifUZO2UZ7/60yelIxhgH5O07zOWv/sia7ft5ZWQ37hvsH+MJFbHC4AGN6kTy7k29GNiuKf+Ytoppy7Y6HckY40V7C49x3Ru/sOvAUT64pTdDO/nPeEJFrDB4SHhoCC+N6EqP5IbcM2Up89cVOB3JGOMFh44Wc9Pbi9i44yCTrkvjzIT6TkdymxUGD4oKD+X1UWmkNKnLre8uZknObqcjGWNq0LHiEu6YvJjFObt57qounJXS2OlIHmGFwcPqRYXz9o09aFIvkhveWkRm3n6nIxljakBJifKXqcuZu7aA8Rd34gI/P3xUlhWGGtCkbhTv3tiL8NAQrnvzl3LzqRtj/J+q8thXGXy2JJd7B7VmRK/q3SPG11lhqCGJjaJ558aeHDhSxLV2EZwxAeXleet5c8FGbjgrmTsHpDgdx+OsMNSgds3r8caoHuTuPsQNby3yyi35jDE16/2fc3hqxlou6RrH3y9sT+nk0YHFCkMN69miIS+P7Maqrfu47d3FNreSMX7s21Xb+dvnKxjQJpYnL+9MiJ9ep1AZKwxecH67pjx+SSd+yNrB87MznY5jjDkNOTsLuWfKMjrFxfDyyO4+c1OdmuCRPRORISKyVkSyROSBCtY/KyJLXY91IrKnzLriMuumeSKPL/pDjwSu6B7PS3Oz+CFzh9NxjDHVcLSohD9+8CsIvDSiG7UiQp2OVKPcLgwiEgpMBIYC7YGrRaR92TaqereqdlHVLsCLwKdlVh86vk5Vh7mbx5f9c3gHWsXW4a6PllKw3wajjfEXT81Yw7Ite3nyss4kNIx2Ok6N80SPoSeQpaobVPUo8CEw/BTtrwY+8MDn+p3oiDAmjujG/sPHGDtlKSUl/neTJGOCzZw1ebz2/Uau7Z3k91NdVJUnCkMcsLnM6y2uZf9DRJKAFsCcMoujRCRdRBaKyMUeyOPT2jSry8PDOvB95g5e+W6903GMMaewbe8h7pmyjHbN6zHuwnZOx/Eab4+eXAVMVdXiMsuSXLeaGwE8JyKtKtpQREa7Ckh6QYF/z0N0VY8ELurcnGdmriN90y6n4xhjKlBUXMKfP1jKkaISJo7oSlR4YI8rlOWJwpALJJR5He9aVpGrOOEwkqrmur5uAOYBXSvaUFUnqWqaqqbFxsa6m9lRIsLjl3Yirn4t/vTBEnYfPOp0JGPMCV6Ynckvm3bx2MUdaRlbx+k4XuWJwrAISBWRFiISQekv//85u0hE2gINgJ/KLGsgIpGu542Bs4DVHsjk8+pGhfPSiK4UHDjCfVOXo2rjDcb4igVZO3hxbhaXd4/n0m7xTsfxOrcLg6oWAWOAGUAGMEVVV4nIIyJS9iyjq4APtfxvwHZAuogsA+YCE1Q1KAoDQOf4+jwwtB2zMvL4z4JNTscxxgAF+49w10dLadm4No8M7+B0HEeEeeJNVPVr4OsTlv3jhNcPV7Ddj0AnT2TwVzeelcxP63fy+PQM0pIb0Dne/+dyN8ZflZQoY6csZd+hY7xzY0+iIzzyK9LvBO6le35CRHjq8s40rhPJnz9cyuFjxZVvZIypEW/9uInvM3fwj9+3p13zek7HcYwVBh/QoHYET17emY07DjJxbpbTcYwJSrl7DvH0t2vp3yaWET0Daxrt6rLC4CPOSY3l0q5xvDJvPWu32819jPEmVeXvn69EFR4d3jEgZ0ytDisMPuRvF7WnblQYD3y63K6KNsaLvlqxjTlr8rlnUOugmPKiMlYYfEjD2hH8/aL2LMnZw+Sfs52OY0xQ2Ft4jIenraZTXAzX9012Oo5PsMLgYy7pGsfZKY154pu1bN972Ok4xgS8Cd9ksLvwKI9f2omwAJ5KuzrsX8HHiAjjL+nIseISHpq20uk4xgS0nzfs5INfNnPT2S3oGBfjdByfYYXBByU1qs1dA1szY1Ue36zc7nQcYwLSkaJiHvxsBfENanHXwFSn4/gUKww+6uZzWtCueT0emraS/YePOR3HmIDz8tz1bCg4yGMXdwzaC9lOxgqDjwoPDWHCpZ3I33+EJ79Z63QcYwJKZt5+Xp6XxfAuZ9C/TROn4/gcKww+7MyE+lzfN5n3fs5mcbZNz22MJ5SUKA9+uoLoiDD+flH7yjcIQlYYfNw9g9rQvF4UD366gqNFJU7HMcbvfbAoh/Ts3Yy7sB2N60Q6HccnWWHwcXUiw3j04o6syzvAa99vcDqOMX4tf/9hJkxfQ5+Wjbiie/BNp11VVhj8wPntmjK4Q1Mmzs0ib59d22DM6Xrm23UcOlrMY5fYtBenYoXBT/z1gnYUFStPz7CBaGNOx+qt+/gofTPX9UmmVZDdka26rDD4iaRGtbn+rGSm/rqFlbl7nY5jjF9RVR77ajUxtcL58/l2zUJlPFIYRGSIiKwVkSwReaCC9deLSIGILHU9bi6zbpSIZLoeozyRJ1CNOS+FBtERPPrf1XYrUGOqYVZGPj+u38ld56cSEx3udByf53ZhEJFQYCIwFGgPXC0iFZ0D9pGqdnE9Xndt2xB4COgF9AQeEpEG7mYKVPWiwrn7d635eeMuZqzKczqOMX7haFEJ//o6g1axtRnZO8npOH7BEz2GnkCWqm5Q1aPAh8DwKm47GJipqrtUdTcwExjigUwB6+oeCbRuWofHp2dwpMju9mZMZd5dmM3GHQcZd2E7wm2SvCrxxL9SHLC5zOstrmUnukxElovIVBFJqOa2xiUsNIRxF7Yne2ch7/xoU3Mbcyq7Dx7l+VnrOCe1MQPsCucq81b5/BJIVtXOlPYK3q7uG4jIaBFJF5H0goICjwf0J/1ax9K/TSwvzMlk54EjTscxxmc9PzuTA0eK+NuF7e301GrwRGHIBRLKvI53LfuNqu5U1eO/wV4Huld12zLvMUlV01Q1LTY21gOx/dvfLmxH4dFinpuV6XQUY3xSVv4B3l2YzVU9E2nTrK7TcfyKJwrDIiBVRFqISARwFTCtbAMRaV7m5TAgw/V8BjBIRBq4Bp0HuZaZSqQ0qcvIXom8/0sOmXl2j2hjTvSvrzOIDg9l7O9aOx3F77hdGFS1CBhD6S/0DGCKqq4SkUdEZJir2Z9EZJWILAP+BFzv2nYX8CilxWUR8IhrmamCuwa2JjoilMe+yqi8sTFB5PvMAuasyefO81JsPqTTIP54PnxaWpqmp6c7HcMnvP79Bh77KoO3buhh0wcbAxQVl3DhCz9QeKyIWWP7ERkW6nQknyEii1U1rbJ2du6Wn7uuTzLJjaJ57KsMiopt9lVjPkrfzNq8/Tw4tJ0VhdNkhcHPRYSF8OAF7cjKP8AHizZXvoExAezAkSKe+XYdPZMbMrRjM6fj+C0rDAFgUPum9EhuwIuzMzl01C56M8HrzR82svPgUR68oK2dnuoGKwwBQES4b3Bb8vcf4Z2fNjkdxxhH7D54lNfmb+B37ZvSNdFm1nGHFYYA0bNFQ/q3ieWV79az7/Axp+MY43WvfreeA0eLuHdQG6ej+D0rDAHk3kFt2FN4jNfn253eTHDJ23eYt37cxMVd4uxiNg+wwhBAOsbFcGGn5rz+w0Z22FQZJoi8OCeT4hLlroF2rwVPsMIQYMYOas3hY8W8PHe901GM8YrsnQf58JfNXNUzgaRGtZ2OExCsMASYVrF1uLx7PO8tzCZ3zyGn4xhT456blUlYqPDH86y34ClWGALQn1y3Lnxxtk2wZwLb2u37+XxpLqP6JtO0XpTTcQKGFYYAFN8gmhG9Evl48RY2FBxwOo4xNebpb9dSJyKM285t5XSUgGKFIUDdOSCFyLAQnpm5zukoxtSIJTm7mbk6j1vObUmD2hFOxwkoVhgCVGzdSG48qwX/Xb6NVVv3Oh3HGI97+tu1NKodwY1nt3A6SsCxwhDAbjm3JTG1wnl6xlqnoxjjUQuydrAgayd3DEihTmSY03ECjhWGABZTK5zb+rVi7toCFm2y21yYwKCqPDljLWfERDGyV6LTcQKSFYYAN6pvErF1I3nqm7X44703jDnRzNV5LNu8hz8PTCUq3KbVrgkeKQwiMkRE1opIlog8UMH6sSKyWkSWi8hsEUkqs65YRJa6HtNO3Na4JzoijD+el8Ivm3YxP3OH03GMcUtJifLMzHW0aFyby7rFOx0nYLldGEQkFJgIDAXaA1eLSPsTmi0B0lS1MzAVeLLMukOq2sX1GIbxuCt7JHBGTBTPz1pnvQbj175ZtZ012/dz18BUwkLtgEdN8cS/bE8gS1U3qOpR4ENgeNkGqjpXVQtdLxcCVuq9KDIslDsGpPBrzh7rNRi/VVKiPD8rk5axtbmo8xlOxwlonigMcUDZW4dtcS07mZuA6WVeR4lIuogsFJGLT7aRiIx2tUsvKChwL3EQuiItnjNionjOeg3GT01fuZ21efv58/mphIbYTXhqklf7YiJyDZAGPFVmcZLr5tQjgOdEpMJLGFV1kqqmqWpabGysF9IGlsiwUO48L4Ul1mswfqikRHl+9jpaWW/BKzxRGHKBhDKv413LyhGRgcA4YJiq/jYntKrmur5uAOYBXT2QyVTgiu4J1mswfmn6yu2syzvAn6y34BWeKAyLgFQRaSEiEcBVQLmzi0SkK/BvSotCfpnlDUQk0vW8MXAWsNoDmUwFIsJCrNdg/I71FrzP7cKgqkXAGGAGkAFMUdVVIvKIiBw/y+gpoA7w8QmnpbYD0kVkGTAXmKCqVhhq0BXdE4irX8t6DcZvWG/B+zxyLbmqfg18fcKyf5R5PvAk2/0IdPJEBlM1EWEh3Dkghb9+toLv1hXQv00TpyMZc1LHewspTepYb8GL7ETgIHR593hXryHTeg3Gp329cpv1FhxghSEIHe81LN28h+/W2am/xjcdv24hpUkdLuzU3Ok4QcUKQ5CyXoPxdV+v3EZmvvUWnGCFIUhZr8H4MustOMsKQxA73mt41noNxsd8taK0t2BXOTvDCkMQiwgLYcx5KSzbvId51mswPqK4RHlhdiapTepwgfUWHGGFIchd1s3GGoxv+XqFjS04zQpDkCvba7CxBuO0khLlxTmlYwvWW3COFQbDZd1KZ159Ybb1GoyzZqwqvcr5j+elWG/BQVYYDBFhIdzuul/Dj+t3Oh3HBKnSq5wzadnY5kRymhUGA8Af0uJpVi+K52dnOh3FBKlZGXms2b6fOwdYb8FpVhgMUHq/hlv7teSXjbtYuMF6Dca7VJUX52SR2DCa4V2st+A0KwzmN1f3TKRxnUhenGO9BuNd89YWsCJ3L2MGpNi9nH2AfQfMb6LCQ7mtX0sWZO0kfdMup+OYIKFaOrYQV78Wl3Q71V2BjbdYYTDljOiVSKPaEbwwJ8vpKCZIfJ+5g6Wb93DHgFaEW2/BJ3jkuyAiQ0RkrYhkicgDFayPFJGPXOt/FpHkMusedC1fKyKDPZHHnL7oiDBuPqcl89cVsHTzHqfjmACnWnqVc/OYKC7vHu90HOPidmEQkVBgIjAUaA9cLSLtT2h2E7BbVVOAZ4EnXNu2p/RWoB2AIcDLrvczDrq2TxIlmd9zdtd2hISEkJyczOTJk52OZQLI5MmTSU5OJjQ0lM8fuJjOR1YQGWb/9X2FJ3oMPYEsVd2gqkeBD4HhJ7QZDrztej4VOF9ExLX8Q1U9oqobgSzX+xkHfTH1I7Z/9TwHd25HVcnOzmb06NFWHIxHTJ48mdGjR5OdnY2qUryvgPeeGmc/Xz7EE4UhDthc5vUW17IK27juEb0XaFTFbY2XjRs3jmNHDpdbVlhYyLhx4xxKZALJuHHjKCwsLLfMfr58i9+M9IjIaBFJF5H0ggKb06cm5eTkVGu5MdVhP1++zxOFIRdIKPM63rWswjYiEgbEADuruC0AqjpJVdNUNS02NtYDsc3JJCYmVmu5MdVhP1++zxOFYRGQKiItRCSC0sHkaSe0mQaMcj2/HJijpbO1TQOucp211AJIBX7xQCbjhvHjxxMdHV1uWVStWowfP96hRCaQjB8/ntCIqHLLoqOj7efLh7hdGFxjBmOAGUAGMEVVV4nIIyIyzNXsDaCRiGQBY4EHXNuuAqYAq4FvgDtVtdjdTMY9I0eOZNKkSSQlJSEihMU04dwbxjFy5Eino5kA0KnfhdQfdCcNmpyBiJCUlMSkSZPs58uHiD9Os5yWlqbp6elOxwgaj0/P4LX5G5g5th+tYus4Hcf4uZvfXsSiTbv54f4B1I0KdzpOUBGRxaqaVlk7vxl8Ns655ZyWRIaFMtGuhjZuWpm7l1kZ+dx0dgsrCj7MCoOpVOM6kVzTO5HPl+ayccdBp+MYP/bC7EzqRoUxqm+y01HMKVhhMFUy+tzSeWwmzrVegzk9q7fu49vVedx4VgtiallvwZdZYTBVEls3kpG9kvhsSS7ZO63XYKrvxTmZ1I0M48azWjgdxVTCCoOpslv7tSQ0RHh57nqnoxg/s3b7fqav3M71ZyUTE229BV9nhcFUWdN6UYzomcgnv25h867CyjcwxuWFOZnUjgi13oKfsMJgquXWfi0JEeHledZrMFWTmbefr1dsY1TfZBrUjnA6jqkCKwymWprH1OLKHglMXbyZ3D2HnI5j/MCLc7KoFR7Kzee0dDqKqSIrDKbabuvfCoBX5tkZSubUsvIP8OXyrVzbJ4mG1lvwG1YYTLXF1a/FFWkJTFm0hW17rddgTm7i3CyiwkK5xXoLfsUKgzktt/drRYkqr9pYgzmJjTsO8sXSXK7pnUjjOpFOxzHVYIXBnJaEhtFc3j2eDxZtJm/f4co3MEHnpTlZhIeGMPrcVk5HMdVkhcGctjv6p1Bcorz6nfUaTHnZOw/y+dJcRvZKIrau9Rb8jRUGc9oSG0Vzadc43v85h3zrNZgyJs7NIjREuK2fjS34IysMxi13DkihqET59/wNTkcxPmLzrkI+/TWXET0TaVIvqvINjM+xwmDckty4Nhd3ieO9hdk21mCA0hlUQ0KE2/rZ2IK/cqswiEhDEZkpIpmurw0qaNNFRH4SkVUislxEriyz7i0R2SgiS12PLu7kMc748/mpFJcoL9vMq0Fv446DfLokl2t6JdEsxnoL/srdHsMDwGxVTQVmu16fqBC4TlU7AEOA50Skfpn196lqF9djqZt5jAMSG0VzRVo8H/xiV0MHu+dnrSM8VLitv40t+DN3C8Nw4G3X87eBi09soKrrVDXT9XwrkA/Euvm5xsfcOSAFRe1+DUEsK38/Xyzbyqg+yTSpa70Ff+ZuYWiqqttcz7cDTU/VWER6AhFA2fMbx7sOMT0rInZem5+KbxDNVT0SmbJos828GqSenZVJdHgot9rYgt+rtDCIyCwRWVnBY3jZdqqqgJ7ifZoD7wI3qGqJa/GDQFugB9AQuP8U248WkXQRSS8oKKh8z4zX3TkghZAQ4YXZmU5HMV6WsW0fXy3fxg1ntbA5kQJApYVBVQeqascKHl8Aea5f+Md/8edX9B4iUg/4ChinqgvLvPc2LXUE+A/Q8xQ5JqlqmqqmxcbakShf1Cwmimt6JfHpErs3dLB5duY66kaGcfM5dr+FQODuoaRpwCjX81HAFyc2EJEI4DPgHVWdesK640VFKB2fWOlmHuOw2/u3IjzUeg3BZMWWvXy7Oo+bzmlB/WjrLQQCdwvDBOB3IpIJDHS9RkTSROR1V5s/AOcC11dwWupkEVkBrAAaA4+5mcc4LLZuJKP6JPP50lyy8vc7Hcd4wbOz1hFTK5wbz7beQqAIc2djVd0JnF/B8nTgZtfz94D3TrL9ee58vvFNt/ZrxXsLs3l2ViYTR3RzOo6pQb/m7GbOmnzuG9yGelF2L+dAYVc+G49rWDuCG85qwVfLt5GxbZ/TcUwNenbmOhrWjuD6vslORzEeZIXB1IhbzmlJ3agwnpu1zukopob8snEX32fu4LZ+Lakd6dbBB+NjrDCYGhETHc7NZ7dkxqo8VubudTqOqQHPzFxLbN1Iru2d7HQU42FWGEyNufHsZGJqhfPMTOs1BJofs3awcMMu7ujfiloRoU7HMR5mhcHUmLpR4Yw+tyVz1uTza85up+MYD1FV/m/mOprVi+LqnolOxzE1wAqDqVHX902mcZ0InvxmDaUXxxt/Nzsjn8XZuxlzXgpR4dZbCERWGEyNqh0Zxp/OT2Xhhl3MW2tTmfi7ouISnvhmDS0b1+bKHglOxzE1xAqDqXFX90wkuVE0E6avobjEeg3+7JNft5CZf4D7BrchPNR+fQQq+86aGhceGsJ9g9uyNm8/n/66xek45jQdOlrMMzPX0TWxPkM6NnM6jqlBVhiMV1zQqRlnJtTnmZnrOHys2Ok45jS8uWAjefuO8ODQdpROb2YClRUG4xUiwoND27Jt72He+nGT03FMNe06eJRX561nYLsm9GzR0Ok4poZZYTBe07tlI85r24SJc7PYffCo03FMNbw4J5ODR4u4f0hbp6MYL7DCYLzq/iFtOXikyG4B6kdydhby3sJsruieQGrTuk7HMV5ghcF4VZtmdbmsWzzv/JRttwD1E09/u5bQEOHu37V2OorxEisMxuvGDmqNCDZVhh9YsWUv05Zt5cazWtAsJsrpOMZLrDAYr2seU4sbzmrB50tzWbXVJtjzVarKhG8yaBAdzm39Wzkdx3iRW4VBRBqKyEwRyXR9bXCSdsVl7t42rczyFiLys4hkichHrtuAmiBwe/9WxNQKZ8L0NU5HMScxP3MHC7J2Mua8VLsJT5Bxt8fwADBbVVOB2a7XFTmkql1cj2Fllj8BPKuqKcBu4CY38xg/EVMrnDEDUvg+cwc/ZO5wOo45QUmJMmH6GuIb1OKa3jZRXrBxtzAMB952PX8buLiqG0rpFTLnAVNPZ3vj/67tk0Rc/Vo8Pj2DEpsqw6d8vjSXjG37uG9wGyLDbKK8YONuYWiqqttcz7cDTU/SLkpE0kVkoYgc/+XfCNijqkWu11uAuJN9kIiMdr1HekGBTcYWCCLDQrlvcBtWbd3Hp0tynY5jXAqPFvH0jLV0jKvH7zuf4XQc44BKC4OIzBKRlRU8hpdtp6VzKp/sz74kVU0DRgDPiUi1R7JUdZKqpqlqWmxsbHU3Nz5q2Jln0DWxPhOmr2Hf4WNOxzHAxLlZbN17mId+34GQEJv6IhhVWhhUdaCqdqzg8QWQJyLNAVxf80/yHrmurxuAeUBXYCdKEzLDAAAQAUlEQVRQX0SO3yw2HrA/G4NMSIjw6PCO7Dx4hGft9FXHbdxxkNfmb+TSrnH0SLapL4KVu4eSpgGjXM9HAV+c2EBEGohIpOt5Y+AsYLWrhzEXuPxU25vA1zEuhhE9E3nnp2zWbN/ndJygpao8PG0VEWEhPDDUpr4IZu4WhgnA70QkExjoeo2IpInI66427YB0EVlGaSGYoKqrXevuB8aKSBalYw5vuJnH+Kn7BrehXlQY//hild3pzSEzV+fx3boC7hqYSpN6djFbMAurvMnJqepO4PwKlqcDN7ue/wh0Osn2G4Ce7mQwgaF+dAT3DW7LXz9bwbRlWxne5aTnIZgacPhYMY/8dzWtm9ZhVN9kp+MYh9mVz8ZnXNkjgc7xMYz/KoMDR4oq38B4zCvz1rNl9yH+Oayj3ZnNWGEwviM0RHhkeEfy9x/hhdmZTscJGjk7C3nlu/X8/swz6NOqkdNxjA+wwmB8SpeE+lyZlsCbP2wkM2+/03GCwiP/XUVYiDDugnZORzE+wgqD8Tl/GdKG6IhQHv7SBqJr2pw1eczKyOdP56fa7KnmN1YYjM9pVCeSewe3YUHWTr5esd3pOAHr8LFi/vnlalrF1ubGs1o4Hcf4ECsMxieN7JVE++b1eOyr1Ry0gega8dr8DWTvLOThYR2ICLNfBeb/s58G45NKB6I7sG3vYbsNaA3YsruQifOyGNqxGeek2hQzpjwrDMZnpSU35NJucbz2/Qa7ItqDVJV/fLEKgL9d1N7hNMYXWWEwPm3cBe2IqRXO2I+WcbSoxOk4AeHj9C3MWZPPXwa3Ja5+LafjGB9khcH4tEZ1Ihl/SSdWb9vHS3Ps2gZ3bdldyCP/XU2vFg253q5wNidhhcH4vMEdmnFptzgmzlvPss17nI7jt0pKlL9MXY6q8vQVZ9qU2uakrDAYv/DQ7zvQpG4k93y8jMPHip2O45feXZjNj+t38reL2pPQMNrpOMaHWWEwfiGmVjhPXNaZrPwDPD1jrdNx/M6GggM8Pj2D/m1iuapHgtNxjI+zwmD8xrmtY7mmdyJvLNjIzxt2Oh3HbxSXKPd8vIzIsFCeuKwzpbdbN+bkrDAYv/Lg0HYkNIjm3qnL7MK3Kpo0fwNLcvbwyPAONLX7LJgqsMJg/ErtyDCevuJMtuw+xL++znA6js9bs30fz85cx9COzRh25hlOxzF+wq3CICINRWSmiGS6vjaooM0AEVla5nFYRC52rXtLRDaWWdfFnTwmOPRs0ZCbz27B5J9z+G5dgdNxfNbRohLumbKMerXCeOzijnYIyVSZuz2GB4DZqpoKzHa9LkdV56pqF1XtApwHFALflmly3/H1qrrUzTwmSNwzqA0pTepw/9Tl7D10zOk4PumluVms2rqP8Zd0olGdSKfjGD/ibmEYDrztev42cHEl7S8HpqtqoZufa4JcVHgoz/zhTAoOHOEfX6y06blPsCRnNxPnZnFptzgGd2jmdBzjZ9wtDE1VdZvr+XagaSXtrwI+OGHZeBFZLiLPishJ/6wRkdEiki4i6QUFdvjAQOf4+vz5/FS+WLqV/yzY5HQcn5G/7zC3vbeYZvWieOj3HZyOY/xQpYVBRGaJyMoKHsPLttPSP9lO+mebiDQHOgEzyix+EGgL9AAaAvefbHtVnaSqaaqaFhtrs0GaUmMGpDCofVPGf53BD5k7nI7juMPHirn1vcXsO1TEa9elEVMr3OlIxg9VWhhUdaCqdqzg8QWQ5/qFf/wXf/4p3uoPwGeq+tsBYVXdpqWOAP8Berq3OybYhIQIz1zZhVaxtbnz/V/ZtOOg05Eco6r87fOVLMnZwzN/OJP2Z9RzOpLxU+4eSpoGjHI9HwV8cYq2V3PCYaQyRUUoHZ9Y6WYeE4TqRIbx+nU9EIFb3kln/+HgHIx+c8Empi7ewp/OT2Vop+ZOxzF+zN3CMAH4nYhkAgNdrxGRNBF5/XgjEUkGEoDvTth+soisAFYAjYHH3MxjglRio2heHtGNDTsOcvdHSykpCa7B6O8zCxj/1WoGd2jKXeenOh3H+Dnxx7M50tLSND093ekYxge9/eMmHpq2ijEDUrh3cBun43jFph0HGT5xAc1jovjk9r7UjgxzOpLxUSKyWFXTKmtnP0EmoFzXJ4mMbft4aW4WbZvX5aLOgX217/7Dx7j5nXRCBF67Ls2KgvEImxLDBBQR4ZHhHUlLasC9Hy9jZe5epyPVmJIS5e6PlrJxx0EmjuxmU2kbj7HCYAJORFgIr1zTnYbREYx+J50dB444HalG/N/MtczKyOeh37enb6vGTscxAcQKgwlIsXUjmXRdGrsKj3L9f35h18GjTkfyqDd/2MjEueu5umcC1/ZOcjqOCTBWGEzA6hgXwyvXdCcz7wBX/vsn8vYddjqS21SVF2Zn8sh/VzO0YzP+OcwmxzOeZ4XBBLQBbZrw9o092brnEFe8+hObd/nvNF2qyuPT1/DMzHVc1i2eF6/uSkSY/Rc2nmc/VSbg9W7ZiMm39GbvoWNc8epPZOUfcDpStZWUKOM+X8mk+Ru4rk8ST13embBQ++9raob9ZJmg0CWhPh/d2puiEuXKf//kV2crHSsuYeyUpbz/cw539G/FP4d1ICTEDh+ZmmOFwQSNts3q8fFtfYgMC+Hq1xayOHu305EqdaSomDsm/8rnS7dy3+A2/GVIWxtTMDXOCoMJKi0a1+bj2/vSqHYE177xMwuyfHdG1sKjRdz8djozV+fxz2EduHNAitORTJCwwmCCTlz9Wky5rQ+JDaO54T+L+GhRjs/d6CdnZyEjXy8tXE9fcSaj+iY7HckEESsMJig1qRvFh6N70y2pPvd/soJr3/iFnJ3On7FUXKK8Nn8Dg577jsy8A7w8shuXd493OpYJMjaJnglqJSXK+7/kMGH6GopLlHsGteaGs1oQ6sDgbsa2fTzwyXKWbdnLwHZNePTijjSPqeX1HCZwVXUSPSsMxgDb9h7ib5+tZPaafM5MqM8Tl3WibTPv3OjmSFExL83J4pV564mpFc7DwzpwUefmNshsPM4KgzHVpKp8uXwb/5y2ir2HjnHHgBTuHNCKyLDQGvvMxdm7uP+TFWTlH+DSrnH8/aL2NKgdUWOfZ4JbVQuDW2MMInKFiKwSkRIROemHicgQEVkrIlki8kCZ5S1E5GfX8o9ExP5HGMeICMPOPIOZY/sx7MwzeGF2Jhc8/z1vLdjI9r2em07jWHEJ360rYOyUpVz+6k8cOlrMWzf04Jkru1hRMD7B3cHnlcClwPyTNRCRUGAiMBRoD1wtIu1dq58AnlXVFGA3cJObeYxxW8PaETxzZRfeuqEHYSEhPPzlano/PpvLXvmR17/fwJbd/ztIPXnyZJKTkwkJCSE5OZnJkyeXW3+kqJjZGXncM2UZ3R+dyag3f+HbVXlc3zeZGXefS/82Tby1e8ZUyiOHkkRkHnCvqv7P8R0R6QM8rKqDXa8fdK2aABQAzVS16MR2p2KHkow3ZeUf4JuV2/h6xXZWb9sHwJnxMQzt1JxzUhvz7bRP+OvYMRw+dOi3baJq1eJfz7xEu7MvYPrKbczOyOfAkSLqRYUxsH1TLujYnLNTGxMVXnOHqYw5kS/dwS0O2Fzm9RagF9AI2KOqRWWWx3khjzHVktKkDmPOS2XMealk7zzI9JXbmb5iGxOmr2HCdNjyygMUlykKAIcPHeK++x8g/vamNIgO56LOzRnSsRl9WzW2ie+Mz6u0MIjILKBZBavGqeoXno900hyjgdEAiYmJ3vpYY8pJalSb2/q14rZ+rdiyu5AlOXsY/mTFV0+X7N/B+7f0omdyQ5vwzviVSguDqg508zNygYQyr+Ndy3YC9UUkzNVrOL78ZDkmAZOg9FCSm5mMcVt8g2jiG0STmJhIdnb2/6xPTEy0O6sZv+SNP2MWAamuM5AigKuAaVo6uDEXuNzVbhTgtR6IMZ4yfvx4oqPL3285Ojqa8ePHO5TIGPe4e7rqJSKyBegDfCUiM1zLzxCRrwFcvYExwAwgA5iiqqtcb3E/MFZEsigdc3jDnTzGOGHkyJFMmjSJpKQkRISkpCQmTZrEyJEjnY5mzGmxC9yMMSZIeOUCN2OMMYHHCoMxxphyrDAYY4wpxwqDMcaYcqwwGGOMKccvz0oSkQLgf68oqprGgO/e6Ldm2D4HB9vnwOfu/iapamxljfyyMLhDRNKrcrpWILF9Dg62z4HPW/trh5KMMcaUY4XBGGNMOcFYGCY5HcABts/BwfY58Hllf4NujMEYY8ypBWOPwRhjzCkEbGEQkSEislZEskTkgQrWR4rIR671P4tIsvdTelYV9nmsiKwWkeUiMltEkpzI6UmV7XOZdpeJiIqIX5/BUpX9FZE/uL7Pq0TkfW9n9LQq/FwnishcEVni+tm+wImcniQib4pIvoisPMl6EZEXXP8my0Wkm0cDqGrAPYBQYD3QEogAlgHtT2hzB/Cq6/lVwEdO5/bCPg8Aol3Pbw+GfXa1qwvMBxYCaU7nruHvcSqwBGjget3E6dxe2OdJwO2u5+2BTU7n9sB+nwt0A1aeZP0FwHRAgN7Az578/EDtMfQEslR1g6oeBT4Ehp/QZjjwtuv5VOB8EREvZvS0SvdZVeeqaqHr5UJK75rnz6ryfQZ4FHgCOOzNcDWgKvt7CzBRVXcDqGq+lzN6WlX2WYF6rucxwFYv5qsRqjof2HWKJsOBd7TUQkrvhtncU58fqIUhDthc5vUW17IK22jpzYT2UnqzIH9VlX0u6yZK/+LwZ5Xus6uLnaCqX3kzWA2pyve4NdBaRBaIyEIRGeK1dDWjKvv8MHCN66ZhXwN/9E40R1X3/3u1VHrPZxN4ROQaIA3o53SWmiQiIcAzwPUOR/GmMEoPJ/WntEc4X0Q6qeoeR1PVrKuBt1T1/0SkD/CuiHRU1RKng/mrQO0x5AIJZV7Hu5ZV2EZEwijtgu70SrqaUZV9RkQGAuOAYap6xEvZakpl+1wX6AjME5FNlB6LnebHA9BV+R5vofSe6sdUdSOwjtJC4a+qss83AVMAVPUnIIrSOYUCWZX+v5+uQC0Mi4BUEWkhIhGUDi5PO6HNNGCU6/nlwBx1jer4qUr3WUS6Av+mtCj4+7FnqGSfVXWvqjZW1WRVTaZ0XGWYqvrrfWGr8nP9OaW9BUSkMaWHljZ4M6SHVWWfc4DzAUSkHaWFocCrKb1vGnCd6+yk3sBeVd3mqTcPyENJqlokImOAGZSe1fCmqq4SkUeAdFWdBrxBaZczi9JBnqucS+y+Ku7zU0Ad4GPXOHuOqg5zLLSbqrjPAaOK+zsDGCQiq4Fi4D5V9duecBX3+R7gNRG5m9KB6Ov9/I88ROQDSgt8Y9fYyUNAOICqvkrpWMoFQBZQCNzg0c/3838/Y4wxHhaoh5KMMcacJisMxhhjyrHCYIwxphwrDMYYY8qxwmCMMaYcKwzGGGPKscJgjDGmHCsMxhhjyvl/4g2liWl3eZ0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "c, d = 0, 1\n",
    "x = np.array([1, 0, -1, 0]).T\n",
    "y = np.fft.fft(x) / np.sqrt(x.size)\n",
    "a = np.array(y.real)\n",
    "b = np.array(y.imag)\n",
    "f = lambda t, a, b, c, d : a[0] / np.sqrt(a.size) + 2 / np.sqrt(a.size) * \\\n",
    "np.sum([a[k] * np.cos(2 * k * np.pi * (t - c) / (d - c)) - \\\n",
    "        b[k] * np.sin(2 * k * np.pi * (t - c) / (d - c)) for k in range(1, int(a.size / 2))]) + \\\n",
    "a[int(a.size / 2)] / np.sqrt(a.size) * np.cos(a.size * np.pi * (t - c) / (d - c))\n",
    "\n",
    "t = np.linspace(c, d, 40)\n",
    "xt = np.linspace(0, 1, 4, False)\n",
    "data = [f(i, a, b, c, d) for i in t]\n",
    "\n",
    "plt.plot(t, data)\n",
    "plt.plot(xt, x, 'o', color='k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Find the trigonometric interpolant for the temperature data : $x = [-2.2, -2.8, -6.1, -3.9, 0.0, 1.1, -0.6, -1.1]$ on the interval $[0,1]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VPW9//HXJ/tKWLJA1mFHCHsQcKcqUEURd41L1Zpqte2ttba3+d3+7PXG2lqXLrYaK26N+16qYFEERbawCARZEkhCAkkmAUI2ss33/pHoBQ1kmcmcWT7PxyOPRzKZnHkfkrw5+Z7v+R4xxqCUUsp3BFgdQCmllGtpsSullI/RYldKKR+jxa6UUj5Gi10ppXyMFrtSSvkYLXallPIxWuxKKeVjtNiVUsrHBFnxorGxscZms1nx0kop5bU2btxYbYyJ6+55lhS7zWYjPz/fipdWSimvJSIlPXmeDsUopZSP0WJXSikfo8WulFI+RotdKaV8jBa7Ukr5GC12pZTyMVrsSinlYyyZx66Uvzvc0MI7W8qJDA1iTEI0o+OjiAzVX0flGvqTpJQbHWpo4elP9/LC58U0tLSf8LmkgeGMSYjqKPqEaCYkDmDTx/8kOzub0tJSUlNTycnJITMz06L0yltosSvlBtX1zTy9ai8vri2hqbWdiycO4645owgLDmR3ZR17KuvYXVnP7so6VhfW0NLuoL5gBbX/foK25mMAlJSUkJWVBaDlrk5JjDHOb0RkMbAAqDLGpHf3/IyMDKNLCih/YK9rJndVEf9YW0pzWzuXTE7k7jmjGJ0QfdKvaWt3UFzTSEb6GGrtB7/1+bS0NIqLi/sxtfJUIrLRGJPR3fNcdcT+HPAX4AUXbU8pr2aM4bHle8hdVURLm4PLpiRx13dGMTIuqtuvDQoMYFR8FEerK7r8fGlpqavjKh/jkmI3xqwSEZsrtqWUL3jkw938ZUUhCyYN42dzxzI8NrLX20hNTaWk5NtrPqWkpLgiovJhOt1RKRd7bvU+/rKikGtnpPDn66b2qdQBcnJyiIiIOOExCQpl7IIsGlvaXBFV+Si3FbuIZIlIvojk2+12d72sUm71zy8O8JslO5g7PoH/uSwdEenztjIzM8nNzSUtLQ0RITU1lZt+8SCFA6aw6InP2Vfd4MLkype45OQpQOdQzBI9ear81ad77Nz63AampgzihdtOJyw4sN9e50cvb8YYeOeuM/v8F4HyPj09eapDMUq5wNayI9zx4kZGxkXx9M0Z/VbqAGePjuPdu84kQCDrhXzqm3VYRp3IJcUuIi8Da4CxIlImIre5YrtKeYN91Q3c8uwGBkWG8PytpxMTHtzvr5k2JJInrp/G3uoG7nl1Cw6Ha/7yVr7BJcVujLnOGDPMGBNsjEk2xjzjiu0q5emqjh7jxmfWYYAXbj2dhAFhbnvtM0bFkn3RaXy4o5I/fbzHba+rPJ8OxSjVR7VNrdy0eD2HGlp47pYZjOjBHHVXu+VMG1dMS+bx5XtYVtD1vHflf7TYleoDYwz3vLqFIns9T904nUnJAy3JISLkLEpncnIM97y6hT2VdZbkUJ5Fi12pPnh/WwUf7aziF/PHcfboOEuzhAUH8uSN0wkPCeL2F/KpbWy1NI+ynha7Ur1U29TK/f8sYGJSDLecOdzqOAAMiwnnyRumUX6kiR+/spl2PZnq17TYleqlh5ftpKa+mQcXTSQwoO8XILlahm0wv7k0nZW77Ty8bJfVcZSFdNlepXphY8lh8taVcssZw5mYHGN1nG+5fmYqBQdqeXJlEeMTB3Dp5ESrIykL6BG7Uj3U2u4g++1tDBsQxs/mjrE6zkn9/0smkJE2iP98cyulNY1Wx1EW0GJXqof+/uk+dlbU8ZuF6R59G7uQoAD+dN1UAgKEn762hbZ2h9WRlJtpsSvVA6U1jfzxo93Mm5DAheMTrI7TrcSB4fzPZelsLDnMkyuLrI6j3EyLXaluGGP4f+9uJ1CE+y+dYHWcHrt0ciILJg3j8eV72Fp2xOo4yo202JXqxpKtB1m1286988YyLCbc6jg9JiLkXDaR2KhQ/uPVLTR94+bZyndpsSt1CrWNrfzmnzuYlBzDTbNtVsfptZiIYB65ejJ77Q389oMvrY6j3ESLXalT+N2ynRxq8Lw5671x5qhYbj1zOC+sKeGTXVVWx1FuoMWu1ElsLDnES+tKufXM4aQned6c9d64b/5YxiREcd8bWznc0GJ1HNXPtNiV6oIxht99sIuEAaH89ELPnbPeU2HBgTx+zVQON7bwq7e34ao7pynPpMWuVBfWFNWwvvgQd80Z5dFz1ntjfOIAfjZ3LB9sr+DNTeVWx1H9SItdqW8wxvDY8t0MHRDG1RkpVsdxqdvPHsHpwwdz/3sF7D+kV6X6Ki12pb7h86IaNhQf5odzRvbrvUutEBggPHr1ZAT4ySubadWrUn2Sq+55Ol9EdolIoYj80hXbVMoKxhge7zxav2aGbx2tfyV5UAQPXj6RTaVH+MOHugqkL3K62EUkEHgC+C4wHrhORMY7u12l3CkvLw+bzUZgYCBv/+IyprZsIzTIt47Wj3fJ5ESun5nKUyv3smKnToH0Na44Yj8dKDTG7DXGtACvAAtdsF2l3CIvL4+srCxKSkowxtB+1M4LD2eTl5dndbR+9esF4xk3NJp7XtvCwdomq+MoF3JFsScB+4/7uKzzMaW8QnZ2No2NJ55IbGxsJDs726JE7hEWHMgTmdNobnPw45c36yqQPsRtJ09FJEtE8kUk3263u+tllepWaWlprx73JSPjonhw0UQ2FB/mseW7rY6jXMQVxV4OHH+WKbnzsRMYY3KNMRnGmIy4OGtv/qvU8VJTU3v1uK+5bGoS12Sk8NdPili1Ww+6fIErin0DMFpEhotICHAt8J4LtquUW+Tk5BAYHHrCYxEREeTk5FiUyP3uv3QCY+Kj+emrW6g8eszqOMpJThe7MaYNuBtYBnwJvGaMKXB2u0q5S8qMuQycdzdDhiYhIqSlpZGbm0tmZqbV0dwmPCSQJzKn0tjSzk9e2Uy7Q5cc8GYuGWM3xrxvjBljjBlpjPGfwxzl9b66ynTMGd+lfH8JDoeD4uJivyr1r4yKj+aBy9JZu/cQf/xoj9VxlBP0ylPl11btqWZz6RHu+s4on5633lNXTk/mimnJ/PnjPTq/3YtpsSu/9dVVpokxYVw13TevMu2LBy6bwPhhA7jjHxv5vLDa6jiqD7TYld/6rPD/jtZDgvRX4SsRIUG8eNtMbEMiue35fNbvO2R1JNVL+tOs/Nazq4uJjQrhyunJVkfxOIMjQ/jH92cybGAYtzy7nk2lh62OpHpBi135peLqBlbsquL6mWk6tn4ScdGhvPT9WcRGh3Lz4vVsK6u1OpLqIS125ZdeWFNCoAiZM/3jIqS+GhoTxku3zyImPJgbF69jx4GjVkdSPaDFrvxOQ3Mbr+fv56KJw0gYEGZ1HI+XNDCcl2+fRXhwIDc+s449lXVWR1Ld0GJXfuetTWXUNbfxvTNtVkfxGimDI3jp9lkEBAjX/30de+31VkdSp6DFrvyKw2F47vNiJifHMDVloNVxvMrw2Ehe+v5MHA7D3LseZFhyCgEBAdhsNp9f4tjbaLErv/JZYTVF9gZuPsOGiFgdx+uMTojm2sEl7HvnUSrKyzDGUFJSQlZWlpa7B9FiV37l+c87pjhePGmY1VG81hMP/w+O1uYTHmtsbORXv/qVRYnUN2mxK79RXN3AxzrF0WknX79+PweO6J2YPIEWu/IbOsXRNU62Tn1QTBwXPrqSxZ/t09UhLabFrvyCTnF0nZycHCIiIk54LCIigkd+/xAzhg/mv5fs4OI/fcpfPylkd2UdxmjJu1uQ1QGUcged4ug6Xy1pnJ2dTWlpKampqeTk5JCZmYkxhiVbD/LUqiJ+v3QXv1+6i5TB4Zw/LoHzT4tn5vAhui6PG4gV/5tmZGSY/Px8t7+u8k/GGC54dCVRoUG8c9eZOhvGTSpqj/HRzko++rKK1YXVNLc5iAoN4pwxscweGcvo+ChGxkURGxWi35MeEpGNxpiM7p6nR+zK5301xfHRqydrgbjR0JgwMmemkTkzjaaWdlYXVn9d9O9vq/j6edFhQYyM6yj5kfGRjIyLYkLiAJIHRZxi6+pUtNiVz3tutU5xtFp4SCAXjE/ggvEJGGM4WHuMIns9RVX1FNkbKLLX81mhnTc3lX39NYkxYcwYPpgZto630fFRBATof8w94VSxi8hVwP3AacDpxhgdX1EepaSmY4rjj74zWqc4eggRIXFgOIkDwzl7dNwJn6s71speewNb9h9hffEh1hTV8O6WAwAMjAgmI20QM2yDuXRKIsNiwq2I7xWcPWLfDlwOPOWCLEq5nE5x9C7RYcFMThnI5JSB3HyGDWMMpYcaWb/vEBuKD7Gh+DDLv6zikQ93c/WMZO44d6QO2XTBqWI3xnwJ6Lil8kiNLW28plMcvZqIkDYkkrQhkVyV0XH7wv2HGnlyZRGvbtjPK+v3c+X0ZH543ihSh2jBf8Vt845EJEtE8kUk3263u+tllR97f1sFdcfauGFWmtVRlAulDI4gZ9FEVv58DpkzU3lrczlzHvmEe1//gsefXIzNZvP7xcm6ne4oIsuBoV18KtsY827ncz4B7u3pGLtOd1TucG3uGipqj7Hi3vP0r0ofVnn0GE+t3MvfFj9H5ZI/Ydr+bx2biIgIcnNzv5577+16Ot2x2yN2Y8wFxpj0Lt7edU1UpVyvtKaRtXsPceX0ZC11H5cwIIxfXzKeoI2vnlDq0LE4WXZ2tkXJrKOXgCmf9OamMkTg8ml6o2p/UV62v8vHT7ZomS9zqthFZJGIlAGzgX+JyDLXxFKq7xwOwxsbyzhrVCyJA3VKnL842eJk0UOGcqy13c1prOVUsRtj3jbGJBtjQo0xCcaYea4KplRfrd1bQ/mRJq6crkfr/qSrxcmCQ8MInnU9Vz+1hnI/WlJYh2KUz3ljYxnRYUHMm9DVOX/lqzIzM8nNzSUtLa1jmmRaGs8+83deeuhe9tkbuOTPn7G6sNrqmG6hi4Apn1J3rJUZOcu5fFoyDy6aaHUc5SH22uv5wYsbKbLXc9/8cfzgnBFeeVLdZbNilPIm7287yLFWhw7DqBOMiIvinbvO5Lvpw3jog5089MFOn14nXhcBUz7l9fwyRsZFMjVloNVRlIeJDA3iz9dNZXBkCE+t2ktLu4NfLxjvlUfu3dFiVz5jr72e/JLD/GL+OJ/8ZVXOCwgQ/nvhBIIChWdXF9PuMNx/yQSfWzVSi135jDc3lREgcPm0JKujKA8mIvx6wXhCAgN4atVeWtsNOZel+1S5a7Ern9DuMLy5sZxzx8Tpgl+qWyLCL787jqBA4YkVRbS1O3joikkE+ki5a7Ern7C6sJqKo8f4rwXjrY6ivISIcO/csQQFBPDHj/bQ7jA8fNVknyh3LXblE17fWEZMeDAXjI+3OoryIiLCTy8cQ3Cg8IcPd9PqMDx29WSCAr17wqAWu/J6tY2tLCuo4NoZKXqXJNUnd39nNEGBATz0wU7aHQ4ev2YqIUHeW+5a7Mrr/XPrAVraHFw1PcXqKMqL3XHuSIIDA3hgyQ5a2jbyl+unERbsnQcK3vtfklKdXt9YxtiEaNKTBlgdRXm5284azgOXpbP8yypufyGfphbvXDxMi115tT2VdXyx/whXZei668o1bpyVxsNXTmJ1YTXfe3Y99c1tVkfqNS125dXe2FhGYICwcIrOXVeuc1VGCo9fO5X8ksPc9Mw6aptarY7UK1rsymu1Owxvby5nztg44qJDrY6jfMylkxN54vppbCuvJfPvaznc0GJ1pB7TYldea92+GqrqmvVoXfWb+elDyb0xg92V9Vz39Frsdc3df5EH0GJXXmvJ1oOEBwdy/mk6d131nznj4nn2ezMoqWnkms4bpHs6Z2+N97CI7BSRrSLytojoknrKLdraHSzdXsH5p8UTEaKzdlX/OnNULC/cdjpVR5u55C+f8dkez75hh7NH7P8G0o0xk4DdwH86H0mp7n1eVMOhhhYWTEq0OoryEzNsg3njztnEhAdz4+J1/G7pTlrbHVbH6pKz9zz90Bjz1VygtYDe3UC5xZKtB4gKDeK8sXFWR1F+ZNzQAfzz7rO4dkYKf/ukiKufWsP+Q41Wx/oWV46x3wp84MLtKdWllraOYZgLxyd47ZWBynuFhwTy28sn8Zfrp1JYVc9Ff/yUJVsPWB3rBN0Wu4gsF5HtXbwtPO452UAbkHeK7WSJSL6I5NvtdtekV37ps0I7R4+1sWDSMKujKD+2YFIi7//4bEYlRHH3S5v55ZtbPeZK1W7POhljLjjV50Xke8AC4HxzipsIGmNygVzouJl172Iq9X+WfHGQAWFBnD1ah2GUtVIGR/DaD2bz2L9387eVReSXHOY/vzuOc8bEEWzhCpHOzoqZD9wHXGqM6deBpry8PGw2GwEBAdhsNvLyTvrHgfJhx1rb+XBHJfMmDPXq1feU7wgODOC++eN48daZ1B1r5bbn85n92494YMkOCg7UWpLJ2XlifwFCgX93rtOx1hhzh9OpviEvL4+srCwaGzv+7ygpKSErKwuAzMxMV7+c8mArd9upb25jwWSdDaM8y1mjY/n0vu+wcredNzeW8cKaYp75bB/jhkZzxbRkFk5NJD7aPXf3klOMnvSbjIwMk5+f3+Pn22w2SkpKvvX4oPhEfvvaKiJCg4gIDiQyNJA2h+FIYyu1Ta0cbWr9+v2v3oZEhWAbEsnw2EiGx0UyIjaSpIHhXr+wvr/40cubWV1YzbpfnW/pn7pKdedwQwtLth7gzU3lbNl/hACBc8bEce/csaQnxfRpmyKy0RiT0d3zvOLKjtLS0i4fP1x1kN9+sPOkXxcaFEBMeDAx4cEMjAhmaEwYNfXNvLO5nLrjVmwLDhRSBkcwKi6KH84ZxZQUvc7KEzW1tPPRl5VcNjVJS115vEGRIdw428aNs20UVtXz1qYy3t5cToAbViH1imJPTU3t8og9NTWFgt/Mo6GljaaWdhqa2wkKFAaGBzMgPPikU+GMMdQ0tLCvuoF99gb2VjdQXN1Afslhrvjb59w9ZxR3f2eUloeH+XhnFY0t7TobRnmdUfFR3Dd/HPfOHUuAG+6p6hXFnpOTc8IYO0BERAQPPvggkaFBRIb2bjdEhNioUGKjQplhG/z147VNrfzmvQL++NEePtlVxaPXTGFkXJTL9kM5Z8nWA8RFhzJz+BCroyjVJ+4odfCSRcAyMzPJzc0lLS0NESEtLY3c3FyXnziNCQ/m0Wum8MT10yg51MjFf/qUF9cUY8V5CHWi+uY2Pt5ZxUXpQ33iLvJK9SevOHlqhcqjx/j5G1tZtdvOOWPiePjKSSQMcM8ZbfVt724p5yevbOH1O2af8FeWUv6kpydPveKI3QoJA8J4/pYZPHBZOuv31TDv8VX8a+tBq2P5rX9+cZChA8KYnjrI6ihKeTwt9lMQEW6clcb7Pz6btMER3PXSJpbvqLQ6lt+pbWpl1W47F08a5rYxSqW8mRZ7D4yIi+K1O2Yzbmg0v3xrKzX13nEXFV/x7x2VtLQ7dDaMUj2kxd5DoUGBPHbNFI42tZH99nY9oepGS7YeIHlQuF5foFQPabH3wmnDBnDP3DEsLajgrU3lVsfxC4cbWvhsTzULJiUibriwQylfoMXeS7efPYIZtkHc/14B5UearI7j85YWVNDmMDoMo1QvaLH3UmCA8MhVU3AYw72vfYHDoUMy/en9bQcZHhvJhMQBVkdRymtosfdB6pAI/mvBeNbsreG5z4utjuOzjjS2sKaohvnpQ3UYRqle0GLvo2tmpHD+uHiy//A3klJSdZ34fvDRl1W0OQzzJwy1OopSXkWLvY9EhIz2Airf/xMHyvZjjPl6nXgtd9dYVlDBsJgwJiX3bYlTpfyVFrsTfvfA/ZjWE+e0NzY2kp2dbVEi39HY0sbK3XbmTdBhGKV6S4vdCSdbJ/5kj6ueW7nLTnObg3k6DKNUr2mxOyE1NbVXj6ueW1pQweDIEGbYdG0YpXrL2ZtZPyAiW0Vki4h8KCJ+dSPKnJwcIiIiTngsNCycnJwcixL5hpY2Bx9/WcUFp8XrLQuV6gNnf2seNsZMMsZMAZYAv3ZBJq/xzXXiQwfGM+Gan3P99ddbHc2rfV5UTV1zG/PTdRhGqb5wqtiNMUeP+zAS8LurdTIzMykuLsbhcPD0++uoGXo6G4oPWx3Lqy0rqCAqNIgzRsZaHUUpr+T037kikiMi+4FMTnHELiJZIpIvIvl2u93Zl/VIV05PYXBkCLmr9lodxWu1Owz/3lHJeWPjTnrPWqXUqXVb7CKyXES2d/G2EMAYk22MSQHygLtPth1jTK4xJsMYkxEXF+e6PfAg4SGB3DgrjeVfVlJYVW91HK+0seQw1fUtOgyjlBO6LXZjzAXGmPQu3t79xlPzgCv6J6b3uGl2GqFBATzzmR6198XS7RWEBAVw3th4q6Mo5bWcnRUz+rgPFwI7nYvj/YZEhXLl9GTe3FSOvU5vyNEbxhiWFVRw9qhYokKDrI6jlNdydoz9oc5hma3AXOAnLsjk9b5/9gha2x08rwuE9UrBgaOUH2ling7DKOUUZ2fFXNE5LDPJGHOJMUbvPgEMj41k7vgEXlxbQmNLm9VxvMbS7RUECFxwWoLVUZTyanr1Rz/JOmcktU2tvLZhv9VRvMayggpmDh/C4MgQq6Mo5dW02PvJ9LRBTE8bxDOr99HW7rA6jscrrKpnT1W9zoZRygW02PtR1jkj2H+oiaUFFVZH8XjLOv+N5k7QYRilnKXF3o8uPC2B4bGR5K7aizF+d1Fur3xYUMHklIEMiwm3OopSXk+LvR8FBAjfP3s4W8tqWbfvkNVxPFb5kSa+KKvVOyUp5SJa7P3simnJDNFlBk7pw85hmHk6DKOUS2ix97Ow4EBumm3j451V7KmsszqOR1pWUMGYhChGxEVZHUUpn6DF7gY3zk4jLDiApz/Vo/ZvqqlvZv2+QzoMo5QLabG7weDIEC6bksSSrQdpaNYLlo63/MtKHAbmarEr5TJa7G5y5fRkGlva+WC7Tn083tLtFSQPCmdC4gCroyjlM7TY3WR62iBsQyJ4c2OZ1VE8Rt2xVlYX1jB/wlBExOo4SvkMLXY3EREun5bMmr01lB1utDqOR1ixy05Lu0OvNlXKxbTY3WjR1CQA3tqka6UBLNteQVx0KNNSB1kdRSmfosXuRimDI5g9YghvbSrz+ytRj7W2s2JXFXPHJxAQoMMwSrmSFrubXTE9meKaRjaW+PcNrz/dU01jS7sOwyjVD7TY3ey76UOJCAnkzU3+eRI1Ly8Pm83G3PRhHHjyVvatXWZ1JKV8jha7m0WGBjE/fShLvjjIsdZ2q+O4VV5eHllZWZSUlIAxtNZWcecdPyAvL8/qaEr5FJcUu4j8TESMiMS6Ynu+7sppydQ1t329VK2/yM7OprHxxBlBjY2NZGdnW5RIKd/kdLGLSAod9zstdT6Of5g1YghJA8N5089mx5SWdv0jcrLHlVJ944oj9seA+wD/nubRCwEBwuXTkvhsj53Ko8esjuM2qampvXpcKdU3ThW7iCwEyo0xX7goj9+4fFoyDgNvb/afo/acnBzCwk+8kUZERAQ5OTkWJVLKN3Vb7CKyXES2d/G2EPgV8OuevJCIZIlIvojk2+12Z3N7veGxkUxPG8SbG/1nTntmZiYLfng/QQPiEBHS0tLIzc0lMzPT6mhK+ZRui90Yc4ExJv2bb8BeYDjwhYgUA8nAJhHpcmKyMSbXGJNhjMmIi4tz5T54rSumJbOnqp5t5bVWR3ELYwwV8RlkPr4Eh8NBcXGxlrpS/aDPQzHGmG3GmHhjjM0YYwPKgGnGGP+a6uGEiycNIyQogDf8ZGGwHQePsv9Qk669rlQ/03nsFooJD2behKG898UBmtt8f077su0VBAhcMF5vgadUf3JZsXceuVe7anv+4oppSRxpbGXFziqro/S7pQUVzLANJjYq1OooSvk0PWK32Nmj44iPDuWNjb49O2avvZ7dlfW6NoxSbqDFbrHAAGHR1CQ+2VVFdX2z1XH6zbKCSgDm6fi6Uv1Oi90DXDE9mTaH4d0tB6yO0m+WFlQwOTmGxIHh3T9ZKeUULXYPMCYhmolJMby92Tdnxxw40sQX+48wT4dhlHILLXYPsWhqEtvLj7K7ss7qKC73YediZzrNUSn30GL3EJdOSSQwQHzytnlLCyoYkxDFiLgoq6Mo5Re02D1EbFQo546J490t5bQ7fGeJgZr6ZtbvO6RH60q5kRa7B7l8WhIHa4+xdm+N1VFcZmlBBQ6Djq8r5UZa7B7kgtMSiA4N8qnhmHc3H2BkXCTjhw2wOopSfkOL3YOEBQdy8aRhfLD9II0tbVbHcVrZ4UbWFx9i0dQkRMTqOEr5DS12D7NoahKNLe182HlBjzf7al7+wilJFidRyr9osXuYGbbBJA0M5y0vvwGHMYZ3NpeTkTaIlMERVsdRyq9osXuY42+bV+XFt83bcfAoe6rqWThVj9aVcjctdg+0aGoSDoNXLzHw7pYDBAUICyYOszqKUn5Hi90DjYiLYkrKQN7c5J1LDLQ7DO9uKee8sXEMigyxOo5SfkeL3UNdPi2JnRV17Dhw1OoovbZubw2VR5u5TIdhlLKEFruHWjApkeBA8cqFwd7eXE5UaBAXnKZ3SlLKClrsHmpwZAjnjY3n3S0HvGqJgWOt7SzdXsH89KGEBQdaHUcpv+RUsYvI/SJSLiJbOt8uclUw1XHbvKq6ZlYXes8dBz/6soq65jYW6TCMUpZxxRH7Y8aYKZ1v77tge6rTnHHxxIQH85YXnUR9Z0s58dGhzBoxxOooSvktHYrxYKFBHUsMLCuopL7Z85cYONLYwie7qljYuQSxUsoarij2u0Vkq4gsFpFBLtieOs4V05Jo6hy39nT/2naQ1najSwgoZbFui11ElovI9i7eFgJ/A0YCU4CDwCOn2E6WiOSLSL7dbnfZDvi6aamDSBsS4RWzY97ZXM7o+CgmJOpKjkpZKai7JxhjLujJhkTkaWAGh+axAAAJ6ElEQVTJKbaTC+QCZGRkeM80D4uJCIumJvHHj/ZwsLaJYTGeeTPo/Yca2VB8mJ/PG6srOSplMWdnxRx/vfgiYLtzcVRXLp+aDMAr6/dbnOTk3vuiY/mDSycnWpxEKeXsGPvvRWSbiGwF5gA/dUEm9Q2pQyL4zth48taVcKy13eo432KM4e3N5cyw6UqOSnkCp4rdGHOjMWaiMWaSMeZSY8xBVwVTJ7r1rOFU17d8fWTsSQoOHKWwql6XEFDKQ+h0Ry9xxsghjBsazeLP9mGMZ52ieHdLOcGBwsW6kqNSHkGL3UuICLeeOZydFXWsKfKcm123OwzvfXGA88bGMzBCV3JUyhNosXuRS6ckMiQyhMWr91kd5WsfFlRQebSZK6bpMIxSnkKL3YuEBQeSOSuNj3ZWsa+6weo4GGN44pNChsdGcuH4oVbHUUp10mL3MjfMSiU4IIBnPeCofdWearaXH+XOc0fqEgJKeRAtdi8THx3GJZMTeT2/jNrGVkuzPLGikGExYTobRikPo8XuhW47azhNre28sqHUsgwbig+xft8hss4ZQUiQ/hgp5Un0N9ILjU8cwOwRQ3j+82La2h2WZPjrikIGR4Zw7YxUS15fKXVyWuxe6tazhnOg9hhLC9y/6mPBgVpW7LJz21nDCQ/RuyQp5Wm02L3U+ePiSRsSweLP3H8S9a8riogODeKGWWluf22lVPe02L1UQIBwyxk2NpUeYXPpYbe9bpG9nve3H+TG2WnEhAe77XWVUj2nxe7FrspIITosiGfceNT+5CdFhAQGcOtZw932mkqp3tFi92KRoUFcOyOFD7ZXcOBIU7+/XvmRJt7eXM51p6cSGxXa76+nlOobLXYvd/MZNowxPL+muN9f6+lVewG4/ZwR/f5aSqm+02L3csmDIpifPpSX15Vy9Fj/XbBUXd/My+tLWTQ1iaSBnnkXJ6VUBy12H3DnuaNoaGnnV29t67clfRd/to+Wdgd3nDeyX7avlHIdLXYfMDE5hnsuHMOSrQd5ZYPrb59X29TKi2tKuCh9GCPjoly+faWUa2mx+4g7zx3JWaNiuf+9AnZV1Ll02/9YW0Jdcxs/nKNH60p5A6eLXUR+JCI7RaRARH7vilCq9wIChEevmUx0WBB3v7SJphbX3Bt1V0UdT60s4ryxcUxIjHHJNpVS/cupYheROcBCYLIxZgLwB5ekUn0SHx3GY9dModBez/3vFTi9vZKaBm54Zh3hIYE8sDDdBQmVUu7g7BH7ncBDxphmAGNMlfORlDPOHh3HneeO5NX8/by7pbzP2zlY20Tm39fR1u7gH7fNJGVwhAtTKqX6k7PFPgY4W0TWichKEZlxsieKSJaI5ItIvt1ud/Jl1ancc+EY4io3cNV50wgICMBms5GXl9fjr6+pb+aGv6+jtrGVF26dyeiE6H5Mq5RytaDuniAiy4Gu7nuW3fn1g4FZwAzgNREZYbqYc2eMyQVyATIyMvpnTp4C4NVXXqbg1d/T2tRxNWpJSQlZWVkAZGZmnvJrjx5r5abF6yk/0sQLt85kYrKOqyvlbcSZec8ishT4nTFmRefHRcAsY8wpD8kzMjJMfn5+n19XnZrNZqOkpORbj6elpVFcXHzSr2tqaeemxevYsv8IT9+UwXlj4/sxpVKqt0RkozEmo7vnOTsU8w4wp/MFxwAhQLWT21ROKi3t+s5KJ3scoLmtnR/8YyMbSw7zx2unaqkr5cW6HYrpxmJgsYhsB1qAm7sahlHulZqa2uURe/CAOK7NXcPo+GhGJ0QxKj6K0fHRDIoI5j9e2cKq3XZ+f+UkLpo4zILUSilXcarYjTEtwA0uyqJcJCcnh6ysLBobG79+LCQ0nAtv/g9a2hy8s6WcumNtX38uPDiQptZ2fr1gPFdnpFgRWSnlQs4esSsP9NUJ0uzsbEpLS0lNTSUnJ+frx40xVNU1s6eynsKqOvZU1TM5ZaCWulI+wqmTp32lJ0+VUqr33HXyVCmllIfRYldKKR+jxa6UUj5Gi10ppXyMFrtSSvkYLXallPIxWuxKKeVjtNiVUsrHWHKBkojYgW8vZtIzsfjfQmO6z/5B99k/OLPPacaYuO6eZEmxO0NE8nty5ZUv0X32D7rP/sEd+6xDMUop5WO02JVSysd4Y7HnWh3AArrP/kH32T/0+z573Ri7UkqpU/PGI3allFKn4LHFLiLzRWSXiBSKyC+7+HyoiLza+fl1ImJzf0rX6sE+3yMiO0Rkq4h8JCJpVuR0pe72+bjnXSEiRkS8egZFT/ZXRK7u/D4XiMhL7s7oaj34uU4VkRUisrnzZ/siK3K6kogsFpGqztuGdvV5EZE/df6bbBWRaS4NYIzxuDcgECgCRtBxg+wvgPHfeM4PgSc7378WeNXq3G7Y5zlAROf7d/rDPnc+LxpYBawFMqzO3c/f49HAZmBQ58fxVud2wz7nAnd2vj8eKLY6twv2+xxgGrD9JJ+/CPgAEGAWsM6Vr++pR+ynA4XGmL2m476qrwALv/GchcDzne+/AZwvIuLGjK7W7T4bY1YYY766kelaINnNGV2tJ99ngAeA3wHH3BmuH/Rkf28HnjDGHAYwxlS5OaOr9WSfDTCg8/0Y4IAb8/ULY8wq4NApnrIQeMF0WAsMFBGX3UXeU4s9Cdh/3MdlnY91+RxjTBtQCwxxS7r+0ZN9Pt5tdPyP78263efOP1FTjDH/cmewftKT7/EYYIyIrBaRtSIy323p+kdP9vl+4AYRKQPeB37knmiW6u3ve6/ozay9kIjcAGQA51qdpT+JSADwKPA9i6O4UxAdwzHn0fEX2SoRmWiMOWJpqv51HfCcMeYREZkNvCgi6cYYh9XBvJWnHrGXAynHfZzc+ViXzxGRIDr+hKtxS7r+0ZN9RkQuALKBS40xzW7K1l+62+doIB34RESK6RiLfM+LT6D25HtcBrxnjGk1xuwDdtNR9N6qJ/t8G/AagDFmDRBGx3oqvqxHv+995anFvgEYLSLDRSSEjpOj733jOe8BN3e+fyXwsek8K+Glut1nEZkKPEVHqXv72Ct0s8/GmFpjTKwxxmaMsdFxXuFSY0y+NXGd1pOf63foOFpHRGLpGJrZ686QLtaTfS4FzgcQkdPoKHa7W1O633vATZ2zY2YBtcaYgy7butVnj09xVvkiOo5WioDszsf+m45fbOj45r8OFALrgRFWZ3bDPi8HKoEtnW/vWZ25v/f5G8/9BC+eFdPD77HQMfy0A9gGXGt1Zjfs83hgNR0zZrYAc63O7IJ9fhk4CLTS8VfYbcAdwB3HfZ+f6Pw32ebqn2u98lQppXyMpw7FKKWU6iMtdqWU8jFa7Eop5WO02JVSysdosSullI/RYldKKR+jxa6UUj5Gi10ppXzM/wIfZpbUqYQi2wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "c, d = 0, 1\n",
    "xt = np.linspace(c, d, 8, False)\n",
    "x = np.array([-2.2, -2.8, -6.1, -3.9, 0, 1.1, -0.6, -1.1])\n",
    "y = np.fft.fft(x) / np.sqrt(x.size)\n",
    "a = np.array(y.real)\n",
    "b = np.array(y.imag)\n",
    "t = np.linspace(c, d, 40)\n",
    "data = [even_trigonometric_interpolant(i, a, b, c, d) for i in t]\n",
    "plt.plot(t, data)\n",
    "plt.plot(xt, x, 'o', color='k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 10.3 The FFT And Signal Processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Orthogonal Function Interpolation Theorem\n",
    "\n",
    "Let f_0(t),\\cdots,\\f_{n-1}(t) be functions of t and $t_0,\\cdots,t_{n-1}$ be real numbers. Assume that the $n \\times n$ matrix\n",
    "\n",
    "### $A = \n",
    "\\begin{bmatrix}\n",
    "f_0(t_0) & f_0(t_1) & \\cdots & f_0(t_{n-1})  \\\\ \n",
    "f_1(t_0) & f_1(t_1) & \\cdots & f_1(t_{n-1})  \\\\ \n",
    "\\vdots & \\vdots &  & \\vdots \\\\ \n",
    "f_{n-1}(t_0) & f_{n-1}(t_1) & \\cdots & f_{n-1}(t_{n-1}) \n",
    "\\end{bmatrix} $\n",
    "\n",
    "is a real $n \\times n$ otrhogonal matrix. If $y = Ax$, the function\n",
    "\n",
    "### $F(t) = \\sum_{k = 0}^{n-1}y_kf_k(t)$\n",
    "\n",
    "interpolates $(t_0, x_0), \\cdots, (t_{n-1}, x_{n-1})$, that is $F(t_j) = x_j$ for $j = 0,\\cdots,n-1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lemma\n",
    "\n",
    "Let $n \\ge 1$ and k, l be integers. Then\n",
    "\n",
    "### $\\sum_{j = 0}^{n-1}\\cos{\\frac{2\\pi jk}{n}}\\cos{\\frac{2\\pi jl}{n}} =  \\left\\{\\begin{matrix}\n",
    "n & \\text{if both (k-l)/n and (k+l)/n are integers} \\\\ \n",
    "\\frac{n}{2} & \\text{if exactly one of (k-l)/n and (k+l)/n is an integer} \\\\ \n",
    "0 & \\text{if neither is an integer}\n",
    "\\end{matrix}\\right.$\n",
    "\n",
    "### $\\sum_{j = 0}^{n-1}\\cos{\\frac{2\\pi jk}{n}}\\sin{\\frac{2\\pi jl}{n}} = 0$\n",
    "\n",
    "### $\\sum_{j = 0}^{n-1}\\sin{\\frac{2\\pi jk}{n}}\\sin{\\frac{2\\pi jl}{n}} =  \\left\\{\\begin{matrix}\n",
    "0 & \\text{if both (k-l)/n and (k+l)/n are integers} \\\\ \n",
    "\\frac{n}{2} & \\text{if (k-l)/n is an integer and (k+l)/n is not} \\\\ \n",
    "-\\frac{n}{2} & \\text{if (k-l)/n is not and (k+l)/n is an integer} \\\\ \n",
    "0 & \\text{if neither is an integer}\n",
    "\\end{matrix}\\right.$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Orthogonal Function Least Squares Approximation Theorem\n",
    "\n",
    "Let $m \\le n$ be integers, and assume that data $(t_0,x_0),\\cdots,(t_{n-1},x_{n-1})$ are given. Set $y = Ax$, where A is an orthogonal matrix of form like Orthogonal Function Interpolation Theorem show. Then the interpolating polynomial for basis functions $f_0(t),\\cdots,f_{n-1}(t)$ is\n",
    "\n",
    "### $F_n(t) = \\sum_{k = 0}^{n - 1}y_kf_k(t)$\n",
    "\n",
    "and the best least squares approximation, using only the functions $f_0,\\cdots,f_{m-1}$, is\n",
    "\n",
    "### $F_m(t) = \\sum_{k = 0}^{m - 1}y_kf_k(t)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Corollary\n",
    "\n",
    "Let $[c,d]$ be an interval, let $m < n$ be even positive integers, $x = (x_0,\\cdots,x_{n-1})$ a vector of n real numbers, and let $t_j = c + j(d - c)/n$ for $j=0,\\cdots,n-1$. Let ${a_0,a_1,b_1,a_2,b_2,\\cdots,a_{n/2-1},b_{n/2-1},a_{n/2}} = F_nx$ be the interpolating coefficients for x so that\n",
    "\n",
    "### $x_j = P_n(t_j) = \\frac{a_0}{\\sqrt{n}} + \\frac{2}{\\sqrt{n}}\\sum_{k=1}^{n/2-1}(a_k\\cos{\\frac{2k\\pi (t_j-c)}{d-c}}-b_k\\sin{\\frac{2k\\pi (t_j-c)}{d-c}})+\\frac{a_{n/2}}{\\sqrt{n}}cos{\\frac{n\\pi (t_j-c)}{d-c}}$\n",
    "\n",
    "for $j = 0,\\cdots,n-1$. Then\n",
    "\n",
    "### $P_m(t) = \\frac{a_0}{\\sqrt{n}} + \\frac{2}{\\sqrt{n}}\\sum_{k=1}^{\\frac{m}{2}-1}(a_k\\cos{\\frac{2k\\pi (t-c)}{d-c}}-b_k\\sin{\\frac{2k\\pi (t-c)}{d-c}})+\\frac{a_{\\frac{m}{2}}}{\\sqrt{n}}cos{\\frac{n\\pi (t-c)}{d-c}}$\n",
    "\n",
    "is the best least squares fit of order m to the data $(t_j,x_j)$ for $j = 0,\\cdots,n - 1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def least_squares_even_trigonmetric(m, t, a, b, c, d):\n",
    "     return a[0] / np.sqrt(a.size) + 2 / np.sqrt(a.size) * \\\n",
    "        np.sum([a[k] * np.cos(2 * k * np.pi * (t - c) / (d - c)) - \\\n",
    "        b[k] * np.sin(2 * k * np.pi * (t - c) / (d - c)) for k in range(1, int(m / 2))]) + \\\n",
    "        a[int(m / 2)] / np.sqrt(a.size) * np.cos(a.size * np.pi * (t - c) / (d - c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Fit the temperature data by least squares trigonometric functions of orders 4 and 6 where $x = [-2.2, -2.8, -6.1, -3.9, 0.0, 1.1, -0.6, -1.1]$ on the interval $[0,1]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xdc1dX/wPHX57K3gILKHi5EBHEPFPdeaaakqZkts12WlX0ts8ysn5mVltuypWZuzL3FrThAZCqCqOx1L+f3x02TRFkXLuM8Hw8eyed+Pue87w3fHs7nfN5HEUIgSZIk1RwqfQcgSZIk6ZZM7JIkSTWMTOySJEk1jEzskiRJNYxM7JIkSTWMTOySJEk1jEzskiRJNYxM7JIkSTWMTOySJEk1jKE+Oq1bt65wd3cv07WZmZlYWFjoNiAdkHGVjoyrdGRcpVNV4yqv48eP3xRC1Cv2RCFEpX8FBgaKstq1a1eZr61IMq7SkXGVjoyrdKpqXOUFhIkS5Fg5FSNJklTDyMQuSZJUw8jELkmSVMPo5eZpUfLz84mPjycnJ+eR59nY2HDhwoVKiqrkKjIuU1NTnJ2dMTIyqpD2JUmqWapMYo+Pj8fKygp3d3cURXnoeenp6VhZWVViZCVTUXEJIUhJSSE+Ph4PDw+dty9JUs1TZaZicnJysLe3f2RSr40URcHe3r7Y32QkSZLuqjKJHZBJ/SHk5yJJUmlUmakYSar1hICsFEi5AimRkH0bWj4BFnX1HZlUzcjEXk7Lli0jLCyM2bNnl6udtLQ0fHx8GDp0KAsWLNBRdFKVJQREhML1U9okfvcrJ7XweXvnQPf3QXjqJ06pWpKJvRTuPdWlKvsMllqtxtDwwY/9/fffJygoqDzhSdVFgQa2vAXHftB+b+MC9l7QYiTYeYG9t/Z7dQ5snQab3yDQ0gu8vweXNvqNXaoWqtQcu77NmzcPX19ffH19+eqrrwCIjo6mSZMmjBs3Dl9fX+Li4li6dCmNGzembdu2HDhw4N71ycnJPPbYY7Rp04Y2bdrce+3DDz9k7NixdOrUibFjxz7Q7/Hjx7lx4wa9e/eunDcq6U9+Nvw6TpvUO74E0xPh1XMw7k8Y8AV0eAEa99YmdsfmMG4DjFiCcd5t+LEn/DkFMm/q+11IVVyVHLH/76/zhF9LK/I1jUaDgYFBqdv0aWjNjEHNH/r68ePHWbp0KUeOHEEIQbt27ejatSu2trZERESwfPly2rdvz/Xr15kxYwbHjx/HxsaG4OBgAgICAHj55Zd59dVX6dy5M7GxsfTp0+fe2vbw8HD279+PmZlZoX4LCgp4/fXXWbVqFTt27Cj1+5Kqkaxb8NMoiD8G/eZAu2eLv0ZRwPcxjiaa0aXgEBxeCBf+gh7vQ+AEUJX+74JU81XJxK4P+/fvZ9iwYfcqwg0fPpx9+/YxePBg3NzcaN++PQBHjhyhW7du1KunLbA2atQoLl++DMCOHTsIDw+/12ZaWhoZGRkADB48+IGkDrBw4UL69++Ps7Nzhb4/Sc9uR8Oqx+BOHDy+HHyGlOpyjaE5dPsI/ENg8xuw6XU49TOMXQem1hUTs1RtVcnE/qiRtT4eUCpp+c+CggIOHz6Mqalpids4dOgQ+/btY+HChWRkZJCXl4elpSWffvppuWKWqpBrJ2H146DJg6c2kF6/OWvOLOb3y79jbmSOVx0vvGy88KzjiZeNF27WbhgZPOQpY4em8NRfcOYXWP8CrJ0MT/wE5bjvU1p3svK4lJhOoJttpfUplU6VTOz60KVLF8aPH8+0adMQQrBu3TpWrlz5wHnt2rXj5ZdfJiUlBWtra3777TdatmwJQO/evfn666958803ATh16hT+/v6P7Hf16tX3/nx3hY1M6jVIxA7tnLq5PbfH/MyqpEP8vP810vPT6dCgAyYGJoSnhLM9ejsCAYCBYoCrtStNbJvgEuHCvI/nERsbi6urK7NmzSIkJES7DDInDba8CbtmaadmKtCNtBy2n09k2/kbHIpKQVMgaN7QmpFumgrtVyobmdj/0apVK8aPH0/btm0BmDRpEgEBAURHRxc6r0GDBnz44Yd06NCBOnXqFErc8+fP58UXX8TPzw+1Wk1QUBDfffddZb4NqSo59TP8+SJJjs1Y3qIXv+1+gRx1Dj3dejKpxSR87H3unZqjziE6LZord65w5c4VolKjWP/reqJ+iKIgrwCAmJgYJk+eDKBN7m2fgRtnYd9c7Y1W3+E6DT/6Zibbziey7XwiJ2LvAOBZ14LJQZ642pkzL/QyMw/nct3oAq/0bIyZsZzvryoUbe32cjaiKEuAgUCSEMK3uPNbt24twsLCCh27cOECzZo1K7av2lYr5q6Sfj7/tXv3brp166b7gMqpxsd1/TTXlvTiR2dv1qlyKBAF9PPox6QWk/Cq41WiJpxdnUmIS3jguJub278DDnUuLB8E18/A09uhgV+5Q09Ky2HyyuOcitMmc18na/o2r0+f5vXxdrC89yR0anY+U3/cyZ54Na525swe3oJO3lXjYaqq+vNVXoqiHBdCtC7uPF2N2JcBC4AVOmpPkqovdR6X/nyOCQ0dyFayGeo9lIm+E3GxcilVM9firxV5PDY29t9vDE1g1CpY1A3WjIFndoFl8TunPcztzDye/PEICbezeW9AM/o0r4+LnXmR59qYGTHB14Tn+rXm3XVnCfnhCCMDnZk+oBl1zI3LHINUfjq54yKE2Avc0kVbklTdRe/6kMnGqZib1mHD0A3M6DCj1EkdwNXVtcjjDg0dCh+wdIAnVkNmsnY+X51XlrDJyFUzfulRolOyWPxUayZ18XxoUr9fBy97trzchRe6ebH2ZAI95+1hy9nrZYpB0g2dTMUAKIriDmx82FSMoiiTgckAjo6OgWvWrCn0uo2NDd7e3sX2U9Z17BWtouOKjIwkNTW1+BP/IyMjA0tLywqIqHxqaly5qWf44ua3ZBkYM8XpHeob1S9zWzt27GDu3Lnk5ubeO6YyVtFwfENC+oXQ26Y3KuXfsZnDjT34XJhHQsO+RDR+vlR95WkE847nEHG7gCkBJgQ4lOyX+f9+XrFpGpaey+NqWgFvtzGlmb1+/q5W1Z+v8goODi7RVIzONqgG3IFzJTm3qM2sw8PDS7SZa1paWonOq2wVHVdJP5//qqqb+tbEuFIyboiBS/xE+yXNxfn4QzqJZ9WqVcLNzU0oiiLc3NzEkuVLxLS904TvMl/xfOjzIledW/iC7e8LMcNaiKM/lLiPPLVGTFx6VLhP2yjWn4wvVXxFfV5ZuWrRdc5O0eWznSIzN79U7elKVf35Ki/kZtaSVHnS89J5bsNIrqNmQfPJ+Di110m7ISEhREdHs3PnTqKjo5kwbgKfdP6E6e2msy9hH/879L+7AyutHjOgUW9tLZroAw9v+B+aAsFrv57m74tJfDzUlyH+TuWO2czYgM8e8yP2VhZztl4qd3tS6cnELknllK3OZsrWiUTkpjDPrDGBbadWaH+KovBE0yd4oeULbLiygWXnl/37osoAHvsBbD3g94nate4PIYTgvfXn+Ov0Nab1a0pIOzedxdjO056nOrix/FA0x6Ll7bfKppPErijKz8AhoImiKPGKojyti3arg2XLljFlypRytREbG0vv3r1p1qwZPj4+D6ydl6qufE0+r+16lZO3LjI7NY+gQYsrre9nWz5Lb7fefHn8S/bE7fn3BVMbGP49ZNzQPrxUBCEEn265yM9HY3mhmxfPdS3ZEszSeKtvU5xtzXjr9zNk58kHmSqTrlbFjBZCNBBCGAkhnIUQP+qi3apGCEFBQUG52lCr1Q8cGzduHG+++SYXLlzg6NGjODg4FHGlVNVoCjS8u/9d9l87wAcpKfTtOadSN8VQKSo+7vwxzeyb8dbet4i4HfHvi06B0GYSHF2kLWnwH9/uucL3e6MY296NN/s0qZD4LEwM+Wy4H1dvZjIvVE7JVCY5FXMffZTtDQ8PR61W06tXLwAsLS0xNy9+iZmkf/OOz2Nr9FZeu53GCJee0HxopcdgZmjG/OD5WBhZ8NLOl7iVc9+0R4/3waIebHxVWwP+Hwcib/L5tksMbtmQ/w1uXqFbL3b0rktIO1d+2H+V4zG3K6wfqbCqWVJgyzRIPFvkS2YaNRiUIez6LaDfw2uw6Kts7+XLl6lTpw7Dhw/n6tWr9OzZk08//bRKLumU/nUy6SQrwlcwSmPKhNwM6P+F3mJxtHDk/4L/jwnbJvDqrlf5ofcP2iJipjbQ5xP442kIWwJtnyElI5dXfzmFZ10LPn2sBSpVxe+n+07/Zuy+lMxbv59m09QumBrJn+2KJkfs/7i/bK+lpeW9sr3AQ8v2GhsbM2rUqHtt7NixgylTpuDv78/gwYNLVLZXrVazb98+5s6dy7Fjx4iKimLZsmUV/4alMsvV5DLj4AwaGlryWlwEDJhbrqc9daFFvRbM7DiTE0kn+PjIx/+ulPF9DDyD4e+ZiLRrvP3HGe5k5fP16FaYG1fOuM7SxJDZw1twJTmTr3ZEFH+BVG5Vc8T+iJF1dg0r2+vs7Iy/vz+ento9LYcOHcrhw4d5+ulac/+52ll0ZhFXU6/yXUom5t69oblui2+VVX/P/kTeiWTx2cV41/FmrM9Y7UYdA76AhR2I/ukVdkQ/xYxBPvg0rNwa7kGN6zGqtQuL9l6hn299WrrUqdT+axs5Yv9Hly5dWL9+PVlZWWRmZrJu3Tq6dOnywHnt2rVjz549pKSkkJ+fz2+//Xbvtbtle+86depUsf22adOGO3fukJycDMDOnTvx8fEp5ipJXy7dusSSs0sYZO5Gp7QU6PGBNnlWEVMCptDDtQdzw+ayL177Gyf2XiQHTMEjcRsvucUwvqO7XmKbPrAZjtamvPHbaXLVcpVMRZKJ/R/3l+1t167dvbK9/3V/2d5OnToVqrg4f/58wsLC8PPzw8fHp0Qlew0MDJg7dy49evSgRYsWCCF45plndPreJN3QFGj48OCHWBtb8taVU+AzFOoXW8y0UqkUFZ90/oTGto2Ztm8aiZmJZOdpGHepIzE05JXc71HUOXqJzdrUiE+GtyAiKYOv/47USwy1RdWcitGT1157jddee63QMXd3d86dO1fo2IQJE5gwYUKhY+np6dStW5dffvnlgXY//PDDR/bbq1cvzpw5U7agpUqz6sIqzqWcY45tG+rknIVu7+g7pCKZG5nzRdcvGPnXSN7Z9w71s17hQnIuqQM+w+3vsbBvHnSfrpfYgps4MCzAiUX7ohjbwQ1H6wenLaXykyN2SSqBuPQ4FpxcQNf67el7dgu0GKHdpq6KcrV2ZXr76YTdCOOPKyt5NsgTvy6DwW8U7P8Ski/rLbZXezZGUyD4bs8VvcVQ08nELknFEEIw89BMDFQGvKex0k5ldJ2m77CKFWjXEzICMKkXSu9W/1SI7P0xGJvDptdAR5VdS8vV3pxhAU78dCSWpHT9TAvVdDKxS1Ix1keu5/D1w7zqM5H6J1aD3xNQt/gS0/qkKRC8+utpCpKH4WDuyPQD08jIy9DWbu/5IUTv026IrSdTgr3J1xSweG+U3mKoyWRil6RHuJl9k7lhc2nl0IqRCZehQA1d39J3WMX6emcER6/eYuagNnzRbQ7XM68z68g/dWNajQen1hD6AeSm6yU+97oWDPV3YuXhGG5m5BZ/gVQqMrFL0iPMPjKbbHU2M1o8h+rEcvAPATsPfYf1SKHhN/hqRwTDA5wY3soJfwd/nmv5HBujNvLXlb9ApYJ+c7RFwvZ+rrc4X+zuTZ66gMX75Khd12Ril6SH2B23m+0x23mu5XN4nvxFOycd9Ia+w3qkiBvpvPrLKfycbfhkeIt7dWCeafEMrRxaMevILOLS48A5UPuP1KGFkKKfm5he9SwZ1LIhKw/FcCuzbNv5SUWTib2cdFG296233qJ58+Y0a9aMqVOnFt44QdILTYGGr45/hYeNBxOcesDJlRD4FNQpeh/SqiA1K59nVoRhamTA92MDC9VkMVQZMrvLbFSomLZ3GvkF+dpNOQxNYdu7eov5pe7eZOdr+EGO2nVKJvZSqIiyvQcPHuTAgQOcOXOGc+fOcezYMfbs2fOQq6XKsiV6C1dSr/Ci/4sY7Z8HigF0eV3fYT2UWlPAlJ9PkHAnm+/HtqKBzYN1iRpaNuSDjh9w5uYZvjv9HVg5Qtc34fJWiNihh6jB28GKAS0asPxgNHey5KhdV2Riv48+yvYqikJOTg55eXnk5uaSn5+Po6Nj5b1p6QH5Bfl8e+pbmtg2oZelJ5z6GVpPBOuG+g7toeZsu8S+iJt8PNSXQDe7h57X170vQ72HsvjMYo4lHoN2z4OdF2ydBmr9JNaXujciM0/Dkv1X9dJ/TVQlnzz97OhnXLx1scjXNBpNmUraNrVryttt337o6/oq29uhQweCg4Np0KABQgimTJlSqEyBVPk2RG4gNj2Wr7t/jWrvXDAwhs6v6jush1p3Mp5Fe6N4qoMbo9oUP1X0Ttt3OJl0krf3vs2vg36lbt/Z8NPj2k05OpZvWrEsmtS3on+L+iw9EM3TnT2xMTeq9BhqGl1tjddXUZRLiqJEKopS9Z/cKIK+yvZGRkZy4cIF4uPjSUhIYOfOnff6lSrP6tWrcXd3p3v37jzZ6UmszlrR1aQBnP0V2k7STltUQWfi7/D2H2dp72nHewNLVjzubsmBtLw0pu2dhsa7J3j3gj2fQUZSBUdctCnBjUjPVbP0oBy160K5R+yKohgA3wC9gHjgmKIoG4QQ4WVt81Ej6/QaVrZ33bp1tG/fHktLSwD69evHoUOHiqwsKVWM1atXM3nyZLKysgDIuZnDiYUn+EnzIiF1TaHTK3qOsGhJ6TlMXnGcepYmfDOmFUYGJR+nNbFrwnvt3+P9A+/zzalvmNp3NixsD3/PhCELKjDqovk0tKa3jyNL9l9lYmcPrE3lqL08dDFibwtECiGihBB5wBpgiA7arVT6Ktvr6urKnj17UKvV5Ofns2fPHjkVU8mmT59+L6nflZudy/SluyDgyUrdx7SkctUanl91gtTsfBaNC8Te0qTUbQz1HsrwRsNZfHYxe3OuQ7vn4OSqIvdIrQxTezQiLUfN8gPReum/JtHFHLsTEHff9/FAu/+epCjKZGAygKOjI7t37y70uo2NDenpxT8Fp9FoSnReaTVq1IjRo0fTunVrQLvBtLe3NzExMRQUFNzr09LSkmnTptGuXTtsbGzw8/MjLy8PjUbDJ598wuuvv46vry9qtZpOnTrx1VdfkZubi5GRUZFx9+nTh23bttG8uXbvyZ49e9KtW7cHzs3JyXngMyuJjIyMMl1X0apSXLGxsUUfTy3gsNKKnCoQ5/2f153cAr4/ncuFWwW84G9C8uWT7C5jTa9OBZ04anSUN3e9yTv1ptDfyJrsNc9xMuDTEtWZ1/X/R/96Bny3+zLeIh4zw7LXua9KP196IYQo1xcwAvjhvu/HAgsedU1gYKD4r/Dw8AeOFSUtLa1E51W2io6rpJ/Pf+3atUu3gehIVYrLzc1NAA98udU113do99z9vA5EJovAj0JFk/c2i9/C4nTSdmxqrOiwuoMY9dcokXtsiRAzrIU4/Uup4tKV03G3hdvbG8WCnRHlaqcq/XzpEhAmSpCXdTEVkwC43Pe98z/HJKlamDVrFsamxoWOmRvBrPfe1FNEDyoQgvl/R/DkD0ewMTPkzxc7MyLQWSdtu1i78HHnjzmfcp45+fHQMOCfOjIZOmm/NPyc6xDcpB4/7r9KVp66+AukIukisR8DGimK4qEoijHwBLBBB+1KUqUYOGIgrhNdsXKwQlEU3GyNWDS2GSFTZ+g7NABuZuTyRVgO80IvM7hlQzZM6UyT+rpdQNDdtTsTmk/gl8u/sjFgOKRfh92zddpHSb3UoxG3MvNYfbjoKTKpeOVO7EIINTAF2AZcAH4VQpwvY1vlDadGkp9LxVp+fjmmbU05duEYp3+dTfRUM0Je/ahK7GV6JCqFAfP3cel2AbOHt+DLUf5YmFTM4ydTW00l0DGQmZdXEdlyBBxeCAnHK6SvR2nlaktn77p8vzeKnHy5N2pZ6GQduxBisxCisRDCSwgxqyxtmJqakpKSIpPYfwghSElJKXIJpVR+t3JuserCKvq696WJXRNc4jaAjQs00+/CrtSsfL7ZFcnoxYcxNzbkg/amjG7req+oV0UwVBnyedDnmBua86pIJNPKEf58SS9PpL7U3ZubGbmsOSpH7WVRZZ48dXZ2Jj4+nuTk5Eeel5OTUyWTXEXGZWpqirOzbuZTpcKWnF1CriaX5/2fh2snqZN6TrvLkEHp/2oIIchVF5CTryEnX/tfUyMDbC2MMDF8+NPSmgLB5RvpnIy9w4nY25yMvc2V5EwABvo1YPbwFhw/fOCh1+tSPfN6fN71cyZtn8RUj6YsOLMbswNfVXoN+nae9rR1t+O7PVGMbuf6yM9PelCVSexGRkZ4eBRf53r37t33HuGvSqpqXNLDJWUlsebSGgZ6DsTTxhN2PIPawAzDVuOKvTY2JYvPtl0kLPrWvSSeq354gThLE0NsLYywszDBzlz7XytTQyKS0jkdl0pGrvZGoZ2FMa1c6zC8lTNt3O1o425boaP0orSp34aPO33MewfeY4qnDwv2zsGs2SBwqNznK17q4c3YH4/y+/F4Qtq5VWrf1V2VSeySVNlWX1hNfkE+z/k9B6kJcH4t1xv2x8XU5qHXZOSq+WZXJD/uu4qBSqFfi/pYmRhiamSAiZEBpkYqTA0NtN8bqshRa7idmUdKZh63M/O4lZVPckYul29kcCcrD/e6FgwLcKKVWx0CXGxxszev9ERelEFegwB4b/97THGsx4I/X8Ts6VBQVd7IubN3Xfxd6vDt7is83tqlVE/W1nYysUu1UlZ+Fr9f/p0erj1wsXbRLu8TBSQ4DSy0dveuggLB7yfi+XzbJZLTcxke4MRbfZtS36bqTQvqyr/JfTpTsuNYcGgBZp1errT+FUVhag9vJi4LY93JBB5vXdT/GakoMrFLtdKGKxtIy0tjrM9Y7XrtsGXQbDA5Zg8W+zoWfYv//XWecwlpBLjWYdHYQAJcbSs/aD34N7m/y5RzC1nQqA9mDk0rrf/gJg74OlmzcFckwwOcMJSj9hKRn5JU6xSIAlZfWI2vvS/+9fy19VFyU6FD4ZK1yem5TPnpBCO/O0RKRh7/94Q/a5/vWGuS+l2DvAbxceCbhJkYMWXLWLLzs4q/SEcURWFKcCOiU7LYeOZ6pfVb3cnELtU6+xP2E50WzVifsSiiQLte26UduLS5d06+poBnV4YRGn6Dl3s04u/XuzLE36lKzH/rwyDfcXzs1Jcwkc2UP0eSrc6utL57+zjSxNGKBbsiKSiQy6FLQiZ2qdZZEb4CB3MHern3goub4E4MdHix0Dlf7bjMidg7zB3Zkld7NcbcWM5aDuoxh49FXcIyYug3rScuri50794dd3d3Vq9eXWH9qlQKU7p7E5mUwZZziRXWT00iE7tUq1y6dYkj148wuulojFRGcOgbqOMGTQfeO+dA5E0W7r7CE21cGNSy6m6HV+lUKgYN/pGeodfZ+/Uh4uPiEUIQExPD5MmTKzS592/RAM96Fny9M0KO2ktAJnapVll1YRVmhmaMbDwS4o9D3GFo//y9ZXxpuYJXfjmFVz1LZgxqrudoqyB7L/74WyDyCifXrKwspk+fXmHdGqgUpgR7czExnR0XblRYPzWFTOxSrXEz+yabojYx2GswNiY22j0+ja20m2mgXdK4+Gwuadn5LBgTgJmxfNqxKLFJd4o8HhMbQ1YF3lgd3LIhbvbmfL0zUpYeKYZM7FKt8eulX8kvyCekWQhk3oTza6HlE2CirZT4w/4ozt7U8P5AH5rWt9ZztFWXq2vRG2Yb2Rkx5M8h7IzdWSH9GhqoeKGbF2cTUtl1ST97s1YXMrFLtUKuJpdfLv1CkHMQHjYecGIFaPKgzSQATsXdYc7WS7R2NCCkXdGJS9KaNWsW5ubmhY6Zm5kxY+YMLI0seXnXyzwb+izrI9dzM/umTvseFuCMR10LPtp4QVZ+fAR5q1+qFTZHbeZWzi3tA0kFGghbCu5dwKEpaTn5TP35JI7WpkzwVWrtksaSCgkJAbR7xcbGxuBqrWLWcFdCJr1CvqExq8NXsyJ8BQevHQSguX1zgpyDCHIOwsfeB5VS9vGksaGKmUOaM/bHo3y/J4qXezbSyXuqaeSIXarxhBCsvLCSxraNaVe/HVzeBqmx0GYSQgjeXXuWhDvZzB/tj4WRTOolERISQnR0NDt37iL6yEZC3G7A2skYoWK873j+Hvk3vw78lSn+UzBQGfDd6e8YvWk0wb8GM33/dNZHrud08mlSc1NL3XeXRvUY6NeAb3ZHEn0zswLeXfUnR+xSjXck8QgRtyOY2XGmdjR+7AewagBNB/BrWBwbz1znzT5NCHSzY/dVfUdbDTXpB30+ga3TtDV3+sxCURSa2TejmX0znm35LLdzbrM/YT/74vexO243G678u8manakd7tbuuFm74W7jjru1O03tmtLAosFDf3t6f6APuy8l88GG8yyf0Eb+lvUf5UrsiqKMBD4EmgFthRBhughKknRpZfhK7Ezt6O/ZH1KuwJW/odu7RN/OY8aG83T2rsvzXb30HWb11u457Wd7aAHYe0HriYVetjW1ZZDXIAZ5DUJdoCYuPY7o1Ghi0mKITovmaupV9sbvZV3kunvXOJo7EuAQgL+DPwEOATS2bYyhSpuyHK1Neb13Y/73VzibzyYywK9Bpb7dqq68I/ZzwHDgex3EIkk6dzdhPN/yeUwMTODYj6AyhMCnWLA1EoAvHm+JSiVHfOWiKND3U7gdDZvegDqu4N2zyFMNVYZ42Hhob2L/R1peGtGp0Zy7eY5TSac4mXySrdFbATA3NMevnh8BDgH09+jP2PZu/BYWz8yN5wlqXBcrU6OKfIfVSrkSuxDiAiB/DZKqrNUXVmOkMuLxJo9DXhacWgXNBpOgsWH9yRM82d4NR+uaW3q3UhkYwoglsKQv/DIOBv0f+I0sVRPWxtb41fPDr54fY5qNAeB6xnVOJp3kZNK+YOHjAAAgAElEQVRJTiWf4vsz3/Pd6e/o7d6bZ3uN4pWVaXwZGsEHg3wq4l1VS3KOXaqxUnNT2XBlAwM8B1DXrC4cXw45qdD2GRbvjQLgmSBPPUdZw5haw5N/wG/jYe0kiNoF/eaAiWWZm2xg2YAGlg20U2loHzRbFb6KNZfWsC16G56+Aaw42Z7HAp1o3vDhm6TUJkpxT3ApirIDqF/ES9OFEH/+c85u4I1HzbErijIZmAzg6OgYuGbNmjIFnJGRgaVl2X9IKoqMq3QqI6696Xv57dZvvFX/LVyMnQk8/iqKKGCn31e8sSebtg0MmdTCpNLjKovqFpdSoMEtZg1uMb+RbdaQcJ83yLDS7T+iWZos9qbvZVf6brIKMjHK9cQvwo4tKzaTnJyMg4MDkyZNomfPoqeEqqPg4ODjQojWxZ4ohCj3F7AbaF3S8wMDA0VZ7dq1q8zXViQZV+lURlwjN4wUIzeM1H4Te0SIGdZCHP1BfLHtonCftlFE3EjTS1xlUW3jitojxNwmQsysK8Shb4UoKNB5DJl5meL1rfNFw6cbCcVYEcC9L3Nzc7Fq1Sqd96kvQJgoQY6V69ilGik8JZwLty4wrNEw7YGji8HEmoymj7HsYDS9fRzxdrDSb5C1gUcQPHcAvLrD1rfh59GQmaLTLsyNzPm89xTurEur9OJkVVW5EruiKMMURYkHOgCbFEXZppuwJKl81kasxVhlTH+P/pCRDOHroeVofjqZQlqOmhe6ees7xNrDwh5Gr9GumrnyN3zXWfuQWIHuSgIoikL27aLrx8TGxuqsn+qiXIldCLFOCOEshDARQjgKIfroKjBJKqscdQ6bozbTy72XtorjieWgySOv1QR+2HeVTt72tHSpo+8waxdF0ZZHfjoUjMzgp8dhng9sfUdbPlkH1RofVpzMtK4psWm1K7nLqRipxgmNCSU9P53h3sNBo9bWhfHoyu8xFiSl58rRuj419IfnD8CIpeDcWvsU8A/dYX4A/P0RJF0sc9OzZs3C7D/FyUzMTHB+3JlRG0fxd8zf5Y2+2pDLHaUaZ13kOpwtnWldvzVc3Axp8Wj6zOb7zVfwc7aho5e9vkOs3YzMwHe49iv7DlzcCGd/h/3zYN9ccPQFt47asg/WDQv/9xHLJu8WJ3vn3XeJi43D0LouU6bN4LWXhvH67td5ZfcrPOXzFC8HvqzdPasGk4ldqlFi02I5lniMqQFTtVUEj/0A1k5syfMnJuUs3z3ZSj5QV5WY1dFudBLwJGQkwfn12jr5Z37RPnPwXybW2gRvbAGGpmBoDAYm9/4bYmFKyDePEXkjjb137Dlz+xrxhw+yPOAN5l5Zy/Lw5Zy9eZb53edrp+lqKJnYpRplXeQ6VIqKwV6D4WYkRO1CBE/nm70xeNWzoLdPUY9kSFWCpQO0m6z9AsjLhPRESLum/Uq/BmnXIf065GeBOlf7NLHmNqjzQJOrPZafjVf2bbwR2gy3X/v1LhBgYc574jgT1gSzyKEHdV3aQwN/sPUAVc2ZmZaJXaox1AVq/oz8k85OnXG0cIT980FlyEGbAVy4Hs3nI/xkTZjqxNhCW1DMvvQF2vbs+ptubf3JvJPI7N8PkJJ0jRfb1qGfbT62yaeZmn6aCQkbWRz2I/U1GjCxgQZ+0KAluLTTLs8sx9Oy+lZz/omSar0DCQdIzk5meKPhkJ8Dp36CJv35v8NpNLQxZYi/k75DlCqLYgAW9lg4NWfacxNJcu7NkCON2WL7JO0fW8X3/Vdw09yGpxq3JK73TGgxQvtbwNHF8OtYmOMJP42CEyu12yhWMzKxSzXG2oi12JvaE+QcBBf+guxbXHYZwdHoW0zq4omxofxxr40sTQxZNqEN/i51eOnnk2w5e50AhwB+6PMjWULNU/F/cqXzi/DMTng3AcZv0pYdvhEOG6bA3EawpB8cXAC3qkfBfvmTLtUIN7Nvsid+D4O9BmtXPBxfCrYefHapPrbmRjzR1kXfIUp6ZGVqxLIJbWjpUofnV5/gqSVHuX3LkSV9liAQTNg6gfCUcDAwAvfO0O9TeOUMPLsPgt6E3DTYPh3m+8OiYO0S2tx0fb+th5KJXaoRNlzZgEZotCUEki9BzAFuNh3D35duMr6jB+bG8nZSbWdlasSKiW15o3djzl9L48kfjzB1RSIhrp9hYmjKpG2TOJV06t8LFEU77x78rnbt/dST0OsjyM+Gja/A3Caw4SWdPWClSzKxS9WeEIJ1Eeto5dBKu3lD2FJQGbEqpxOGKoUx7Yp+IlGqfSxMDJnSvRH73w5mzmN+qDUFfLw+mdQrk1EJKyaHTubI9SNFX2znCZ2mwguHtE/QNh+mXX//Q3dtmYQjiyD7duW+oYeQiV2q9k4knSA6LVo7Ws/PhtM/oWk6kJVnsujZzJF6VibFNyLVKqZGBjzexoVtrwSxdEIbvO2cSQifQFamDZO3P8/s3b8TdyvrbvXawhQFXNrC0G/g9UswYB6oDGDLm/BFU+0oPvlS5b+p+8jfT6Vqb23EWiyMLOjt1lv7gEtOKsfsh5CSmccoObcuPYJKpRDcxIHgJg6cS2jGt3sbsiftU1ZHf8SPBy9iT3tau9nR2t2W1m52NGtghaHBfeNhU2to87T269op7b2d02vgxApo3Bc6vgRunbT/GFQimdilai09L53QmFAGeA7A3Mhc+xfL3puF0Q1paJNBUKN6+g5RqiZ8nWz4ZnQQabmtmLztRc4rv+JsYMLJ2FZsOnsdAHNjA4KbOjCpswcBrraFG2joDw3/D7q/r33i+egiWDYAGgZoE3yzIdrtAytB9ZqKSTyH7a1TxZ8n1Rpbrm4hW52tLfh14zzEHeGOTwj7Im8ysrULBvKBJKmUrE0sWT5gEd2cu3JZs4KJA69ycFp35o8OYHgrJ/ZeTmbYwoOM+PYgW88loin4z3SNRV3oNg1ePQ8Dv4ScNPh9orbQ2eFvK2U1TfVK7Afn0/LMDFg+GBJO6DsaqQpYH7meRraN8K3rq71pamDCT7mdARjZ2lnP0UnVlYmBCfOC5zHAcwD/d+L/+OXKdwzya8DHQ1tw6J0efDDQh8S0HJ5bdZzuX+xm+cFosvLUhRsxMtOuh58SBk/8BDZOsHUaRO2p8Pir11TM4K+JyLKi0bV1sDhYe1e6+/tleuRYqv6iU6M5e/Msb7R+AyU/C878QoHPEFaeTqdLo3o425oX34gkPYSRyohPOn+CpZElS84tIT0vnentpmNpYsjEzh6M6+DG9vAbLN4XxYwN55kXepmQdq48G+SFjfl91SNVKmg6QPuVcFxbm6aCVa8Ru6EJCc6DYOopCHoLLm+Hb9rCxle1xYKkWmXL1S0oKPR17wvn1kJuGqcchnI9NYcn2sibplL5qRQV09tNZ1KLSfx2+Tfe2f8O+QX5ABgaqOjfogHrXujEH893oIOnPd/uucKYHw6TmpVfdINOgdoVNBUdd3kuVhTlc0VRLiqKckZRlHWKolTOtjSm1tB9Orx8CgInaO9Azw+Av2cWXepTqnGEEGy+upk29dtoC34dXwr1mrLoqiP2Fsb0bOao7xClGkJRFF5u9TKvtHqFLVe38GzosyRmFh5IBrrZ8d3YQJZNaEvEjQzGLjlCavZDknslKO+IPRTwFUL4AZeBd8ofUilYOsCAufDiUWjSD/Z9Ad+0h6jdlRqGVPku3LpAdFo0/Tz6wfUzkHCcDN+x7LiYxGOBzrIujKRzT7d4mlmdZ3Hu5jke2/AYoTGhD5zTtXE9vhvbigvX0xi35ChpOfpJ7uXd83S7EOLuHYPDgH7uVtl7wYglMGmnttTniiGw9V1thT+pRtoctRlDlSG93HppR+uGpvya3wl1geDx1nIaRqoYg70G89ug33C1cuW13a/xwYEPyMrPKnRO96aOLAwJ5HxCKuOXHCUjV/2Q1iqOLoc1E4EtOmyv9JwD4dm90OYZOPyN9gZr4jm9hiTpXoEoYEv0Fjo37IwNKjjzK6L5MFaeSqWtux3eDtW3jrZU9blZu7Gi/wqeafEMf175kxF/jeBM8plC5/TycWTBmABOx6cyYelRMis5uStFPjJ7/wmKsgMoatuZ6UKIP/85ZzrQGhguHtKgoiiTgckAjo6OgWvWrClTwBkZGVhaFv8X1y4ljKYXv8ZQnUGU5zjinQeBUnG/npc0rspWE+OKyIlg/o35jK87noGpN2lyeSHrvWbzynk3nmlhTCensu9nWRM/r4pU2+OKzIlkxc0VpGpS6WfTj942vbVbMv7j6HU1357OpYmdilcDTTExKN9zFcHBwceFEK2LPVEIUa4vYDxwCDAv6TWBgYGirHbt2lXykzOShfh5jBAzrIVYNlCIO3Fl7lencVWimhjXhwc/FG1WtRGZeZlCfNdFiIUdxSs/nxC+H2wVWblqvcVVkWRcpVOZcaXmpoq39rwlfJf5irGbx4rTSadFQUHBvdfXn4wXHtM2ijGLD4nsvPL9fAJhogQ5tryrYvoCbwGDhRBZxZ1fHqtXr8bJxZXg4O64u7uzevXq4i+yqAujVsHgr7WlNb/tqF0iKVVb+Zp8QmNCCXYJxjz5Elw/TXaLsWw+l8iQgIaYGVf8UjJJup+1sTWfBX3G7C6zibgdQcjmEAauG8i3p74lNi2WIf5OfD6iJQevpPDMijBy8jUVHlN55yYWAFZAqKIopxRF+U4HMT1g9erVTJ48mWvxcYAgJiaGyZMnlyy5Kwq0GgfP7wdbd1gzBi5srIgwpUpw8NpBUnNTGeA5AI4vB0Mz1hd0JFddwBNtZHleSX8Geg5k+4jtzOw4kwYWDfj29LcMWDeAkM0h5Fns44PBbuyPvMmui0kVHku5njwVQnjrKpBHmT59OllZhX8hyMrK4qU3XyK3ZS7mhuaYGZphZmiGRmhIy00jLU/7lZ6Xrv1zbhppDepha+yG646puF7fi1vjgbhaueJk5aTddUeq8jZf3YyNiQ0d7FrA2TGI5kNZcTIVXydrfJ1s9B2eVMtZGVsxrNEwhjUaRmJmIluubmFj1EY+OfIJhoohQV3a4tbQA2hQoXFUi5ICsbGxRR6/ff02Xx7/8qHXmRiYYGVshbWxNdbG1tibO3ALhTNKPhlxWyBOu4jHQDGggUUDPGw8eLrF0wQ6BlbI+5DKJys/i11xuxjoORCji39BXjpXXB7jwpE0Phrqq+/wJKmQ+hb1meA7gQm+E7h8+zKbojaxKWpToZurFaVaJHZXV1diYmIeOO7i6sqRMUfIUmeRrc4mW52NgWKgTeQm1pgYFL3BgshJ5/bPI4i9cYqYdpOIsXEkLj2Ok0knmbB1Ak/6PMnUgKmYGppW9FuTSmF33G6y1dnah5L+mgb1mrIk1hFTowSG+DfUd3iS9FCNbRvTOLAxL7d6GYWKrzhaLR7PmzVrFubmhQs6KYYmdH/yZcyNzKlrVhcXKxca2zbGq44X9czrPTSpAyimVtiFrMW/QVuG7PmGqaq6fN71czYM3cDjTR5nZfhKRv41ktPJpyv6rUmlsOXqFhzMHQgUJpAQRq7fWDacvk7/Fg2wNpVTaVLVp1JUKJWw6Ua1SOwhISEsWrQINzc3FEXBzc2NbpOmc0Tlw/XU7LI1amwBo38Bz26w/gU4sQJzI3Pea/8ei3otIleTy7gt4/jy+JfkafJ0+XakMriTc4f9Cfvp79Ef1cmVYGDMJiWIjFw1IXJPU0kqpFokdtAm9+joaHbu3El0dDTLZ71OgRB8svli2Rs1NofRa8C7h3afwmM/AtChYQfWDl7LUO+hLDm3hFEbR3E+5byO3olUFqGxoaiFmn4u3eHML+AzhOWn0mjsaEmr/+5kI0m1XLVJ7P/lYmfOs129+Ov0NY5EpZS9ISNTbRH8Rn1g02tw7g8ALI0t+V/H//FNj29Iy00jZFMIC08tLHpzW6nCbbm6BXdrd5pdvwg5qVx1HcHp+FTGtHWtlF9tJak6qbaJHeD5rl441TFjxobzqDUFZW/I0ET7IJNLO9jwMtyKuvdSkHMQa4espa9HX749/S0/XfxJB5FLpXEj8wZhiWH09+iPcmI52HmxOK4hJoYqhrWSuyRJ0n9V68RuZmzA9AHNuJiYzs9Hi14SWWKGxvDYD9rdTn6fCOp/59VtTGyY3Xk2Qc5BzAubR8TtiHJGLpXG1uitCAT9bJpC7CHyWo7lz1PXGOjXEBszedNUkv6rWid2gH6+9enoZc/c7Ze5lVnOm5x1XGHIQrh2EnZ8WOglRVGY2XEmlsaWvL3vbXI1ueXrSyqxzVc309y+Oe6XQkFlxEZVNzLzNIyRN00lqUjVPrErisKHg5uTkavmi+2Xyt9gs4HQdrK27O+lrYVesjez56NOHxFxO4Kvjn9V/r6kYkWnRhOeEk4/t15w+mdoOoBlpzNp4mhFK9fK2bBLkqqbap/YARo7WjGugxs/HY3lXIIOtsbr9RHUbwHrn4PUhEIvBTkHMbrpaFZdWMXBhIPl70t6pHv7muYpkH2Lq24jOBOfyph28qapJD1MjUjsAK/0bIyduTGzt1wof2NGpjBimXae/Y9JoClcJP+1wNfwsvFi+oHp3M65Xf7+pIfaGr2VQMdAHM+uhTpuLI53xdRIxdAAJ32HJklVVo1J7DZmRkwO8uRAZIpuRu11vWHgPIg9CHvnFHrJ1NCUz4I+IzU3lRkHZ8glkBXkyp0rRKVG0adeK4jeR27LJ/nz9HV501SSilFjEjvA6HauWJoYsmhvVPEnl0TLJ6DlGNgzB67uLfRSE7smvNLqFdb9sg4HZwe6dy9FnXipRLbHbEdBoUdyHCgGbFL1kDdNJakEalRitzY14ok2Lmw6e5342zra96P/52DvDX88AxnJhV5SnVJxffl1bl67iRClrBMvFSs0JpSAev7UO7sWGvdl6Zlsmta3IsBF3jSVpEepUYkdYGJnDxRg6YFo3TRoYgkjl0L2be3N1IJ/H4R6b/p7aHIL74aSlZXF9OnTddN3LXY19SoRtyPobdYQMpOJdh/J2QR501SSSqLGJfaGdcwY6NeANUdjSc3O102j9VtAn1kQueNeyQF4eJ34hx2XSm5HzA4AesSdB2snFl3zwNRIxRB/edNUkopT3j1PP1IU5cw/2+JtVxSlShTFfibIk8w8DT8d0WGCbf00NGgJO2ZAXiagrRNflIcdl0ouNCaUlrZNqR+1j1y/J/nzdCKD5E1TSSqR8o7YPxdC+Akh/IGNwAc6iKncmje0obN3XZYeuEqeuhw1ZO6nUkHfTyEtAQ5+DRRdJ97I1IhZs2bpps9aKi4tjgu3LtCrwAQUhc2G8qapJJVGuRK7ECLtvm8tgCqz7u+ZIE+S0nP581RC8SeXlFtHaD4M9n8FqQkP1Im3cbTBebwz/R/rr7s+a6HQ2FAAekUdA+9e/Hg2j6b1rfCXN00lqUTKPceuKMosRVHigBCqyIgdIKhRXZrWt2LxvijdrjPv+T8QBfdqydxfJ/7ExRNYtreUFSDLKTQ6FF8LZxqmJRLtPoJzCWmEyJumklRiSnFJT1GUHUD9Il6aLoT4877z3gFMhRAzHtLOZGAygKOjY+CaNWvKFHBGRgaWlpYlOvdAQj6Lz+bxWqAJfvV0t72rR9RK3GJ/50TAHNJsmhSKa3HSYiJzI5npNBMT1cO356sspfm8KtPD4rqlvsWMhBk8m23BM7cSedZqIQcTBV91M8fcqOITe3X7vPRNxlW5goODjwshWhd7ohBCJ1+AK3CuJOcGBgaKstq1a1eJz83N14i2s0LF6EWHytxfkXLShfi8kRCLugtRUFAortNJp4XvMl+x/Nxy3fZZRqX5vCrTw+Jadm6Z8F3mK2I/thd3Nr4vvN/dJKavO6P3uPRNxlU6VTWu8gLCRAlybHlXxTS679shQDn2qdM9Y0MVEzp5cPCKjsoM3GViCT1mQEIYnP2t0Et+9fxoW78ty8OXk6/R0XLLWiQ0JpRmxra45Kv5NrUjiqLwYrC3vsOSpGqlvHPsnyqKck5RlDNAb+BlHcSkU2P+KTOweJ+Oygzc1XI0NPCH0H+XP971tO/TJGUlsTFqo277rOESMxM5nXyaXndSyHYJYvG5Asa0daWBjZm+Q5OkaqW8q2IeE0L4Cu2Sx0FCCB0uQdGNu2UGNp7RYZkB+Hf5Y/q1e8sf7+rQsAPN7Jqx5NwSCoSOllvWAn/H/g1Ar5REfinogZGBwgvBXnqOSpKqnxr35GlRJnT2AHRYZuAutw73lj+a5PxbR0ZRFCa2mEh0WjS7Ynfpts8abHv0dhoppjgb1WF2lDtj27vhYGWq77AkqdqpFYndqY4Zg3RdZuCuXjNBFOAZtbLwYddeuFi58OO5H2VZ3xJIzkrmZNJJet1KYpdZLwyMTHiuqxytS1JZ1IrEDjCpSwWUGQDtPqkdX8IxaQ/EHbt32EBlwPjm4zl78yzHEo89ogEJtNMwAkHvzAxmXW/D+I7u2Fvqf7moJFVHtSax+zrZ0NHLnhWHosnX6Hjeu/Or5BrbwrZ34L7R+RDvIdib2vPjuR91218NFBq9HU8N5Bu24KaxM8908dR3SJJUbdWaxA4wsZMH11Nz2HouUbcNm1gS7T4G4o9B5N//HjYwYazPWA5eO0h4Srhu+6xBUrJTCLsRRq+0VBakdmZiZw9sLYz1HZYkVVu1KrF3b+qAu705P+6/qvO2E+sHg42Ldhu9+0btjzd5HEsjS5acW6LzPmuKnXE7KUDQKVvhkHF7nv7nZrckSWVTqxK7SqUwoZMHp+LucCJWt5tQC5URdHoZ4o4U2kbPytiKUU1GERoTSkxajE77rClCr2zCLV/N0cyOjO/SRJbmlaRyqlWJHWBEoDNWpoYsqYBROwFjwbI+7P280OEnfZ7EUDFk2flluu+zmruTc4ejSSfolZnJVqPe95amSpJUdrUusVuYGDK6rStbziWScCdbt40bmWpH7dH7IObgvcN1zeoywHMAm6I2kZWvw4ekaoBdsX+jQdAg3ZE+3YKwNNFdsTZJqq1qXWIHGNfBDSEEKw5F677xwPFgUQ/2zCl0eHij4WSrs9kWvU33fVZj2y7+ilO+mjMFfRjXwU3f4UhSjVArE7uzrTl9fevz85FYsvLUum3c2Bw6TIGoXRAfdu9wy3otcbd2Z33ket32V42l5qZy+FY4QZn5eAWNwdxYjtYlSRdqZWIHeLqzB2k5av44Hq/7xts8DWa2hUbtiqIwrNEwTiSdIDo1Wvd9VkPrz/yORgGTHD9Gd2qi73AkqcaotYm9lastLZ1tWHogmoICHT/yb2IF7V+EiG1w7dS9w4M8B2GgGMhRO5CZL1h7YjkN1GqGDfkAUyMDfYckSTVGrU3siqIwsbMHUTcz2XM5ufgLSqvdZDCxKbRCpp55Pbo4dWHDlQ2oC3Q8BVSN5Ko1zD+RQqzJLYKwxbOpv75DkqQapdYmdoD+LRrgaG1SIQ8sYWoD7Z6Fixvhxvl7h4c2GkpydjIHrx18xMU1V0GB4I3fzmAhNqNWFAb6PanvkCSpxqnVid3IQMW4Du7sj7zJpcR03XfQ/nkwtoS9c+8dCnIOws7UjnUR63TfXzXwxLS5fPtcH9a9s5TI1y9x7qKFvkOSpBqnVid2gDFtXTE1UrH0QAWM2s3toM0kOL8Oki8DYKQyYpDnIHbH7eZWzi3d91mFPffhV/z+5fto0rRTXzkp+Tz73AusXr1az5FJUs2ik8SuKMrriqIIRVHq6qK9ymRrYczwVs6sPZlASkau7jvoMAUMTWHfF/cODfUeilqo2XilcrbO0xQIYlOyiErOIDIpncs30rlwPY3z11I5G5/K6bg73EjLqdAYtp9P5IcvP0GoC3/GWVlZTJ8+vUL7lqTaptwLhxVFcUG736mOC51Xnomd3PnpSCw/HYnlpR6Nir+gNCzrQeuJcOQ76PY22HnibeuNX10/1kWuY6zPWBRF0W2f/xBCcCxRzcwv9xCVnFns+X7ONvRs5kjPZo40a2Cls7hOxt5m6pqTaNJuFvl6bGy1/dGRpCpJF0+EfAm8Bfypg7b0wtvBiq6N67HicAzPdvXC2FDHM1SdpsKxH2DfPBiyANDeRJ15aCbnU87jW9dXt/0B+yNuMmfbRc7E59LIwYiPhvpiZWKISqWgUsBAUf75s/b7SzfS2RF+gy93XGZe6GWc6pjRs5kDvXzq09bDrkyfSVJ6DqHhN/hi+2UcrEzJdHEhPu7BJO7q6qqLtyxJ0j/KldgVRRkCJAghTlfUqLOyTOzswVNLjrL2RDxPtNVxorGqD63GwvHlEPwuWDekr3tf5hydw7qIdTpN7Kfj7jBn20UORKbgVMeMSS2MeWd0EAaqR///6dHMkRe6eZOUnsOui0mEhifxS1gcyw/FYGViSBsPO/xd6tDSpQ7+znWwMS+6AmPcrSy2nU9k2/lEwmJuIwR4O1iyeFxrDjp9wtMTx5Gb9+9GJ+bm5syaNUtn71+SJFCK249TUZQdQP0iXpoOvAv0FkKkKooSDbQWQhT5+7aiKJOByQCOjo6Ba9asKVPAGRkZWFpalunaRxFC8PHhHG7lCD4NMsPEoHT/UBUXl2n2DdodeY5458Fc8Z4AwIqbKzibdZZZzrMwVpVvY4lrGQWsjcgj7IYGKyMY5GVMsKshuVmZZf68cjWC8BQNp5I0RNzRcD1DcPenxdFcwbOOCk8bA5wsVUTe0XD8hoaYNG3SdrVSEehoQKCjIU6WCoqiYJUWwftr32XnxlTybqlxcHBg0qRJ9OzZs1zvXZcq6uervGRcpVNV4yqv4ODg40KI1sWdV+yIXQhR5N86RVFaAB7A3dG6M3BCUZS2QogHtigSQiwCFgG0bt1adOvWrbiui7R7927Kem1xzN1SGLXoMFEGrjzfrXQbKZcorqxQXC5twWXMl2Buh0WiBRO3TSTPPY/enr3LHPfKQ9HM2HYeMyMDXu3ZmKe7eNyrkljez6vPfZKGWrIAABUDSURBVH9Oy8nnXHwqJ+PucDruDqfi7nDo2r83QwPdbHmyc336NK+Pq735A21lr/uZm10deW/S83RWd6+w/4/lUZE/X+Uh4yqdqhpXZSnzVIwQ4izgcPf74kbs1UE7T3t6NHVg4e5IRrd1oY65jrdn6/QKnP0Njv0IXd8k0DEQZ0tn1kesZ6DnwDI1ufj/27vz6CjKdI/j3zcrIYEESEgiEEgIEHYhbBIUguDIooCKjKIsDgIXB+8oOKIZz8w4Fy9H8SKM43gHh2GEK7igEJBVFjcWZYnIEoEEhCCrQEjI3v3cP7rFyJKQpDvVHZ7POXW6qe5669eV5unqt6ve+iyT6Sv3c2dCQ15+oINbLwBdt5Y/PePD6RnvOPhJRDh5sYD0kzm0ja5Lw7q1rr9w3jm+yFxFfngo/ZsPIv87Fw+ZrJS67KY/jv1Kv787gdzCEt7YlOH6xqPaQYu7HEfIFOXhY3wYGj+UbSe3kZVT8cHI/rr+INNX7mdQ+2jefDTRrUX9WowxRIcGkdyqYdlFHSDtHdbV8qOefx0SIxOrJ6BSNymXFXYRaebNe+s/aRVVh/s7N2b+5iOuvxAHQK+nIO8spDlOyhkSPwSDYVnGjR9UJCLMXPMdr647wLBOjZj961vx9/Xgz2i7nYLtb/FpcAh3xv4KPx8dnlcpd/LgamCdp/q3BGDWugOubzzmNmjSHb6cA7ZiooKj6HlLT5YeWorNbit3cRHhpZX7eX3jIX7dtQkzh3fEz5OLOkDmRr4sOEmeEfo37W91GqVqPA+vCNZoFBbEmJ7NWLIzi/STF13buDGOvfbso46hBnCciXry0km2ndxW5qJ2u/DH1L3M/fwwo25rykvD2pd7GKNH+PqfrKtbj9CAULpGdbU6jVI1nhb265jUpzkhgX68svo71zfe4lcQ0Rq+mAUiJMckUzegLksPXn+cdptdeO7Db3l7y/eMvyOOP9/bFh9vKOrZWRQdWMWm2rW4s+md+Ptc+/h3pZTraGG/jrDaAUzqE8/69NNsy/zRtY37+ECv38HpfXBwLYG+gQyMHciGYxu4WHT1NwSbXZj6/je8u/0Yk/vG89yABLcNQ+By295kc1AQl6REu2GUqiZa2MswNqkZUXVrMWN1OuWdyFVh7e6H0CaOvXYc3TGFtsKrLnZttwvTluzmo13HmdK/JVPuauU9RT3/Amyfz7pGLakTUIfuUd2tTqTUTUELexlq+fvydP+W7Dp6gTV7T7m2cV9/6DkZjm6Bo1tp06AN8WHxv7hsnojw4op9vL8jiyfvbOH6Acrcbfs/KS7KYSP59G3SF39f7YZRqjpoYS/HfZ0b0aJhCC+vSafEZi9/gYro9CjUbgBfvIYxhiHNh7D7zG4yszMBx1E58zcf4bGkWJ7q52VFvbgAtr7Jlrju5JTkaTeMUtVIC3s5/Hx9+P3dCWSeucT7Oyp+ElGZAmpD94lwYBWc2seguEH4Gl9SD6Xyj88ymLPhECO6NOGFwa29p/vlJ98sgkunWd4gmtDAUHre0tPqRErdNLSw34B+rRvSpWk9Zq07wPlLRa5tvOs48A+GL2cTUTuCpEZJvJu+lJdW7mNQh2heuq+99xV1uw02z+HiLR3ZcG4PA5oN0G4YpaqRFvYbYIzhhcFtuJBfzKh5X3GxoNh1jdeuD13GOsaQuXCURr53kFvyI51bnWHWg7d6x3HqV9q/HM5lsq7lHRTZixgSP8TqRErdVLSw36COTcJ485HOpJ+8yNh/fc2lwhLXNd5jEhgfji1/iX+uC8JXgomN3e/6C35UBxH48jWoH0dq/jFiQ2Np26Ct1amUuql4YeWwTt+ESOb8uhNpxy4w7t/bKSgufwiAGxLaiBPNhxN56D36RpYwrOVgNl3nmHaPd/gz+GEXx7o8ys7TO7m3+b3e15WklJfTwl5BA9pH8+rwjmw9/CMTFuygsKRqxd1mF978NIMH9/UC48PrjdbyQKthFNmLWH14tYtSV6MvZ0NwBCtq+WMwDIodZHUipW46WtgrYWinRsy4rz2fHjjDk4t2UVzJwyCPnL3EiP/dwoxV6bRNSMCe+BsC97xLG7s/8WHxLDvkZZeRPbEbMtYj3SaQengl3aK6ER0SbXUqpW46WtgraUTXGP50TxvW7D3FlPe+wV6BM1NFhAVbv2fA7M/57lQOs0Z05O+PdKZW8lTwC8J8OoOh8UPZfXY3mRcy3fgqXOzL2RAQQlpcd7Jys7in+T1WJ1LqpqSFvQrGJMUybUACqd/8wL/2FGG3l1/cT2TnM2reV7ywdA9dmtVj7VN3MKxTY0c/dEgE9JgIe5YwKKQ5vsa3QuO0W+r8Edj7ISSOITVrI0F+QXpSklIW0SseVNHE3s3JL7Ixe/1BJi7cQeem9YgICSSiTiDhztv6wQH4GPho13H+mLqXEpvwl6HteKR7zNU/LPacDF+9RfjmN7i90e2syFjBk52exNfH15oXeKM2vw7Gl8Ku41izeiT9YvpR2//q654qpdyvSoXdGPMn4HHgjHPW8yKysqqhvM3v+rXgyJEjrD14lrX7rh5TxsdAaJA/5/OKSWxaj1eHd6RZePC1GwuqB0mTYcN/MaTVy2zK2sSWE1vo1aiXm19FFVw6C7sWQocH2XTxADnFOdoNo5SFXLHHPktEZrqgHa9ljGFYiwBmP96HS4UlnMkp5ExuIWedt2dyHFNCVB0eva1Z+ScddZ8IW/9O729XEhYYxrJDyzy7sH/1DyjJh6T/ZPnuv9KwdkO6RXWzOpVSNy3tinGx4EA/ggP9rr9HfiMC60Cvp/Ffm8LAXqP44OgGsguzCQ0MdV1QVynMcRT2VgM5G9KAL45/wei2oz2/60ipGsxUZZxxZ1fMGOAisB2YIiLnr/Pc8cB4gMjIyMTFixdXap25ubmEhIRUall3cnUuH1sh3bdNJC0kgnGh+YyoP4JedSq+1+7u7RWXMZ+YYx+xo/NMUsniw/Mf8nz080QHlH2Y483yd3QVzVUxnpqrqpKTk3eISJdynygiZU7AJ8Cea0xDgEjAF8fRNdOBeeW1JyIkJiZKZW3cuLHSy7qTW3J9NVfsf6wrw97rLw+teKhSTbh1e51OF/lzfZGlT4iIyPDU4TJi+Qjrc1WB5qoYzVW9gO1yAzW23MMdRaSfiLS7xrRMRE6JiE1E7MBcQDtWXanTKExYDEMuZvPt2W/JuJBhdaKficCq30NAMPT7EwfOH2D/uf36o6lSHqBKx7EbY0p/3x6GY09euYpfAPSexuDjB/A3vixKX2R1op/tT4XMTZD8BwgOZ0XGCvyMHwNiB1idTKmbXlVPUHrZGPOtMWY3kAw85YJMqrQOI2hQrzmDiwzLDi3jfME1f8KoXkWXYPXzENkOujyGzW5jReYKejXuRf1a9a1Op9RNr0qFXUQeFZH2ItJBRO4VkROuCqacfP0g+TlGnTpGga2Ad7971+pE8Pn/wMUsGDgTfP3YemIrZ/LPcG/ze61OppRChxTwDm2GEd/4Nm7PL2LRvoUU2gqty/JjBmyeAx1GQNPbAEjNSKVOQB16N+5tXS6l1GVa2L2Bjw8M/Ttj8oo4V5TNioMWjR8jAqungW8g9H8RgJyiHDYc3cCAZgMI8A2wJpdS6he0sHuL0MZ07fcyrQuL+PfO2dilckMFV8mB1XBwLfSZBnWiAHhn/zsU2Ap4oOUD1Z9HKXVNWti9iOkwnNF1W3O4+CKff/Ov6l15cQGsehbCW0H3CQBcKr7Egv0L6N24N60btK7ePEqp69LC7mXuuuctomzC/B2zHUenVJfNc+DC9zDwZfD1B2Bx+mKyC7OZ0GFC9eVQSpVLC7uX8Q8O55EW97PdT9i78snqWen57+HzV6HNUIjrA0BecR5v73ubpFuSaB/RvnpyKKVuiBZ2L3R/92cIMX78+/hGOLDGvSuz2x1dMMYHfjX98uz3D7zPuYJzTOw40b3rV0pVmBZ2LxQSEMIDCQ+xNjiY48snO8ZDdwe7HZY/CQdWQd8/QGhjAApKCpi/dz7do7pza8Nb3bNupVSlaWH3UiPbjsL4+LIwoARSn3QciuhKdjssnwy7FkDvZ6HHpMsPLTm4hLP5Z5nQUfvWlfJEWti9VFRwFHfHDmBJaCjZB1c6CrCr2O2QOtlxVaTe0yD5eXBewq/IVsS8PfNIjEyka1RX161TKeUyWti92Oi2o8mXEj5o2gFWTYO9S6veqN0Gqb+FtJ+K+nO/eHjpoaWczjutR8Io5cG0sHuxhPoJ9IjuwTu1fCiOaAHvj4Ylj0N+JQcKs9tg2ROQ9n/Q57mrinqxrZi3vn2LjhEd6RHdwwWvQCnlDlrYvdzotqM5XfAjq5J/5yjGe5bAGz0hY0PFGrLbYOkk+GYRJKc4zi69wvLM5Zy4dIIJHSZgTDnXbVVKWUYLu5dLuiWJ+LB45u9fgPR+FsZ9AoEhsGAYfDwVH1tB+Y38VNR3L3aMr97791c9pcRewtzdc2nboK1nX1hbKaUXs/Z2xhjGthtLyhcpLEpfxMOtH4YJn8H6F2HrG3QJ+hhaNoAmzh86ReDicTiTDqfT4cx++CENTu1xHNJ4xzPXXM/KwyvJys3i2W7P6t66Uh5OC3sNMDhuMGuPrOWV7a/QIaID7cLbwd3/DS3vxufdx2DeXZAwCHJOwpnvoPDizwvXDoeGrWHwLOjy2DXbt9ltzN09l4T6CTo0r1JeQAt7DeBjfJjeazrDlw9nyqYpvHfPe4QGhkJcb77uOofbc1dCxnpoEO8YRz2ilaOYRyRAcHi57a85soYjF48wq88s3VtXygtUubAbYyYDTwA24GMRubqDVrldaGAor/Z+lVGrR5HyRQpz+s7Bx/hg8wuGoX+rdLvZhdm88c0bxIfF0zemrwsTK6XcpaoXs04GhgAdRaQtMNMlqVSltI9oz9QuU/k061Pm751f5fZyinKYuG4iP+T+wPPdn8fH6G/tSnmDqu6x/wcwQ0QKAUTkdNUjqap4OOFhdp7ayZydc+gY0bHS7eQV5zHpk0mkn0vnteTX9CxTpbyIkSqMMWKMSQOWAXcDBcBUEfn6Os8dD4wHiIyMTFy8eHGl1pmbm0tISEjlAruRJ+XKt+cz5YMpZL6XSeG5QiIbRjJu3Dj69et3Q8sX2Yt48/SbZBRmMDZ8LLcGu36gL0/aXqVprorRXNUrOTl5h4h0KfeJIlLmBHwC7LnGNMR5+1fAAN2Awzg/LMqaEhMTpbI2btxY6WXdyZNyLVy4UGoF1RLg8lS7dm1ZuHBhucsWlBTI42sel/bz28vHGR+7LaMnba/SNFfFaK7qBWyXcuqriJTfxy4i/USk3TWmZUAW8KFznV8BdqD8wyyUW6WkpFCQ/8sTk/Ly8khJSSlzuWJbMU9vepotJ7bwYtKLDIwb6M6YSik3qeqvYUuBZABjTEsgAHDT4ODqRh09erRC8wGK7cU889kzfJb1GS/0eIGh8UPdFU8p5WZV/fF0HjDPGLMHKAJGO78uKAvFxMTw/fffXzU/sEEgE9ZNIC40jtjQWGJDY4kLjSMsMIyUz1NYf3Q9z3Z9lgdbPWhBaqWUq1SpsItIEfCIi7IoF5k+fTrjx48nLy/v8ryAWgH0m9iP8wXnWXJ6Cfkl+ZcfC/ILIr8kn6cSn+KRNvrnVMrb6ZmnNdDIkSMBR1/70aNHiYmJYfr06Zfn28XOqUunyMzO5HD2YQ5nH6ZteFvua3GflbGVUi6ihb2GGjlyJCNHjmTTpk306dPnF4/5GB+iQ6KJDokmqVGSNQGVUm6jpxIqpVQNo4VdKaVqGC3sSilVw2hhV0qpGkYLu1JK1TBa2JVSqobRwq6UUjWMFnallKphqjQee6VXaswZ4OrBTG5MOJ450JjmqhjNVTGaq2I8NVdVNRWRiPKeZElhrwpjzHa5kYHmq5nmqhjNVTGaq2I8NVd10a4YpZSqYbSwK6VUDeONhf0fVge4Ds1VMZqrYjRXxXhqrmrhdX3sSimlyuaNe+xKKaXK4NGF3RgTZoyZVOrfo40xB53TaA/KtdoYc8EYs8KqTFfmMsbcaozZYozZa4zZbYwZ4SG5mhpjdhpj0pzZJlqV68pspebVNcZkGWNe94RMxhibc3ulGWNSrchURrYYY8xaY8x+Y8w+Y0wzq3MZY5JLba80Y0yBMebmuoiviHjsBDQD9jjv1wcynbf1nPfrWZ3L+e87gXuAFR60vVoCLZz3bwFOAGEekCsACHTeDwGOALd4wjYrNW828A7wuidkAnKtfF+Vk20T0L/U37O2J+QqNb8+cM6qXFZNHr3HDswAmhtj0oBtwDoROSci54F1wN1W5zLGvCIi64Eci7KUVnp7PS4iBwFE5AfgNFDuiQ3VkGu6iBQ65wdi/bfGX/wtjTGJQCSw1lMyWZjjWkpnex/wE5F1ACKSKyJ5ZS9eLblKb7MHgFUW5rKG1Z8sN/opDEwF/lDqsReAqVbnKjWvDx60x37F/G7AfsDHE3IBTYDdQB7whKdsMxwfMpuAxsAYPGePvQTYDmwFhnrQ9hoKrAA+BHYBrwC+Vue6Yv4GYLCV28yKyeq9JeVmxphoYAEwVkTsVucBEJFjItIBiAdGG2Mirc7kNAlYKSJZVge5QlNxnEX5MPCaMaa51YGc/IDbcex0dQXicHwgegTne789sMbqLNXNmwr7cRx7ej9p7JynrsMYUxf4GEgRka1W57mSOLqI9uAoDp7gNuC3xpgjwExglDFmhrWRQESOO28zcXyj6GRpoJ9lAWkikikiJcBSoLPFmUp7EPhIRIqtDlLdPL2w5wB1nPfXAHcZY+oZY+oBd2HdJ3HpXJ7kci5jTADwEfC2iHxgaapf5mpsjAly3q8H9AK+84RsIjJSRGJEpBmOvdC3RWSalZmc7/dA5/1wIAnYZ0Gmq7IBXwNhxpiffrvpi3XZrvV/8iFgkQVZLOdndYCyiMiPxpgvjTF7gFXAX3C8mQBeFJFzHpKrB5AAhBhjsoDfiEi1f+hckSsYx7eaBsaYMc6njBGRNItzFQO+xhgBDDBTRL6t7kzXybZKRJ6xKst1Mp0Aoowxdhw7YjNExLLCfo33/lRgvTHGADuAuR6S6284vuF/akUeq+mZp0opVcN4eleMUkqpCtLCrpRSNYwWdqWUqmG0sCulVA2jhV0ppWoYLexKKVXDaGFXSqkaRgu7UkrVMP8PFj3VzYsHrD4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "m1, m2, m3 = 4, 6, 8\n",
    "c, d = 0, 1\n",
    "xt = np.linspace(c, d, 8, False)\n",
    "x = np.array([-2.2, -2.8, -6.1, -3.9, 0, 1.1, -0.6, -1.1])\n",
    "y = np.fft.fft(x) / np.sqrt(x.size)\n",
    "a = np.array(y.real)\n",
    "b = np.array(y.imag)\n",
    "t = np.linspace(c, d, 40)\n",
    "data1 = [least_squares_even_trigonmetric(m1, i, a, b, c, d) for i in t]\n",
    "data2 = [least_squares_even_trigonmetric(m2, i, a, b, c, d) for i in t]\n",
    "data3 = [least_squares_even_trigonmetric(m3, i, a, b, c, d) for i in t]\n",
    "plt.xticks(xt, ['t0', 't1', 't2', 't3', 't4', 't5', 't6', 't7'])\n",
    "plt.plot(t, data1, label='order 4')\n",
    "plt.plot(t, data2, label='order 6')\n",
    "plt.plot(t, data3, label='order 8')\n",
    "plt.plot(xt, x, 'o', color='k')\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}