{
 "metadata": {
  "name": "",
  "signature": "sha256:c06499ba1fa9759a13111e7fff4f82eab9d43a509c8d4950f93349c03c046203"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Benchmarks\n",
      "\n",
      "This script shows benchmarks of the Non-Uniform Fast Fourier Transform (NUFFT). There is a Fortran implementation of the NUFFT (python wrappers at http://github.com/dfm/python-nufft/) and the pure-Python implementation of NUFFT (http://github.com/jakevdp/nufftpy/).\n",
      "Both are $O[N\\log N]$ for $N$ observations and $N$ frequencies, but the fortran version is about two times faster than the pure Python version."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "import seaborn; seaborn.set()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import nufft\n",
      "help(nufft.nufft1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Help on function nufft1 in module nufft.nufft:\n",
        "\n",
        "nufft1(x, y, ms, df=1.0, eps=1e-15, iflag=1, direct=False)\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import nufftpy\n",
      "help(nufftpy.nufft1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Help on function nufft1 in module nufftpy.nufft:\n",
        "\n",
        "nufft1(x, c, M, df=1.0, eps=1e-15, iflag=1, direct=False, fast_gridding=True, use_numba=True)\n",
        "    Fast Non-Uniform Fourier Transform (Type 1: uniform frequency grid)\n",
        "    \n",
        "    Compute the non-uniform FFT of one-dimensional points x with complex\n",
        "    values c. Result is computed at frequencies (df * m)\n",
        "    for integer m in the range -M/2 < m < M/2.\n",
        "    \n",
        "    Parameters\n",
        "    ----------\n",
        "    x, c : array_like\n",
        "        real locations x and complex values c of the points to be transformed.\n",
        "    M, df : int & float\n",
        "        Parameters specifying the desired frequency grid. Transform will be\n",
        "        computed at frequencies df * (-(M//2) + arange(M))\n",
        "    eps : float\n",
        "        The desired approximate error for the FFT result. Must be in range\n",
        "        1E-33 < eps < 1E-1, though be aware that the errors are only well\n",
        "        calibrated near the range 1E-12 ~ 1E-6. eps is not referenced if\n",
        "        direct = True.\n",
        "    iflag : float\n",
        "        if iflag<0, compute the transform with a negative exponent.\n",
        "        if iflag>=0, compute the transform with a positive exponent.\n",
        "    direct : bool (default = False)\n",
        "        If True, then use the slower (but more straightforward)\n",
        "        direct Fourier transform to compute the result.\n",
        "    fast_gridding : bool (default = True)\n",
        "        If True, use the fast Gaussian grid algorithm of Greengard & Lee (2004)\n",
        "        Otherwise, use a more naive gridding approach\n",
        "    use_numba : bool (default = True)\n",
        "        If True, use numba to compute the result. If False or if numba is not\n",
        "        installed, default to the numpy version, which is ~5x slower.\n",
        "    \n",
        "    Returns\n",
        "    -------\n",
        "    Fk : ndarray\n",
        "        The complex discrete Fourier transform\n",
        "    \n",
        "    See Also\n",
        "    --------\n",
        "    nufftfreqs : compute the frequencies of the nufft results\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M = 100000\n",
      "x = 100 * np.random.random(M)\n",
      "c = np.exp(1j * x)\n",
      "\n",
      "kwds = dict(eps=1E-8, iflag=-1, direct=False)\n",
      "\n",
      "k1 = nufft.nufft1freqs(M)\n",
      "F1 = nufft.nufft1(x, c, M, **kwds)\n",
      "\n",
      "k2 = nufftpy.nufftfreqs(M)\n",
      "F2 = nufftpy.nufft1(x, c, M, **kwds)\n",
      "\n",
      "print(np.allclose(k1, k2))\n",
      "print(np.allclose(F1, F2, atol=1E-8))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "True\n",
        "True\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Mrange = (2 ** np.arange(3, 21)).astype(int)\n",
      "kwds = dict(eps=1E-8, iflag=-1, direct=False)\n",
      "\n",
      "nufft_times = []\n",
      "nufftpy_times = []\n",
      "\n",
      "for M in Mrange:\n",
      "    x = 100 * np.random.random(M)\n",
      "    c = np.exp(1j * x)\n",
      "    \n",
      "    t1 = %timeit -oq nufft.nufft1(x, c, M, **kwds)\n",
      "    t2 = %timeit -oq nufftpy.nufft1(x, c, M, **kwds)\n",
      "    \n",
      "    nufft_times.append(t1.best)\n",
      "    nufftpy_times.append(t2.best)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.loglog(Mrange, nufftpy_times, label='nufft python')\n",
      "plt.loglog(Mrange, nufft_times, label='nufft fortran')\n",
      "plt.legend(loc='upper left')\n",
      "plt.xlabel('Number of Elements')\n",
      "plt.ylabel('Execution Time (s)');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAFsCAYAAAAtwdttAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VVW+//F3eiEhBBJ6CSVsIJAE7IgFHAtKk24DUeoo\njjOWe39T79w79965M3YQFQRRUUAQEOwNRRgbkEaATQ+EFBLSe845+/dHgoMOwiE5JeXzeh6enJK9\n1ncnJJ/svddey8eyLERERKR18PV2ASIiIuI5Cn4REZFWRMEvIiLSiij4RUREWhEFv4iISCui4BcR\nEWlF/L1dwLkYhjEcmFv/9FemaRZ7sx4REZGWoqke8c+hLviXA9O8XIuIiEiL0VSD3880zRogG+ji\n7WJERERaCo+f6jcM4wrgr6ZpjjQMwxdYAsQD1cBs0zQPAxWGYQQCXYEcT9coIiLSUnn0iN8wjMeB\nZUBQ/UsTgEDTNIcD/w48Wf/6UuAl6k75v+7JGkVERFoyTx/xHwIm8s8wHwF8CGCa5reGYVxa/3g3\nMMvDtYmIiLR4Hj3iN01zA2A766VwoOSs5/b60/8iIiLiBt6+na+EuvA/w9c0TcfFNGBZluXj4+Pa\nqkRERJq2Bgeft4N/BzAWWGcYxpVA6sU24OPjQ15eqcsLay6io8O1/9p/b5fhFa1530H7r/0Pv/An\n/QxvBb9V/3EjcKNhGDvqn+u6voiIiBt5PPhN0zwGDK9/bAELPF2DiIhIa6WBdCIiIq2Igl9ERKQV\nUfCLiIi0Igp+ERGRVkTBLyIizc6O1CyWbknHZr+oqV8EBX+TsGTJc8yceQc7d37HwoXzWLDgfkpK\nivnkkw8b3GZNTQ3vvrsJgOXLX2LTprddVa6IiFd9v/8Uf3vte9IOn1bwN4CCvwn44ovPePHFFXTv\n3pOKigpeeGE5hw4dZPv2bQ1u8/TpfLZseQeom+RIRKQlSD6Uz9LN6QQF+vPrqYkEB3p7Hrrmp8V/\nxd76/BDf7z/l0jYvG9CRqaP6/ez777+/ha+/3kF1dTVZWZncdddMRo8ew4MPzuXxx39Hz5692LRp\nPQUFBfj6+pKfn89jj/0KX19fMjOP8/e//w9ZWSc5dOggW7ZsYuzYCT+0PX/+ffTs2YsTJ47Trl0k\nixc/y5///HtuuukWrrpqBMeOHWXJkmfp0CGaY8eOsHLlywBs3/4lW7d+RklJEbNnL+Dqq6/h448/\nYN261QQEBNK9ew8ef/x3fPzxB+esXUTE29KPFrBkYxp+fj78afaVdAwP9HZJzZKO+N2kvLycv/3t\naf7616dYtWol8NMjbx98fHy4997ZdOjQgaeffp7f/vZPxMT04bHHfsuMGfcxbNilPwp9gIKC00yb\ndhcvvLCcbt26s2bNGsaOncAHH7wHwHvvbWbMmAnMnHkfMTF9uPfe2ViWRXR0J559dgkPPfQImzat\np6SkmBUrlvLccy+xZMnLhIeH8847G/Dx8Tln7SIi3mQeL2TR26mADwsnxRPXp4O3S2q2WvwR/9RR\n/c57dO4OPj4+xMb2ByA6uiM1NTX/8jmWZZ33tXO9DxAZ2Z6+fev2Jz4+kdTUndx22ySeeebvFBUV\n8f333zJ//oOcOpX7o3oMYwAA7dt3oKqqiqysk/Tu3YeQkBAAEhKG8d133xAXN/iCtYuIeNLhrGKe\nWZ+K3WHx4MQhxMW093ZJzZqO+N3kXNfVAwODyM/PA+DAgf3n3d7X1/ec4V9cXER2dhYAaWkp9O9f\nF9I333wrTz/9Ny6//Er8/Pzw8fHF4agb9HKudrp06crRo0epqqoCIClpFz179vrZ2kVEvOF4bilP\nr02hptbOvHFxJPSL8nZJzZ6C301+elofYPLkaTz11P/xm98s/CGUz37/7O26devOkSOHWLduzY/a\n9fPz48UXF7Ngwf31p/2nAXDrrWPZtm0rY8aMByAyMhKbrZYXXliEj4/Pj+rx8fEhIqId998/l4UL\n5zFv3ixKSoqZMGHSz9YuIuJpJ/PLeWJNMpXVNmbfNohLB3T0dkktgs/PnVJuRqzWtDTjjBnTeO21\ntT88P7M0ZX5+Pn/5yx955pklXqzO87Q0Z+vd/9a879Dy9z+3oIK/vrGb4vIaZt5icF1itx+939L3\n/0Kio8MbfFSmI/5m5lyn4b/88nMeeeRBZs+e74WKRERcK7+okr+vSaK4vIY7fhH7L6EvjdPiB/e1\nNK++uuZfXrvuulFcd90oL1QjIuJahaXV/H1NEgUl1Uy+vi83XtrD2yW1ODriFxGRJqG4vIa/r04i\nr6iKcVfHcOuVvbxdUouk4BcREa8rq6zlyTVJ5BRUcMsVPRk/ore3S2qxFPwiIuJVFVU2nlqbTGZe\nOaOGdWPK9X0veFvxrqw0Vu1bh91h91CVLYeCX0REvKaqxsYz61I4llPKiPgu3Hlj/wuG/p78fTyx\n/UWS89KoddR6qNKWQ8HfBFzM6nxffrmV6dMn8vbba8/R0r/68sut5Ofnu7pkEZFGq6m1s+jtNA6d\nLOaKQZ2495YB+F4g9A8WHuHlPa/j5+vH/PhZBPsHe6jalkPB3wRczOp8O3ZsY+HCXzNp0jSn2l6/\nfg0VFWWuLllEpFFqbQ6e37iHfRmFDOsfzf23DcTX9/yhf7wkkxdTX8FhWTx69Tz6tdM4gIZo8bfz\nbTj0Lkmn0lza5tCOQ5jY7+dXrHPX6nzbt3/Jt9/+A9PcT0REO7KyMtm48S18fPx+tLree+9txrIs\n7rnnXg4ePMBf/vIf/OEP/8nvfvcYERHtuOqqqxk4MI6VK1/G4XBQWVnJn/70F/z9/fmP//gdnTp1\n5uTJTAYOjOPRR//dpV87ERGb3cFLm9NJO3KaIX06MG9cHP5+5z8OzS7PZXHKy1Tba7hv8F0kdolr\n1RP4NEaLD35vKS8v56mnFpGZeYJ/+7dfM3r0mJ9dne+99zbz9NPPc/p0Pn/602957LHfkpS0i02b\n3v7R6nwjRlzHtm1f8Itf3EyPHj34y1/+yJYtmykvt7No0VO8884GQkNDadu2Lf/7v08CEBvbn8ce\n+y3+/v4UFBSwYsUb+Pv7s3Hjev7wh/8iKiqK119/ha1bP+Wmm0aTmXmcZ55ZQlBQEFOnjqewsIDI\nSC2IISKu4XBYLH9vH7sP5DGgZzseuH0wAf7nD/38ygIWJS2jvLaCuwZMZljHeA9V2zK1+OCf2G/M\neY/O3cGdq/OdeS87O4vevfsQGhpKeXnpj1bX69Hj3Pe+dunSFX//um95VFQUzzzzd0JDQ8nLO0V8\nfCIA3br1+GHFvg4doqiu1up8IuIaDsti5Qf7+XZvLv26RfDQ5HgCA/zOu01xdQmLkpZSXFPCxH5j\nGN71cg9V23LpGr+buGt1vjNtn1ldr7KyEvjx6nq+vr4/aufMgkBnv/63v/0Pv/vdf/Db3/6JqKjo\nHz5HK/OJiDs4LIs3Pj7A9rRsenUO5+EpCQQHnv/Ys6y2nEXJy8ivKmB0zC+4oee1Hqq2ZVPwu4m7\nVuc748zqejNmzDjv6nqDB8fz3//9J0pLS370+k03jeaBB2bz+OMPExHRjtOn889Rt4hI4zksi1Uf\nmWxNOkn36DAemZZIaPD5Q7/KVsWSlBVkl+dyfferua33jR6qtuXT6nzNnFao0v631v1vzfsOzWf/\nHZbFax/uZ1tKNj07hvHI9ETCQwPPu02NvZYlKcs5WHSEKztfyl0DJ+Pr8+Pj1Oay/+7SmNX5Wvw1\nfhER8Q6Hw+KVD/axIy2HXp3CeWR6ImEhAefdxu6wsyJ9FQeLjpAYPZg7B0z6l9CXxlHwi4iIy50Z\nvf91eg69u4Tzm2mJtAk+f+g7LAev7VtLWv4+BkTGcm/cnfj5nn/wn1w8Bb+IiLiU3eFg+bv7+GZv\nLn26tuU3UxMIvUDoW5bF2gOb2JmbTJ+IXsyNn0mAryLKHfRVFRERl7HZHSzbspfv95+iX7cIfj01\ngZCgC0fN5iMfsv3kN3QP68qC+PsI8jv/OABpOAW/iIi4xJkZ+XaZecR2j+DhKc6F/sfHtvJxxlY6\nhkbxYOJsQgNCPFBt66XgFxGRRrPZHbywaQ9JB/MxerTjV1PiL3ifPsC2zK9558gHRAa1Y2HiHMID\nwzxQbeum4BcRkUaptdWFfvKhfAb2iuShSfEEBV54UN53Obt568AmwgPCWDh0Du2DIz1QrSj4RUSk\nwWptdp7fuIfUw6eJi4nkwUnxBF1gGl6AtPy9vL7vLYL9g3gwcTadQqM9UK2Agl9ERBqoptbO4g1p\n7DlawOA+7Xnw9iEXnHsfwCw4xMt7VuHv48cvE+6je3hXD1QrZyj4RUTkolXX2ln0dip7jxUS37dD\n/Sp7Fw79o8XHeTFtJVgWcxPupU9EjNtrlR9T8IuIyEWprrHz7PoU9h8vIrFfFAsmXHhpXYCsshyW\npCyn1l7L7CH3MLB9fw9UKz+l4BcREadV1dh4dl0q5okihvWPZv74OPz9Lhz6pyryWZS8jApbJfcM\nnEpi9GAPVCvnouAXERGnVFbbeGZdCgczi7nUiGbuOOdCv7CqiEXJyyipKWVK7Hiu7HKpB6qVn6Pg\nFxGRC6qstvHUW8kcPlnC5QM7MmfsIPx8Lxz6pTVlPJe8lIKqQsb2uZnre1ztgWrlfBT8IiJyXhVV\ntTz1VgpHskq4clAn7h8z0KnQr6itYFHyMk5V5HNjz+u5udcoD1QrF6LgFxGRn1VeVcuTa5I5llPK\nVXGduf+2gfj6Xngp+CpbNUtSXuFkWTYjul3J+L6j8fFp8BLy4kIKfhEROaeKKhtPrE4mI7eUEUO6\ncO/oAU6Ffq29lqVpr3K0JIPLOg1jWv8JCv0mRMEvIiL/oqbWznPrU+pCP74+9J0Ib7vDzor0NzEL\nDxEfFcc9A6fg63PhywLiOfpuiIjIj5xZcOdAZjGXDujIvbc4F/oOy8Fr+9aSmp/OgMhY7ou7Ez/f\nC0/qI56l4BcRkR84LItX3t9HSv3c+3PGDHLq9L5lWaw1N7IzN5k+Eb2YGz+TAL8AD1QsF0vBLyIi\nQF14r/n0IF+n59K3a1semDjEqRn5LMti0+H32Z71Ld3DurIg/j6C/AI9ULE0hIJfREQA2LLjGJ/u\nyqRbVBt+NSWB4EDnhoF9lPE5nx7/kk6h0TyYOJvQgBA3VyqNoeAXERE+25XJpu1HiYoI5jfTEgkL\nce40/dYT29ly5CPaB0eyMHEO4YFhbq5UGkvBLyLSyn2TnsMbnxygbZtAHpmeSGR4kFPbfZ31PesP\nbqZtYDgLE+cQGdzOzZWKKzTZ4DcMY5RhGMu8XYeISEuWejif5e/tIyTIn99MTaBTZKhT2+0+lcob\n+9fTxj+UhYlz6Bga5eZKxVWaZPAbhtEXSASCvV2LiEhLdeBEEc9v3IOfrw8PT4mnZ6dwp7ZLP72f\nlemrCfIL5IHE++ka1tnNlYorNcngN03zsGmaT3m7DhGRlup4binPrk/B4bD45e1DiO3u3Gn6g4VH\nWJb2Gr4+PsyPn0Wvtj3cXKm4msdm7jMM4wrgr6ZpjjQMwxdYAsQD1cBs0zQPG4bxX0A/YIFpmkWe\nqk1EpDXJLajgqbXJVFXbmTNuEPF9Ozi1XUbJCV5MfQWHZTEvfiaxkX3cXKm4g0eC3zCMx4G7gbL6\nlyYAgaZpDq//g+BJYIJpmn/wRD0iIq1VYWk1T6xJpqSilrtv6s+Vg5w7TZ9VlsPzycupttdw3+C7\niOswwM2Virt46lT/IWAicGb6pxHAhwCmaX4LXHqujUzTvMcj1YmItAJllbU8uTaZ0yVVTLimN6OG\ndXdqu1MV+SxKXka5rYK7BkxmWMd4N1cq7uSRI37TNDcYhhFz1kvhQMlZz+2GYfiapuloSPvR0c4N\nSGmptP/a/9aqNe87XNz+V1bb+L83k8jKL2fcNX24b/xgp1bMO11RyJJvXqakppR7h07h1v6jGlOy\nS7X2739DeWt1vhLqwv+MBoc+QF5eaeMraqaio8O1/9p/b5fhFa153+Hi9r/W5uDZ9SmYxwu5Kq4z\n44b3Ij+/7ILbVdRW8MSuJeRVFDCm981cFnlZk/ma6/vf8D96vDWqfwdwK4BhGFcCqV6qQ0SkRXM4\nLJZtSWfvsUIS+0Ux61bnl9ddtmcVuRWnGNXjGm6JaTpH+tI4nj7it+o/bgRuNAxjR/3zWR6uQ0Sk\nxbMsi9c+2s9OMw+jRzsWTIjD38+5RXfWmBs5UHiIhKg4bu93m1OXBaR58Fjwm6Z5DBhe/9gCFniq\nbxGR1mj9l4fZlpJNr07hPDQ5ngB/P6e2++zENv6R/R09wrsxM+4OfH2a5JQv0kD6boqItEAffJPB\nB98cp1P7UH49NYGQIOeO81Ly9rDp0PtEBLZlfvy9Wl63BVLwi4i0MNtSslj3xWEiw4N4dFoibds4\nF97HSzJZmb6aAF9/FiTMol1QhJsrFW9Q8IuItCC7zFO8+uF+wkICeGRaIh0inFvypLCqiBdTX6HW\nYWNW3J30CO/m5krFWxT8IiItRPqxAl7anE5ggB+/nppA16g2Tm1XZavmxdSVFNeUMqHfrcRHx7m5\nUvEmBb+ISAtwJKuExW+nAfDQxCH07tLWqe0cloOVe1eTWZbF1V2v4IYe17qzTGkCFPwiIs1cVn45\nz6xLocZmZ964wQyMae/0thsPvUda/l4GRMYyrf8E3bbXCij4RUSasdPFVTy5NpmyylruvWUAlxjR\nTm/71clv+PzEV3QK7cj9g+/Gz9e52/2keVPwi4g0UyUVNTyxNpnC0mqmjOzLNQldnd52X8EB3jqw\nibCANvwyYRahASFurFSaEgW/iEgzVFFVy9NrU8gtqGD0lT0ZfUUvp7fNKc9l+Z5V+OLD3CEziQrp\n4MZKpalR8IuINDO1Njt/WfEdGbmlXBPfhcnX9XV629KaMpakvEKlrYq7Bk6hb7sY9xUqTZKCX0Sk\nGbE7HLz4Tjpph/O5pH80M24xnB6QV2uvZWnaq5yuKmB0zC+4vPMwN1crTZGCX0SkmbAsi1c/MEk6\nmE9CbBRzx8Xh5+vcr3HLsli1fx1HijO4pGMCt/W+0c3VSlOl4BcRaQYsy+KtrYfYnpZN7y7h/Pbe\nywnwd/5X+PvHPmVnbjK92/binoFTddteK6bgFxFpBt7/JoOPvjtBlw6hPDwlgdDgAKe33ZmTxPtH\nP6FDcCTz4mcS4Of8ttLyKPhFRJq4L5NP8vaXR2jfNohHpiUSHur8inlHio/x+v51BPsFMz9+FuGB\nYW6sVJoDBb+ISBO2c/8pXvvI/GHRnfZtnVt0ByC/soCXUl/FYTmYPfhuuoZ1dmOl0lwo+EVEmqj0\nYwUs3VK36M5vpiXQpYNzi+4AVNRW8kLKCspqy5kSO56BHfq7sVJpThT8IiJN0D8X3fHhoUnxxHR2\nbtEdALvDzvI9q8ipOMXIHiO4tvtV7itUmh0Fv4hIE3Myv5yn30qmxmZn/vg4BvaKdHpby7J46+A7\n7C88yOAOA5nYb4wbK5XmSMEvItKE5BdX8tTaZMqrbNw7egDD+ju/6A7A9qxv2H7yG7qFdWFW3B34\n+ujXvPyY/keIiDQRJeU1PLmmbtGdqSP7cU2884vuABwpzmDdgc2EBbRh3pB7CfZ3fiCgtB4KfhGR\nJqCy2sZTbyWTW1jJrVf24pYrel7U9sXVpbyc9joOy8GsuDvpEOL85QFpXRT8IiJeVmuz89z6VI7n\nlnFtQlcmXdfnora3O+ysSF9FcU0J4/uOZkD7WDdVKi2Bgl9ExItqbXYWbUjDPFHEJUY0M252ftGd\nMzYefo9DRUdJjB7CL3pe56ZKpaVQ8IuIeEmtzcHzG/ew50gB8X07MHdsHL6+Fxf63+cksfXEdjqH\nduSegVM0B79ckIJfRMQL6kI/jdTDpxnSpwMP3D74ohbdAThZls0b+9cT7BfE3CEzNJhPnKLgFxHx\nMJvdwQub9pB6+DSDe7fnwYmDCfD3u6g2ymrKWZr6KrWOWmYMmkanNh3dVK20NP7eLkBEpDU5E/rJ\nh/KJi4nkwYlDLjr0HZaDRd+8Sn5VAbf0GkVC9GA3VSstkY74RUQ8xGZ38OI76SQdzGdgr0gWToon\nMODiQh/g/aOfkpSdzsD2/bmtz01uqFRaMgW/iIgH2OwOXtqczu4DeQzo2Y6HJjcs9NPy9/LBsU/p\n2KYDs+Lu1Mx8ctF0ql9ExM3sDgdLt+xll1kX+r+anEBQA0L/VEUeK9PXEODrz6NXz6ONLdQN1UpL\npz8VRUTcyO5wsGzLXnbuP0X/HvWhH3jxoV9lq2ZZ2utU2au4c8BkYiJ7uKFaaQ0U/CIibuJwWCx/\ndx/f7TtFbPcIHp4S36DQtyyLN/evJ6s8h+u6X83lnYe5oVppLRT8IiJu4HBYLH9vL9/szaVftwge\nnpJAcGDDrq5+fuIrdp1KoU9EDBP73ebiSqW1UfCLiLiYw2Gx4v19fJ2eS9+ubfn11ARCghoW+gcK\nD7Hp8Pu0DQxn9uC78ffV0CxpHAW/iIgLOSyLVz7Yxz/25NC7S1t+PTWxwaFfWFXE8j1vADB78D1E\nBLV1ZanSSin4RURcxGFZvPrBfnak5dC7SziPTEsgNLhhoV9rr2VZ2uuU1ZYzOXYcfdvFuLZYabUU\n/CIiLuCwLF770OSr1Gx6dQ7nkWmJhAYHNLi9dQffIaP0BFd0voRru13lwkqltVPwi4g0ksOyWPWR\nybaULHp1CufR6Y0L/R0nv2VH1nd0D+vKdGOiVtwTl1Lwi4g0gmVZvPHxAb5IzqJnxzAemZ5Im0aE\n/rGS47x1YBNt/EOZM2QGgX4Nb0vkXBT8IiINZFkWb35ykK1JJ+nRMYxH7xhKWEjDg7q0poxlaa9j\ntxzMiruTqJD2LqxWpI6CX0SkASzLYvWnB/lsdybdo9vw6PTERoW+3WFnxZ43KKouZmyfmxnYob8L\nqxX5JwW/iMhFsiyLNZ8d4tNdmXSLbsOjdwwlPDSwUW2+c+QDDhQdJiF6MDf1GumiSkX+1QXvMzEM\now8wBogFHMBBYItpmhlurk1EpMk5E/qf7DxB16g2PDZ9KG0bGfr/yPqez45vo1NoNPcMnKrBfOJW\nPxv8hmF0BZ4GYoDt1AV+LdAHeMswjGPAI6ZpZrq9ShGRJsCyLFZ/dpBPd2bSNaoNj98xlLZtGhf6\n3+Xs5s3962kTEMrcITMI8Q92UbUi53a+I/7/Bf5smubec71pGEYC8FfgbncUJiLSlJwd+t2i2vCY\nC0J/96lUXtu7lmD/YBYmzqFzm04uqlbk5/1s8JumOfN8G5qmmYJCX0RaAXeEfkpeOq+kv0mQXyAP\nJt5Pj/BuLqpW5PycucZ/BTACWAxsAYYB803TXO/m2kREvO7M6P0zA/kem9740N+Tv4/le1bh7+vP\nLxPuJ6ZtTxdVK3Jhzozqfw7YCUwCKqkL/n93Z1EiIk2BO0J/X8EBlu15HV8fXxbEz9Ic/OJxzgS/\nr2maXwK3AW+bpnkc8HNvWSIi3mVZFm+eHfouOL1/sPAwL6W+CsC8+Jn0j+zrilJFLoozy0ZVGIbx\nKHADsNAwjF8Bpe4qyDCMG4BpQCjwN9M0U93Vl4jIuZwJ/c/ODv1G3rJ3uOgYS1JfwWE5mDtkBgPb\na4Ie8Q5njvjvoi6EJ5qmWQB0Bu50Y00hpmnOBZ4AbnJjPyIi/+LMNLyuDP1jJcdZkrIcm8PG/YPv\nYnDUQBdVK3Lxzncf/1jTNLfU36f/n2deN03z/531OeNN03zHlQWZpvmuYRhtgIeAx13ZtojI+fwQ\n+mem4XVB6B8vzWRx8nKq7TXcN/guEqIHu6hakYY536n+3oZhfAKsA7YBmYAN6AWMAqYDG53ppP7O\ngL+apjnSMAxfYAkQD1QDs03TPGwYxn8B/YBfUTc/wB9N08xv2G6JiFwcy7J445MDfL77pMtC/2RZ\nNouTXqbKVsWMQdMY1jHeRdWKNNzPnuo3TfM56u7T7w6sBnKAU8AaoAsw1TTNZy7UgWEYjwPLgKD6\nlyYAgaZpDqfu7oAn6/v7g2madwB/BzoB/2sYxqQG7peIiNPcEfo55bk8l7SUclsFdw6YzOWdh7mo\nWpHGOe/gPtM0c4E/1v9rqEPAROD1+ucjgA/r2//WMIxLf9LneScOEhFxpZ+G/mMuWHDnVEUezyUt\npay2nOnG7QzvepmLqhVpPLevzmea5gbqLhGcEQ6UnPXcXn/6X0TEoyzLYtUPoR/mktDPryzg2aSl\nFNeUMjl2HNd0u8pF1Yq4hjO387laCXXhf4avaZqOxjQYHR1+4U9qwbT/2v/WqjH7blkWL25IZevu\nk8R0actf5g8nIizowhueR355AYu/XUZRdTF3J9zOuAHuvTGpNX/vQfvfUN4I/h3AWGCdYRhXAo2+\nTz8vz23TCjR50dHh2n/tv7fL8IrG7LtlWaz6+ABbk07So2MYv54ST01lDXmVNQ2up6i6mKd3v0h+\n5WnG9L6Jqzpc5dbvTWv+3oP2vzF/9DgzV3974P+oG3E/Ffgb8BvTNAsvsi+r/uNG4EbDMHbUP591\nke2IiDSYw7J446zQf3R6YqNP75fUlPJc0jLyK09zS69RjO79CxdVK+J6zhzxLwM+Bq6gbsa+k8Aq\n6qbwdYppmseA4fWPLWDBxRYqItJYPw39x+4YSlhIQKPaLKspZ1HSMnIrTnFDz2sZ0+dmF1Ur4h7O\nDKrrbZrmS4DdNM0q0zR/D/Rwc10iIi51ZvS+K0O/vLaCRcnLyCrP4fruV3N739vw8fFxUcUi7uFM\n8NcahhFx5olhGLGA3X0liYi4lmVZvLX1EFt3uy70K22VLE5+mcyyLEZ0vYLJseMU+tIsOHOq/0/A\nF0BPwzDeAa4C7nNnUSIirvTO9qN89N0JunQI5ZHpiY0OfbvDztK01zlemsmVXS5lmnG7Ql+ajQsG\nv2maHxqGsQu4nLrleOfWT+wjItLkffBNBpt3HCO6XTCPTm/8jHyWZfHWgU0cKDxEfFQcdw2YjK+P\npiKR5sOZUf0dqZuXP7L+paGGYVimaf7neTYTEfG6z3Zlsu6Lw7RvG8Rj04cSGd64+/QBvsjcwfas\nb+kW1oVTOyU9AAAgAElEQVSZg6Yr9KXZceZ/7PtA4lnPfer/iYg0WV+lZvHGJwdo2yaQR6cPJapd\nSKPbTD+9n7cPbqFtYDgL4mcR7N/4PyREPM2Za/yWaZq6pi8izcZ3+3JZ+cF+wkICeHR6Ip3bhza6\nzayyHFbseQN/Xz/mxc8kMridCyoV8Txngn+TYRhzgM84a8590zSPu60qEZEGSjqYx7ItewkO9OM3\n0xLoHh3W6DZLa8p4MfUVquzV3Bd3JzFte7qgUhHvcCb4I6hbPjf/J6/3dn05IiINt+foaV7YtAc/\nPx8enpJATOe2jW6z1mFjadprnK4q5NbeN3JJp8QLbyTShDkT/JOBjqZpVrq7GBGRhjKPF7L47TTA\nh4cmxRPbvfGn4i3L4s396zlSfIxLOiZwa4ym4pXmz5nBfYeB9u4uRESkoY5klfDM+lTsDosHbh/M\noBjX/Mr6JOMLvsvZTUzbntw9cKru1ZcWwdnV+fYahrEHOLN0lWWa5ig31SQi4rTjuaU8tTaZmlo7\nC8YPJqFflEvaTc7bwztHPiAyqB1zh8wk0K9xk/6INBXOBP9/n+M16xyviYh4VPbpcp5cm0xFtY3Z\nYwZy6YCOLmn3ROlJXk1fTaBfIPPj7yUiSOu+S8vxs6f6DcO4pP6hBTjO+meh4BcRL8s5Xc7fVydR\nWlHLPTcbDB/cxSXtFlUX82LqSmodNmYNuoPu4V1d0q5IU3G+I/75wBzgz5w76Ee6pSIRkQsoKKni\nb2uSKSqrYdqofowc2s0l7dbYa3gp9VWKqouZ0PdW4qPjXNKuSFNyvuA/CmCa5vWeKUVE5MKKy2v4\n+5pkThVUMOGa3tx8uWvuqXdYDl7b99YPC+/8oud1LmlXpKk536j+KR6rQkTECWWVtTyxJoncggom\njezH2OExLmv7/aOfkHQqlb4RvbnDmKgR/NJiaXUJEWkWKqpsPLk2mZN55dwwrDszbxvksnD+PieJ\nD459RlRwe+YOmYG/r7M3PIk0P+f7351gGIbjZ96zTNP0c0dBIiI/VV1j55n1KWTklDJiSBfuuDHW\nZaF/pDiDVfvXEewXzPyEWYQFtnFJuyJN1fmCP8U0zaEeq0RE5BxqbXYWbUjlUGYxlw/syL2jB+Dr\notA/XVnI0tRXsTvszEuYSZc2nVzSrkhTplP9ItJk1drsPL9xD3uPFTI0NorZYwbh6+ua0K+yVfFS\n2kpKa8uY3H8cgzoYLmlXpKk73xH/Oo9VISLyE1U1Nha9nca+jELierdn/vjB+Pu55ljFYTlYuXc1\nJ8uyubbbVVzf/WqXtCvSHPxs8Jum+T+eLERE5IyKqlqeWZfKoZPFJPaLYsGEOAL8XXeCctPh90nL\n38eAyFgmx45zWbsizYGGropIk1JSUcNTa5I5fqqMKwd14r7bBrrsSB/gH1nf8dnxbXQKjeb+wXfj\n56txytK6KPhFpMkoLK3miTVJZJ+u4LrErtxzk+Gya/oABwsPs9rcQBv/UObHzyI0IMRlbYs0FxcM\nfsMwYoAHqVua98xPoGWa5n1urEtEWplTRZU8sTqJ/OIqbr68B1NH9nPpJDpHio/xYupKAGYPuYeO\noa5ZxU+kuXHmiP8tYFv9vzO0SI+IuMzJ/HKeWJNEcVkNE0b0ZuzVMS4N/UNFR1mSsrxu4Z24O+kf\n2ddlbYs0N84Ev79pmo+6vRIRaZUyckp5cm0yZZW1TB/Vj5tcNPf+GQcLD7Mk9RVsDhv3x91FYsch\nLm1fpLlxZsTMdsMwxhmGEej2akSkVTmYWcTfVu+mvLKWmbcYLg/9A4WHWJKyArvDzuzB9yj0RXDu\niH8Kddf4MYwfJrjQlL0i0ijpRwtYtCEVu91izrhBXDmos0vb319wkBdTV2JZDuYMuYchUYNc2r5I\nc3XB4DdNs4snChGR1iPpQB4vvLMH8OGB24eQGOvagXb7Th/gpbSVWMCcITMYHDXQpe2LNGfOjOpv\nA/wJuKH+8z8Hfm+aZrmbaxORFujr9ByWv7sPf38fHpoUz6CY9i5tP/30fpamvQbAvCEzNRWvyE84\nc41/MRAKzAJmAoHAi+4sSkRapi+STvLylr0EBfrx6LShLg/9tPy9LE19FR9gfvy9Cn2Rc3DmGv8l\npmnGn/X8AcMw9rmrIBFpmT789jhvbT1EWEgAj0xLpFfncJe2n5qXzst7VuHr48uC+FkY7fu5tH2R\nlsKZI34fwzAizzypf1zrvpJEpCWxLItNXx3hra2HiAwP4v/dPczloZ+ct4dle17Hz8eXXybcp9AX\nOQ9njvifAr4zDGMzdTP3jQP+161ViUiLYFkWaz47xCc7TxAVEcxjdwwlup1rp8ndfSqVV9LfxN/X\nn1/G30dsZB+Xti/S0jgzqv8VwzB2AtdSd4bgdtM009xemYg0aw6HxWsf7WdbSjZdOoTy6PShRIYH\nubSPXbnJrNy7hgBffx5ImE3fdjEubV+kJfrZU/2GYYyt/zgTGAqUASXAMMMwZnimPBFpjmx2B0u3\npLMtJZtencL597uGuTz0t2d8xyvpqwn0DeDBxDkKfREnne+I/1JgCzCSc8/N/5pbKhKRZq3WZmfJ\nxj2kHD5Nv+4RPDw5gdBg1y4E+m32Ll7f/xZBfkE8mDib3hGunfFPpCX72Z9G0zT/VP/wTdM0Pz77\nPcMwJrm1KhFplqpr7CzakMreY4UMiolk4cR4ggJdO8nn19k7eWPfOkIDgnkgYTa92vZwafsiLd3P\nBr9hGNOBIODPhmH88ay3AoDfAm+7uTYRaUYqq208sy6Fg5nFJPaLYsGEOAL8XRv6/8j6jjf3v02I\nfzB/uP5hwu2RF95IRH7kfOff2gLD6z+OPOt1G3XBLyICQFllLU+/lczR7FIuHdCRuWMH4e/nzN3C\nztt+8htWmxtoExDKwsS59Gnfk7y8Upf2IdIanO9U/1JgqWEYN5im+ZkHaxKRZqSkvIYn1yZz4lQZ\nwwd3ZtatA/DzdW3ob8v8mrUHNhIW0IaHhs6lW5iWEBFpKGdG3PzeMIzf/+Q1yzTNUe4oSESaj8LS\nap5Yk0T26QquH9qNu2/qj6+Pj0v7+CJzB+sOvENYQBt+NXQeXcNcu4qfSGvjTPD/+azHAcB4oNA9\n5YhIc5FfXMkTq5M5VVTJTZf1YNqofvi4OPQ/yfiCTYffJzwwjF8NnUeXNp1c2r5Ia+TMBD5f/OSl\nTwzD+A74g1sqEpEmL7ewgidWJ3G6pJoxw3tx+zV9XBr6lmXx/rFPef/oJ7QLiuChxDl0atPRZe2L\ntGbOLMt79g2yPsBgwLVLaolIs3Eyv5wn1iRRXFbDxGv7MGZ4jEvbtyyLTYff59PjX9IhOJKHhs4j\nKkS/ckRcxZlT/dv45wQ+FpAPLHRbRSLSZB3PLeWJNcmUVdYy/YZYbrrMtffQOywH6w5sZtvJf9Ax\nNIqHEucSGdzOpX2ItHbOnOqPMQwjwDTNWsMwAoFA0zTLPFCbiDQhR7JKeGptMpXVNmbcbHD90G4u\nbd9hOXhj/3q+yd5J1zadWTh0Dm0DXbuKn4g4sSyvYRhTgd31T3sC+w3DmODWqkSkSTlwoogn1iRR\nWWPjvtsGujz07Q47K9NX8032TnqGd+dXw+Yp9EXcxJlT/X8AfgFgmuYhwzCGAZ8Am9xRkGEYlwAP\nUjee4HHTNE+5ox8RcU76sQIWrU/F7rCYP34wlw1w7SC7WoeNFXveIDU/nT4RMfwyYRYh/q5duldE\n/smZWTYCTNPMPfPEA0EcBDwMvAdc5ea+ROQ8kg/l8+y6VByWxQO3D3F56NfYa3gpdSWp+en0j+zH\ng4mzFfoibubMEf8OwzBWA29QdxQ+FfjaXQWZpvkPwzCuAh6t70tEvGDn/lO8tDkdP18fFk5KIK63\na0fWV9mqeCH1FQ4VHSWuwwBmD76HQL8Al/YhIv/KmeB/gLpR/POAWupG+S+5mE4Mw7gC+KtpmiMN\nw/Ct3z4eqAZmm6Z52DCM/wRigaeAncBo4E/Ary6mLxFpvK/35PDye3sJDPDj4cnxGD1duxhORW0F\nz6es4FjJcRKjhzAr7g78fV27dK+InJszo/qrDMNYD+wDPgJ6mKZZ42wHhmE8DtwNnLkTYAJ1dwYM\nr/+D4Elggmmaf6z//JHACqAGeOlidkZEGu+L5JO8/qFJSJA/v56WQN+uES5tv7SmjMXJL5NZlsVl\nnYZxz8Ap+Pm6dhU/Efl5zkzgMx34HRAKXE3dqf/HTdN83ck+DgETgTOfPwL4EMA0zW8Nw7j07E82\nTXMrsNXJtkXEhT75/gSrPztIWEgAj05PpGcn146sL64u4bnkZeSU53J11yuYbtyOr49rF/QRkfNz\n5ifu36gL/BLTNHOAYcD/c7YD0zQ3ULeU7xnhQMlZz+31p/9FxIve+/oYqz87SESbQP7trmEuD/2C\nqkKe3v0COeW5jOwxgjuMiQp9ES9w5qKa3TTNEsMwADBNM9swDHsj+iyhLvzP8DVN09GI9oiObt33\n+2r/tf+N4XBYvPreXjZ8eYSodiH89/zhdI0Oc1F1dXJKT/HsNy+RX1nAxEG3MG3wOJfM7a/vvfZf\nLp4zwZ9uGMZCINAwjETgl0ByI/rcAYwF1hmGcSWQ2oi2AMjLK21sE81WdHS49l/73+DtbXYHK97f\nxzfpuXRqH8ojUxMIwHLp1zS7PJdFSUsprillbJ9buKHzKPLzGz/5p7732v/Wvv8N5eyo/t8DldQN\nuvsceKQBfZ2Z738jcKNhGDvqn89qQFsi0kiV1TYWb0hjX0Yhfbu25aHJ8YSHBrq0jxOlWSxOXkZZ\nbTmTYscyqsc1Lm1fRC6eM8HfxzTNfz/7BcMwJgPrne3ENM1jwPD6xxaw4CJqFBEXKyyt5pl1KZw4\nVUZivyjmjY8jKMC1I+uPFh/n+ZTlVNmquMOYyIhuV7q0fRFpGGdG1myuvyUPwzA6GIaxlrpR/iLS\nDGWfLud/Xt/FiVNlXJ/YlQcmDnZ56B8sPMKi5KVU2aqYMWiaQl+kCXEm+IcB8YZhfA18C3wHXHr+\nTUSkKTqUWcz/vL6L0yVV3H5Nb+652cDP17Uj6w8UHub5lOXYHHbuH3w3l3ce5tL2RaRxnDnV70vd\njH2h1E3ZawcaNQpfRDxv94E8Xtqcjt1uMevWAVwT39XlfRwqOsoLKStwWA7mDLmHIVGDXN6HiDSO\nM3/q7wEygEuAK6m7Vv+dO4sSEdfaujuT5zem4evjw0OT490S+keKj7EkZTk2y87swXcr9EWaKGeO\n+G81TXN3/eM8YKphGFPcWJOIuIhlWWzYdoT3vs4gPDSAh6ck0LtLW5f3c6zkOM8nr6DWYeO+uLuI\nj45zeR8i4ho/e8RvGMYCANM0dxuG8dOf4hFurUpEGs1md7DivX2893UGHSND+N09l7gl9I+XZrI4\neTnV9mruHTSdoR2HuLwPEXGd853qn3vW41U/ee9aN9QiIi5SWW3jufWp7NiTQ+8u4fz27kvoGBnq\n8n4yS7NYlLTsh9H7l3RKdHkfIuJaWgdTpIUpLq/hmXUpZOSUEt+3AwvGDyYo0PWr32WV5bAoeRmV\ntiruHjhFo/dFmgkFv0gLklNQwVNrk8kvruKa+C7MuMX1t+sB5JTn8lzSUspqy7lzwCSu7KI7fEWa\nCwW/SAtxOKuYZ9elUlZZy7irYxg/ordLFsL5qdyKPJ5NWkppbRnT+t/O1V2vcHkfIuI+5wv+OMMw\njtY/7nrWYwDX3wskIg2WfCifFzftodbuYOYtBtcldnNLP3kVp3kuaSklNaVMiR3Ptd2vcks/IuI+\n5wv+/h6rQkQa7Mvkk7z2kUmAny8LJ8aTGBvlln7yKwt4NukliqqLmdhvDNf3uNot/YiIe/1s8Ncv\nrCMiTZRlWbz50X5Wf2wSFhLAr6bE07drhFv6Kqgq5LmklyisLmJ839Hc0FM39og0V7rGL9IM1drs\nvPqhyT/25BAVEcxvpiXSub3rb9cDKKou5tmkpZyuKmRM75u5qddIt/QjIp6h4BdpZgpLq1m8IY2j\n2SXE9mjHL8fHEREW5Ja+iqtLeHb3S+RXnmZ0zA2M7n2DW/oREc9R8Is0I4dPFrN4YxrFZTVcFdeZ\nR++5lOKiCrf0VVJTyrNJSzlVmc9NvUZyW++b3NKPiHiWgl+kmdiems1rH+3H7rCYPqofN17Wg8AA\n10/MA1BaU8ZzSUvJrTjFDT2uZVyfW9xya6CIeJ6CX6SJszscrP38EJ/uzCQ0yJ/5E+IY3LuD2/or\nr61gUfIysstzGdl9BLf3u02hL9KCKPhFmrCyylpe2LSHfRmFdOkQykOT4+nkhjn3z6ioD/2TZdlc\n0+0qJsWOVeiLtDAKfpEmKjOvjEVvp5JXVEVivyjmjB1ESJD7fmQrbZUsTl7OidKTDO9yOVP7j1fo\ni7RACn6RJmiXmcfL7+6lutbO2OExjL+mN75uDOEqWxXPJ68go/QEV3a+lDsGTMTXx/Vz/IuI9yn4\nRZoQh2Xx7o5jbNp+lMAAXxZMGMxlAzq6tc+ymnKeT1nO8dJMLus0jLsGTlboi7RgCn6RJqKqxsby\nd/ex60AeHdoGs3DSEHp2Cndrn4VVRSxOfpmcilNc1eUy7jB0pC/S0in4RZqAU0WVLH47lcy8cgb0\nbMf8CYNpGxro3j4r8lmUvIyCqkJG9biGif3G6Jq+SCug4Bfxsn3HCliyaQ/lVTZuGNadaTf0w9/P\nvUfdmaVZLE55mdKaMsb2uZmbe41S6Iu0Egp+ES+xLIvPdmWy5rND+Pjg1uV0z3ak+BhLUl6h0lbJ\ntP4TuLb7cLf3KSJNh4JfxAtqbQ5e/9hke2o2bUMDeGDiEGK7t3N7v3tPmyxNew27ZWfmoOlc3nmY\n2/sUkaZFwS/iYUVl1Ty/MY3DJ0vo1TmchROH0L5tsNv73X0qlZXpq/H18WHukBkMiRrk9j5FpOlR\n8It40JGsEp7fmEZhaTVXDurEvaMHuG2+/bPtyPqW1fs3EOQXyPz4e4mN7Ov2PkWkaVLwi3iAw2Hx\nwbcZbPrqKA6HxZSRfbnl8p4eGVD3ScYXbDr8PmEBbXgg4X56tu3u9j5FpOlS8Iu4WX5xJS+/u48D\nJ4poFxbI7DGDGBTT3u39WpbF5iMf8nHGVtoFRbAwcQ6d27h3MiARafoU/CJu9M3eHF7/6ACV1TYu\n6R/NzNEDCAsJcHu/DsvB2gOb2H7yGzqGRPFg4hw6hES6vV8RafoU/CJuUFFlY9UnJt+k5xIU4Mes\n0QMYEd/FI6f27Q47r+1by87cZLqHdeWBxPtpG+jeGQBFpPlQ8Iu42IETRSzbspfTJVX07tKWueMG\nuXUp3bPV2Gt4ec8q0k/vp09EDAviZxEaEOKRvkWkeVDwi7iIze5g845jvPf1MQDGDo9h7NUxbp+F\n74xKWyUvpKzkcPFRBrU3mDPkHgL93Dvtr4g0Pwp+ERfILahg6Za9HM0uISoimNljBtG/h/sn5Dmj\ntKaM55Nf5kRZFpd0TGDGoGn4++rHW0T+lX4ziDSCZVl8lZrN6k8PUl1r56q4Ttx1o0FosOd+tAqq\nClmUvIxTFflc3fUKphu3a4U9EflZCn6RBiqrrGXlB/vZfSCPkCB/5o2L44pBnTxaQ1ZJDk/teoHC\n6iJu7Hk94/uO1mI7InJeCn6RBkg/VsDyd/dSVFZD/x7tmDNmEB0i3D/t7tmOl2TyQtoKSqrLGN93\nNDf1GunR/kWkeVLwi1yEWpudt788wsffn8DP14dJ1/Vh9BW98PX17FF2al46r6S/Sa3Dxh3GREZ0\nu9Kj/YtI86XgF3HSybwyXtq8l8y8Mjq1D2Xu2EH07tLWozVYlsXWzO1sOPguAb7+PDpiHjGBfTxa\ng4g0bwp+kQuwLIvPdmXy1tbD2OwOrkvsyvRRsQQFun9xnbPZHXbWH9zMtpNfExEYzvz4WVzSbSB5\neaUerUNEmjcFv8h5ZJ8u581PDpB+rJCwkABmjY5jaP9oj9dRZatiefob7D1t0rVNZ36ZcB+RwZ67\nXVBEWg4Fv8g5VFTZ2LzjKJ/tysTusBjcpz333TqQdmFBHq+lsKqIF1Jf4WRZNoPaG9w3+C5C/D07\nkFBEWg4Fv8hZHJbF9tRsNnx5mJKKWqIigpk2KpZh/aO8cpvc8ZJMXkx9heKaUq7tdhWTY8fh5+vZ\nSwwi0rIo+EXqHcos5o1PD5CRU0pggC+3X9uHWy7vQYC/d4I2JS+dlfUj9yfFjmVk9xG6R19EGk3B\nL61eYWk16744xDfpuQBcOagTk6/vS/u23jmdblkWW098xYZD7xHg68+cITNIiI7zSi0i0vIo+KXV\nqrXZ+ei7E7z3dQbVtXZ6dQrnzhtjie3uvUFz5xq537Ntd6/VIyItj4JfWh3Lskg6mM/azw+SV1RF\neGgAd/wilhFDunh8Ip6zVdqqWFE/cr9bWBcWxM/SyH0RcTkFv7QqJ/PLWf3pAfYeK8TP14ebLuvB\nuKt7e3RRnXMprCpiScoKsspzGNTB4P64uwjWyH0RcQMFv7QK5VW1vPPVUT7ffRKHVXd73h03xNKl\nQxtvl6aR+yLiUU0y+A3D6AS8a5rmZd6uRZo3h8NiW0oWG7Ydoayylo6RIUy/IZaEvh2axAh5jdwX\nEU9rksEPPAYc83YR0ryZxwt589ODnDhVRlCgH1Ou78svLu1BgL/316q3LIvPT3zFxvqR+3OHzCBe\nI/dFxAOaXPAbhrEAWAU84u1apHk6ml3CB99ksNPMA+DqwZ2ZdH1fr8y6dy52h511BzfzlUbui4gX\neCT4DcO4AviraZojDcPwBZYA8UA1MNs0zcOGYfwnEAt0rH/vcsMwJpmm+bYnapTmrdbm4Pv9uXy2\n6yRHs0sA6N2lLXfeGEvfrhFeru6fKm1VrNjzBnsLNHJfRLzD7cFvGMbjwN1AWf1LE4BA0zSH1/9B\n8CQwwTTNP/5ku9cU+nIh+UWVbNh2mC+TsyitqMUHSOwXxQ2XdGdgTCS+Teh6eXF1Cc+nLK+bc18j\n90XESzxxxH8ImAi8Xv98BPAhgGma3xqGcem5NjJNc4YHapNmyLIszONFfLY7k6SD+TgcFm2C/bnl\nip6MHNqN6HYh3i7xX5yqyGNx8sucripkRNcrmNp/gkbui4hXuD34TdPcYBhGzFkvhQMlZz23G4bh\na5qmo6F9REeHN3TTFqG17H9VtY2tuzN5b/sRMnLq1qDv0zWCMSN6c83QbgQHNrkhKwAcLsjg6aQX\nKKkuY0rcbUyOu82lI/dby/f/XFrzvoP2v7Xvf0N54zdlCXXhf0ajQh8gL6+0cRU1Y9HR4S1+/3ML\nK9i6+yRfpWZTWW3Dz9eHywd25IZLunNVYnfy88soLa6kKX4V9hUcYFnaa9TYa5lu3M41na4iP7/s\nwhs6qTV8/39Oa9530P5r/xv+R483gn8HMBZYZxjGlUCqF2qQJs5hWew5cprPdp0k7chpACLaBHLj\npTFcP7TbDyP0m/I97ztzk3lt71p8gNmD7yax4xBvlyQi4tHgt+o/bgRuNAxjR/3zWR6sQZq4iqpa\ntqdm8/nuk5wqqgSgX/cIbhjWnUuMaPz9vH8PvjO2ntjO+oObCfYLZl78TPpH9vV2SSIigIeC3zTN\nY8Dw+scWsMAT/UrzUFlt48CJIpIO5vPN3hxqah0E+PsyIr4LNwzrTq/Ozec6nmVZbD7yIR9nbKVt\nYDgPJNxP9/Cu3i5LROQHTXM0lLRotTYHh08WszejkH0ZBRzNKsVh1Z0QiooIZuSwblwT35WwkAAv\nV3px7A47q80NfJ39PdEhHXgwcQ5RIe29XZaIyI8o+MXtHA6LjNxS9h4rYH9GIQczi6mx1Y3n9PGp\nm2hnYK9IBvWKxOgZ6dWlcRuqxl7D8j1vsOf0PnqGd+OXCfcTHhjm7bJERP6Fgl9czrIssk5XsD+j\nkL3HCjCPF1FRbfvh/W7RbeqDvj39e7Tz+pK4jVVeW8GLqSs5UnyMAZGxzBlyjybmEZEmq3n/xpUm\n43RxFXszCtiXUci+jEKKy2p+eC8qIphLB0QzoFckA3u1J6JNoBcrda3CqiIWpywnpzyXSzslcs/A\nqfj76sdKRJou/YaSi2azO8g+XUFGTilHsuqu1Z8qrPzh/bahAVw+sCODYtozsFdkk5xJzxVyynNZ\nnLycwuoiru9+NZNix+Lr0zzuOhCR1kvBL+dVa7OTmVdORm4px3NKycgt5cSpcmz2f865FBLkR2K/\nKAb2imRgTCTdoto06fvrXeFocQYvpLxCua2C8X1Gc2Ov61v8PotIy6Dglx9U19g5kVdGRn3AH88p\n5WR+OXaH9cPn+Pn60D06jF6dw+jVKZxendvSq3MYfr6t50h3T/4+Xt6zCrtl564BUxje9TJvlyQi\n4jQFfytVUWXjxKnSH0I+I7eM7NPlWP/MeAL9fYnpHE7PzuF1Id8pnG7RbZrNJDru8G32LlbtX4ef\njy9zh8xgSNQgb5ckInJRFPzNnN1hUV5VS2WVjYpqG5XVNirOfnzWa2ee5xdX/eiaPEBwoB+x3dvV\nH8XXHc137hDaqo7kL+STjC/YdPh9QvxDWBA/i77tYrxdkojIRVPwNwM70rJJP1bwo/A+E+ZVNfaL\nbq9NsD+DYiLrQ77uSD46MqRJrV3flDgsBxsPvcfnJ76iXVAEDyTcT9ewzt4uS0SkQRT8zcAn35/g\n+Km6Fd18fCA0yJ+QIH86tguhbXgQAb4+hNS/dua90OD6x/Ufz7weEuRPgL+O4p1VY6/hjf3r2Zmb\nTKfQjjyYeD/tgyO9XZaISIMp+JuB395zCWWVtYQE+RMc6Pej0eOtfWlKd8opz+XlPavILs+ld9ue\nzE+YRVhAG2+XJSLSKAr+ZiAwwI/2AX7eLqNV+S5nN6vNDdTYa7iu+3Bu7zeGAE3MIyItgH6TiZyl\nxl7L+oPvsCPrO4L9grh/8N0M6xjv7bJERFxGwS9SL7cij+V7VnGyLJvuYV25f/BddAyN9nZZIiIu\npVV6nAEAAA5YSURBVOAXAXblJvPG/vVU22sY0fUKJseOI8CveS0LLCLy/9u79+iqyjOP499cCYQA\nuSAJNKAiPnKJoIAIXgALWjuteCl1vLV1Rpxr1Vlj7cz8YafOmjXTi652XMtltSqtONNpdZSpWEdl\noN4obRExhPQpoFxEEUICCZBILmf+2DuYRgInI+fsnLN/n7VccvbZ77uf5+Qkz37fs89+k6HCL7HW\n3tnOU1ue5ZVdayjMK+TmSdcxo/KcqMMSEUkZFX6Jrb2H9/FI3TJ2tuxidHElfzrlRiqLT4k6LBGR\nlFLhl1hav6eWZfU/o62zjTlVM1l85iIK87JnuWARkb6o8EusdHR18PSWFax+9zUKcwv40sRrmVU1\nPeqwRETSRoVfYmNfayOPbHyC7S07qRxyCrfU3ERV8aiowxIRSSsVfomFt/bW8eP6n9La0cqsyulc\na1cxSFP7IhJDKvyS1Tq7Olm+9Res3PkyBbn53HDWYmZXzfiD2x6LiMSJCr9kraa2/Tyy8Qnead7O\nKUMquGXKTYwZWhV1WCIikVLhl6y0saGeH2/6Tw51HGbGqGlcZ1dTlF8UdVgiIpFT4Zes0tS2n6e3\nrGDdng3k5+bzx3Y1F46epal9EZGQCr9khfbOdlbufIX/2baSI13tjCup5rqzrqG6ZHTUoYmIDCgq\n/JLREokEtQ2beHLzz2lo3cfQgmIWn3kl51dNJzcnN+rwREQGHBV+yVgfHNrDw5uW8ubuTeTm5HJJ\n9UVcfuoChhQMjjo0EZEBS4VfMk5rRxvPb1vJqp2v0pnoxErPYPGZi3QzHhGRJKjwS8boSnTxm93r\neWbrczQfaaG8qJSbp3+RUwtP18V7IiJJUuGXjLCj+V1++vvlvNO8nYLcfP7otIUsGDuPMZVl7N3b\nEnV4IiIZQ4VfBrSWIwf5+dvP8/p7vyFBgnNG1nDVGZ+jfHBp1KGJiGQkFX4ZkDq7Onl51xpWvPMC\nrR1tVBWPYvGERVjZGVGHJiKS0VT4ZcDxxi38bPNy3j/0AYPzi/jChCu4eMxs8nLzog5NRCTjqfDL\ngLGvtYmntzzL+r215JDDnKrzuGL8ZygpHBp1aCIiWUOFXyK3r7WJlTtf5vX31tLe1cFpw8ay+MxF\njBtWHXVoIiJZR4VfIrPr4Pu8uH016/ZsoCvRRemgEXz+9MuYWXmO7ronIpIiKvySVolEgs373+bF\n7avZ1OgAjC6uZMHYucwYNU2f44uIpJgKv6RFV6KLt/bW8cKO1Wxv3gnA+OGncem4eUwuP0s34BER\nSRMVfkmp9q4Ofr17HS/t+CV7DjcAcHbFZBaOm8fpw8dFHJ2ISPyo8EtKtHa08squX7Fq56s0H2kh\nLyeP2VUzWTB2LpXFp0QdnohIbKnwy0l14MNmVu18lVd2/Yq2zjaK8gbx6bEXc0n1RYwYNDzq8ERE\nYk+FX06KDw7v5aXtv+TXu9fRkeikpHAol427nAvHnK9lckVEBhAVfvlE3jmwg5d2rGbD3joSJBg5\nuJwFY+cyq3I6BXkFUYcnIiK9qPBLv3Qlunj7wHY2NtRTu6+e3Yc+AGBcSTULx81j6sjJ+g6+iMgA\npsIvJ9Ta0UZ94++pbdhE3b7fcaj9MAAFuQWcXTGZ+dUXMGHEeH0lT0QkA6jwyzE1tDZS27CJjQ31\nbN7/Np2JTgCGFw7jgtGzqKmYiJVOoFDT+SIiGUWFX4BgCn9b8w5qG+qpbdjE++EUPkB1yRhqyidS\nUzGJ6pIxGtmLiGQwFf4Ya+too75x89Ep/IPthwAoyM1nSvlEaiomMqVior6GJyKSRQZc4TezqcD9\nwFbgR+6+OtqIskMikeBg+yEa25p4p3lHMIXftJWOo1P4JVww+jxqKiZhpWdQmFcYccQiIpIKA67w\nA+cB7wMdQF3EsWSMrkQXzUdaaGxrorG1iX1tTcG/2/aH/2/iSFf7H7SpHjqaKRWTqKmYSHXJGF2N\nLyISAwOx8L8K/ASoBO4Evh5tOANDZ1cnTR8eoLGtZ1Fv4uDGFna3NNDUtv/oBXi9FecPYdSQkZQV\nlVJWVEpl8SlMLj+L0qIRac5CRESilpbCb2azgH919/lmlgs8AJwNfAjc4u5bzeweYALw3wQj/v3p\nim+ge3TjE7yx5y0SJI75/LDCEqpLxlBWNILyojLKikYcLfJlRSMoyi9Kc8QiIjJQpbywmtldwI3A\nwXDTlUChu88JTwjuBa5097vD/WcTfMbfDnwz1fFlguGDhnH68FMpH/xRMS8rKqW8qJQJn6rmQGNb\n1CGKiEiGSMeIegtwNfB4+PhC4HkAd19rZjN67uzua4A1aYgrY1wz4fN9Phd8j16FX0REkpPyq7nc\n/b8ILtTrVgI093jcGU7/i4iISIpF8Rl6M0Hx75br7l2foL+ckSNLTrxXFlP+yj+u4pw7KP+45///\nFcVI+zXgswBmdj7wVgQxiIiIxFI6R/zdl6Q/DSw0s9fCxzenMQYREZFYy0kkjv0VMREREck+uqhO\nREQkRlT4RUREYkSFX0REJEay7pa4ZjYHuDV8eLu7H4gyniiY2SXAde6+JOpY0s3MPg1cCwwBvu3u\nsfnWiJlNB/4ayAHucvc9EYeUdmY2CnjW3WdGHUu6xX1lUzObBNwOFALfdffYLPJmZrcD0whue7/M\n3R883v7ZOOJfQlD4HyEoALFiZuMJ3gBxvUH/YHe/FfgucGnUwaTZIOAOYAUwO+JY0s7McoCvAdsi\nDiUqcV/Z9BbgXYJbmW6LNpT0cvfvE9S9uhMVfcjOwp/n7kcIfgGqog4m3dx9q7vfF3UcUXH3Z82s\nGLgNWBpxOGnl7q8DkwhWtXwz4nCi8OfAMuJ7D+tXCYrftwneA3EznmDG40ngSxHHEoXrgaeS2TGj\npvqTWeUPOGxmhcBoYHd00Z58SeaftZJc5bGC4A/f3e7eEGG4J1WSuc8EfgtcDnyDYNozKyT53l8Q\nbjvPzK5x96T+CGaCJPOfRpaubJpk/nuAw0ATWTSo7cff/Yvc/ZZk+syYFydc5e9hgulM6LHKH/B3\nBKv8ATwE/IBgyv/x3v1kqn7kn5X6kf+9wCjgX8zsmrQHmgL9yH0o8CjwHeCJdMeZKsnm7+7XuPtf\nAGuzrOgn+/PfRjDi/Rbwb2kOM2X6kf+D4X53AP+e7jhToZ9/94ck228mnRUmtcqfu79Bdt4NsL+r\nHN6U3vBSLtmf/5ejCS+lks19FbAqkghTq7/v/Wyb5k3255+tK5smm/86INt+/5N+77v79cl2mjEj\n/riv8qf845t/nHMH5a/845t/qnLP5BfrZK/yl2mUf3zzj3PuoPyVf3zzPym5Z3Lhj/sqf8o/vvnH\nOXdQ/so/vvmflNwz6TP+bnFf5U/5B+KYf5xzB+Wv/ANxzP+k5q7V+URERGIkk6f6RUREpJ9U+EVE\nRGJEhV9ERCRGVPhFRERiRIVfREQkRlT4RUREYkSFX0REJEZU+EVERGIkE+/cJ5LRzOxU4G3gUnd/\nqcf2bcDF7r7jE/a/DTjX3Rs/ST8nOMZY4AWgBZjv7gfD7V8B7gO292ryZwTLhn7D3eenKq6+mNmq\nKI4rMhCp8ItEox142MxquosmH92W85NKADknqa++zAPWufsNxzj2M+7+J70bmNm8FMd0PHMjPLbI\ngKLCLxKN9whGzPcSjIaPCgvk0ZGxmS0FVgGrgeXAVqAG+G247StAKXCVu/8u7OZbZnYu0AoscfdN\nZjYKeBCoBrqAv3f3lWb2j8D54fb73f3BHrGcCTwU9n8IuI3gpOWfgKFm9oC7/2Wv3E540mFmZwAP\nAOXAYeCr7v5mmOtBgnXHRwB3ADcBUwlOKO40szzgOwTFPA9Y6u7fC1+3fwjjnAjUAteHrzFmtga4\nCHgMmByG8oC7//BE8YpkE33GLxKdO4HLzGzBCfZL8NEovga4BzBgJjDO3ecA/wHc2qNNnbufC/wz\nsDTc9n3gUXefASwCfmBmQ8PnCt19cs+iH1oGfM/dpwJ/AzwJ1AN3A8v7KPpXmNn6Hv+tOUZOPwLu\ncvfpBCc+P+nxXJW7TwuP8Vj4/DRgiZkNA5YAibDtLGCRmV0Ytp0N/BVB4R9L8HHKbQDuPhu4ACgN\nX5sF4WORWNGIXyQi7t5iZksIp/yTbLbb3TcAmNm7wMpw+w7gtB77/TA8xnNm9nhYMBcEzeyecJ98\nYDzBScXa3gcKTwrGu/szYV9rzayR4KQjh2OP7BMEJwQfm+rv0W8xwUnLY2bWvbnYzMrC9r/okdNG\nd28I2zUSzDwsAKaa2SXdbYEpBCckG939vXD/eqCs1+Frw9fgeeA54Ot9xSmSrTTiF4mQu78IvEhw\nQVy33p/RF/T495FeXXT00XVnr8ftBL/v8939HHc/h2C0Wxs+33aMPnL5eHHPIZheP971CCea6s8D\nWrvjCGOZ0+NixPYe+x4rv1zga73yWBoet2ceH7vWITzGZOB+ghOYN8xs+AniFckqKvwi0ftb4FJg\ndPi4ATjdzAaFo+CL+tlfDnADgJldBdS7eyvwvwTT4JjZZGADwZX2xyzU7t4MbA37wMzOB0YBG/tq\nc5ztvfvdbGbdMS4kuFYhqfZhHreaWb6ZlQCvAOedoE2nmeWZ2eeAZe6+Arid4HqCTyVxTJGsoal+\nkWgcHTH3mPJ/PnxcZ2YrgDpgG/ByjzZ9jbR7PpcAppjZeuAA8OVw+1eBh8xsA+HJgbsfNLPj9Xsj\n8KCZfZNgNH21u3ccp02C8DP+XtvvI5i6725zQ9jvXcCHwBf7yKP3MRIEFyhOANYT/A17xN1fNrO5\nx8ljOfAmwTUAXzCzujCfp9y9ro82IlkpJ5E4Wd8gEhERkYFOU/0iIiIxosIvIiISIyr8IiIiMaLC\nLyIiEiMq/CIiIjGiwi8iIhIjKvwiIiIxosIvIiISI/8Hu8Wc8ARUrJAAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10da06550>"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For large inputs, the pure-Python version is less than a factor of 2 slower than the optimized Fortran equivalent. Pretty good!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "  "
     ]
    }
   ],
   "metadata": {}
  }
 ]
}