{ "cells": [ { "cell_type": "markdown", "id": "3d0a47e3-7c72-43fc-8563-d44cf886a24b", "metadata": {}, "source": [ "---\n", "title: \"A5: Fast Fourier Transform\"\n", "subtitle: \"Assignment 5: Fast Fourier Transform\"\n", "date: 2026-04-06\n", "format:\n", " html:\n", " other-links:\n", " - text: This notebook\n", " href: A5.ipynb\n", "---\n", "\n", ":::{.callout-note}\n", "\n", "* Complete the following and submit to Canvas before April 15, 23:59\n", "* Late work will recieve 0%,\n", "* Each assignment is worth the same, \n", "* Please get in contact with the TA/instructor in plenty of time if you need help,\n", "* Before submitting your work, make sure to check everything runs as expected. Click **Kernel > Restart Kernel and Run All Cells**.\n", "* Feel free to add more cells to experiment or test your answers,\n", "* I encourage you to discuss the course material and assignment questions with your classmates. However, unless otherwise explicitly stated on the assignment, you must complete and write up your solutions on your own,\n", "* The use of GenAI is prohibited as outlined in the course syllabus. If I suspect you of cheating, you may be asked to complete a written or oral exam on the content of this assignment.\n", "\n", ":::" ] }, { "cell_type": "code", "execution_count": 2, "id": "1f4163eb", "metadata": {}, "outputs": [], "source": [ "using LinearAlgebra, Plots, LaTeXStrings" ] }, { "cell_type": "markdown", "id": "7f848e4a", "metadata": {}, "source": [ "::: {.box}\n", "\n", "***Name:***\n", "\n", "***Approx. time spent on this assignment:***\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "f6aaa7db", "metadata": {}, "source": [ "Recall that we are interested in solving the periodic problem \n", "\n", "\\begin{align}\n", " &-u'' + p u' + q u = f \\qquad \\text{on } (0,1) \\nonumber\\\\\n", " %\n", " &u(0) = u(1) \\nonumber\\\\\n", " &u'(0) = u'(1). \n", "\\end{align}\n", "\n", "We will take $p = 0$ and $q > 0$. We approximate this problem with $\\mathcal L_h U = F$ where $U_0 = U_n$, \n", "\n", "\\begin{align}\n", " (\\mathcal L_h U)_j := -\\frac{U_{j-1} - 2 U_j + U_{j+1}}{h^2}\n", " %\n", " + q U_j, \n", "\\end{align}\n", "\n", "and $F_j := f(x_j)$. By solving the linear system of equations $\\mathcal L_h U = F$, we find an approximation to the continuous problem. \n", "\n", "Recall that we may write $U_j = \\sum_{k=0}^{n-1} \\hat{U}_k (\\bm e_k)_j$ where $(\\bm e_k)_j := \\frac{1}{\\sqrt{n}} e^{\\frac{2\\pi i j k}{n}}$ is the Fourier basis (note that you will see different normalisations in different texts) and $\\hat{U}_k := \\sum_{j=0}^{n-1} \\overline{ (e_{\\bm k})_j } U_j$." ] }, { "cell_type": "markdown", "id": "b1a2dfdd", "metadata": {}, "source": [ "**Exercise 1.** Write $\\mathcal L_h U = F$ in terms of the Fourier coefficients $(\\hat{U}_k, \\hat{F}_k)$ and hence explain how one may solve the linear system $\\mathcal L_h U = F$ by implementing the discrete Fourier transform and the inverse discrete Fourier transform." ] }, { "cell_type": "markdown", "id": "b0a84771", "metadata": {}, "source": [ "::: {.box} \n", "\n", "***Your answer here***\n", "\n", "\n", "\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "cc187657", "metadata": {}, "source": [ "**Exercise 2.** Write down a matrix $\\mathscr F$ such that $\\hat{U} = \\mathscr F U$ and $U = \\overline{\\mathscr F}^T \\hat{U}$.\n", "\n", "(The conjugate transpose is normally written $\\mathscr F^\\star := \\overline{\\mathscr F}^T$)" ] }, { "cell_type": "markdown", "id": "a17d59a4", "metadata": {}, "source": [ "::: {.box} \n", "\n", "***Your answer here***\n", "\n", "\n", "\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "89419e2d", "metadata": {}, "source": [ "**Exercise 3.** Take $f(x) = x^2$ and $q = 1$. Use your previous two answers to fill in the gaps in the following code sample to approximate the solution to $-u''(x) + u(x) = x^2$ with periodic boundary conditions. " ] }, { "cell_type": "code", "execution_count": null, "id": "ce58bda9", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd0AT5/sA8DcJey/ZCihTqOICZQgqKsVFAXHPqjirtlXrrqtOnHXiqLgQcKKgIkuGIMoGQaYgU/ZKQsb9/rh+80uDYsDAQfJ8/iLHm7vnjuOee8e9R8IwDAEAAACiikx0AAAAAACRIBECAAAB2traiA4B/KtPJkIWi1VQUFBQUFBZWdn+t1VVVQUFBc3NzT0fWJcxGAw7O7sFCxYIfM0lJSV2dnabNm0S+Jrbi4uLs7OzO336tKBW+O7dOzs7u6NHjwpqhSKrra3Nx8eH//KXL19mMpndF8/3YDKZly5d4r/81atXGQxG98XTBffu3Tt16pS3t7eTk9Pr16+7YxNTpkz58ccfv1msurp64sSJGzZs6MIm6HR6AZePHz9+8apbUVHBXay4uJifP8fHjx8fP3587do1Pz+/pKSkbj8bsT6orKwMD15ZWbm2tpbnt/Pnz0cI3blzh5DYuoZGoyGEzMzMBL7m7OxshNCPP/4o8DW39/jxY4TQ+vXrO/vF8PDwFStWvHjxgmf5y5cvEULLly8XUIDfVldXt2LFisOHD/fYFrvGy8vLxcWFzWbjH6lUqr29/fbt279YmM1mr1ixIikpif/1p6amLl++nMViCSBWgWKz2atWrUpMTOT/K9nZ2YsWLWIymd0XVackJSX5+fnhP4eEhCgrK1dWVgp8KyoqKoqKipyPHz9+XLFixfnz53mKffr0CSFka2vbhU3Ex8e3TyhGRkbHjh1jMBicYlOmTOEpQyaT9fX1d+zY0dra2n61T58+HT58OM9XNDQ0duzY0dTU1IU4+SEm2LTaw+rq6ry9vffv3090IN+LQqF4eXlpaWkRHQgxMjIyLl26ZGxsPHHiRO7lurq6Xl5e9vb2PRZJa2vrpUuXbG1tN2/e3GMb7YKwsLDCwkI2m02hUBBCTU1N0dHRX7tr3rt37+jRo4cNG8b/+ocMGTJy5Mhz586tXbtWMBELyNGjR3/44YeRI0fy/xUTE5Px48cfP368Z9pFvik2NnbHjh0eHh4UCmXChAmNjY1xcXGurq6C3cqSJUvYbDbnY3V19aVLl9zc3FauXCnYDampqf38888IISqVmpmZGR4e/vvvvycnJ9+8eZO72IwZM0xNTfFiRUVFz549279//9u3b4ODg0kkEqfY/v37d+3ahRCaMmWKi4uLrq5uc3Nzenq6v7///v37paWlt23bJtj4/9VNCbZb4TVCJSUleXl5WVnZiooK7t/2xRph9+kTNUK8NfXYsWPdEVWnlJaWoq7eIPckBwcHXV1dzkc2my0lJeXp6dm+ZFZW1tChQzl1R/6x2Wxra+uSkpLOfrGxsbELm+NHfn6+ubl51+p2NjY2+fn5Ag+pCz5//hwUFMT5GSEUFRXV3Rt99+4dQsjNzY1n+ffXCM3NzbkXPn78GE9s8fHx+BK8RsipBOPS09NlZWURQq9fv+Ys9PPzQwjJyso+e/aMZ1ttbW3Hjx8/cuRIF+LkRx+uEcrLyy9dunTPnj379+8/c+ZMByVDQkKam5tnzpzJvbCmpiY8PFxPT8/KygpfkpSUlJ+fb29vr6CgcOfOndTUVElJSRcXl3HjxuEFoqKinj17Vl9fP2TIkMWLF0tLS/NsiMlkPnnyJC4urqGhQVNT08XFxdramrtAYmJiUVGRo6OjjIyMv79/amoqi8U6ffo0hmF3795VVFR0cXHhLs9ms58/fx4TE1NdXa2qqmpqajpt2jRlZWX8txiGJSUlxcTEFBUVUanUAQMGuLi4WFpadv5YIoRQTExMTExMaWkpmUxWV1e3tra2s7OTkpLiFGCxWE+fPo2JicH3bvLkyTY2Nh2v89WrV5WVlS4uLvhJj6PRaEFBQWpqaviBDQsLS0lJQQilpaUFBATgZezt7TU1NSsrK8PDww0NDUeNGsW92qKionv37hUUFIiJiVlYWMycOVNJSYm7wNOnT6lUqoeHR3l5+c2bNwsLC+Xl5adPn25ra9tBtNnZ2a9evUIIVVdXcyIxMDDg1D8wDHv27FlUVFRdXZ26urqTk5ODg0PHR6Cb9O/fn06ncz6SSCQdHR1dXd32Jfft27dx40bum24+kUikn3/+edeuXVevXu3UF/38/FxdXfv169fZLX7TgQMH1q1bh1eCO2vlypU7d+68deuWwKPqLDU1talTp+I/X7582dbW1s7OroPyjx49wjCMu8qYmZmZlZWlqqo6fvx4zsLXr19/+vTJyckJvz7cv38fwzB3d3eEUEpKCt7L8OnTJ86JbWZmZmFhwb2hjIyMgICAiooKHR2duXPnGhoadmHvpk2bZmdnFx0dHRERwXP142ZhYWFjYxMaGpqXlzd69GiEEIPBwKvsJ06cmDx5Mk95cXHxjRs3cg8vKi0tffToUWFhYXNzs5qaGt6epKmp2YWYEerLNcL+/fs3NTWpq6uLi4vn5eVxftu+Roj/RXnuUuPi4hBC8+fP5yxZtWoVQujKlSt4FZ7jzz//ZLFYS5cu5V5oaWnZ3NzMvcLMzEwTExP8t2Tyv6OQlixZwt1cvnjxYoTQtWvXBg0axFkVk8n8Yh9hTk7ODz/8gJfh/PMPGDCAU2DatGk8v0UIrVu3jnsl/NQI2Ww2ftDwE05cXBz/+eTJk5wypaWlnIZ7zuY8PT2pVCqnTPsaoZOTE0KooKCAe3P4n8/Ozg7/yJPkcMHBwdhX+ghPnjwpISHBHYaqqurLly+5ywwYMIBMJoeFhfEkyP3793dwHA4ePNg+kmXLluG/raqq4iR+zqanTp3aff0WHdi6dauHhwf3EkdHR29vb55iVVVVsrKyPCcq/+rr6+Xk5Orr6zv1LR8fH55GGoHAg6mrq+va11taWuTk5D5//izYqL5HRkbGyJEjq6qqOi6G32wVFxdzluBJUV5envvaYmpqKiEhwTkbufsIvzgWZseOHRhXjfDkyZPclxFJSUn8f7ADX6wRYhg2Z84chNAff/yBf/xijRDDMPz+kjMy4MmTJwghbW1t7p36msDAQE49REZGBv/Bysrqm1/8mj45apRDTk7ujz/+YDAYe/bsEdQ6f//9dxMTk4yMjKampsePHysqKu7bt2/VqlXh4eEhISENDQ3Z2dk2NjYpKSknTpzgfKuqqsrJySk/P3/Xrl2lpaUMBiM9PX3s2LHXrl1r34W5cePGwYMHx8bGfv78OSEh4Yt369XV1U5OTunp6cuWLSssLGQymbW1tS9evJg+fTqnjJGR0cWLF/Pz86lUKpVKDQ0NNTMzO3PmTGdvex8+fHjz5k1bW9vMzMy2tjYajVZcXHzx4kUjIyO8AIvF+umnn5KSkmbNmlVaWspkMtPS0kaMGOHv79+18WbcwsLCDh8+jBDat29f7f/gGbS9R48ebdiwQU5O7sGDB3Q6vaWl5ejRo3V1da6urgUFBdwlMQybNWvW2rVr8/PzGxoabt68KSMjs3v37tzc3K9FsmHDhqysLISQtbU1J5JTp07ha/P09IyLi5s+ffrHjx+ZTOb79+9tbW2fPHni5eX1nUeAR2VlZXR0NM/IOp7+v/79+/PU/9ovQQg9e/bMwsKCuzrO0draGhMTU1FRgX/Mzs6Oj4/n7lVCCCkqKg4cOBC/HelWbW1t0dHRVVVV3AvxoTqcj6GhoYMGDeK5s8FRqdTY2FjOGLrc3NzXr1+zWCzuMjIyMoMHD37+/Hk3hM8Lw7DMzMzk5GQ8/rq6uujo6JqaGu4yNTU1R48eDQkJ+WbVecKECQihsLAw/COLxYqKisJzXmJiIr6wrKwsOzvbyspKTk6u/RoOHDgQGRmJEJo6dSrnxObubMvMzNyzZ8+lS5cqKioqKyt3795Np9OXL1/ehbGaGIZlZGQghDqonFGp1EOHDr19+3bgwIGcNpXo6GiE0Pjx48XEvtFOiccmJiZ2//79tra2lpaW2tra58+ff1c/a5dTKIE4NUIMw+h0uoGBAZlMxk877LtrhEOGDOHuhNi6dStCiEKhpKWlcRbiTXk2NjacJRs3bkQI7dmzh3sT9fX1Ghoa8vLynFtyvEZobGzMc9fTvka4ZcsWhNCiRYs6dWQ+fPiAEHJ0dOQs4adGuH37doTQgwcPvlbg/v377cMuLS2VlJSkUChFRUX4kq7VCLGv9xG2rxHiVeRbt25xF1u9ejVPsQEDBiCEfv31V+5iv/76K0Lo7NmzXz8SX+0jDA0NxU85Go3GWVhdXS0vL08ikTIyMjpYZ6dER0cfOXLk9OnTkydP5iz89OmTvLx8eno6Z0lQUBDP4dq2bVtsbCzP2pYvX75x48b2W3n8+PHhw4fj4uKmTJly7969Y8eO3blzZ9euXRMnTuQp6eXl9csvv3RqFzpbI6yrq9u8eXNERISenh53vWfixIl//vkn5+Mvv/yyatWq9l8PCQk5ePBgXFzcjBkzbt++ferUqZs3bx44cMDe3p6n5Pr163tgBHJZWdmWLVueP39+8uRJDw+P169f792798mTJ4aGhpzzpKWlZdeuXfiYyaqqqvbjpbnFxsZyX6nevHmDEMIHMXFaOHx9fRFC3IeLZ9Rox32ECKGAgADu5WPGjEEIdTzSuH2NkMFg4ENdKBRKbm4uvhCvEaqrqw8cOHDgwIFaWlpkMplMJi9YsID7zz137lyE0O7duzvYIg6/W3V3d/9mSf714T5CnISExM6dO5cuXbpr1y78Qvydli5dyt1EMGLECISQlZUVp5USITRkyBAJCYni4mLOkoCAABKJhF+RORQVFadOnXrlypXExERHR0fO8hUrVnzzrufevXsIIT7HuTU2Nn769IlKpSKEVFVVU1NT+fkWh7a2NkLI399/woQJ8vLy7QsEBwcjhNauXcsdtra29syZM2/evBkSEiLwoWhfVFxcnJ6erqWlNWvWLO7lGzZsOHfu3NOnT3nK89TVbGxsjh8//vHjxy5sGj8CK1eulJSU5CxUVVWdP3/++fPng4ODzc3Nu7BaHnQ6/eHDh8eOHbtx40Z4eDiLxcJPxdDQ0NbWVj09PU5JExMTni7qUaNGGRgY8Kzw/fv3P/30E8/CtLS0oqIifFish4fHkiVLzp49O3v2bAsLi6amJp7COjo6MTEx379rHTh58uTmzZslJCRKS0sTExP79++PEGpoaAgPD583bx73vnD/E+Gys7MzMzP/+OMPhNCcOXOWL19+/PjxefPmjRo1Cr/Icje36OjoBAUFfS0MNpt98uRJfh5xMzQ0xPve2sMw7NChQ4cOHZKWlnZycvrtt99YLNb9+/e9vb3z8vIqKyvxwT5eXl5WVlY3btyg0Wjv3r3D4/8aKysrBQWFsLAwfHfwquGqVauCgoLCwsLwu1h8IXeXYaeoqqp6eHhwL7GxsXn9+vXHjx+/Odi4qKgIH+zd0tKSk5NTW1tLIpH279/P08UoJyeHd14qKytLSEiUlJSEhYVNmDBh0aJFeAH83Pti6wUPTU1NEokUHx+flZU1ePDgzuzoV/X5RIgQWrhw4bFjx4KCgmJjYzseDcEPnquJqqpq+4UkEklZWbm6uhr/2NDQ8OnTJ0lJSZ5EiBDKzMxECHFuu3CcJsevYTAYeXl5YmJixsbGHRQrKyvbunXro0ePGhoaeMLjuQR0bO7cuYcPH75z587jx48dHBzGjRs3ZcoUMzMzTgG8ObH9OYdf/TtobBQsvL5rYmLCM1zC0NBQSkqqrKyspaWF849EIpF4/mpqamoIIc5frVN65ggEBwfjwyhevHhhbW3N2c3IyMhhw4Zx36MYGRnxnEVfbBcqLy/nDK3iuHbtGt4WjRAqKytjMpmenp4IIT8/P0VFRZ7CKioqnObT7oC36quqqj5+/JjFYuHjJhBCMTExLBZr7NixnJJf3JcrV65wukVKS0vpdDreQeXr6ystLc3zL9DxvpDJZLzN4HuEhIRMmjQJv0epqKjA+xQQQsuWLRszZgzex/z8+fOKigrOXTuJROp4WIqYmJidnV1wcPD79+8HDx4cFhampaVlZmY2bty427dv4+d8ZGSkjIwMZ9xfZw0cOJBnCf//LHguRwhJSEhoaGj8+OOPK1as4P7D4f766y/u+9eMjAxXV9fFixfLysriORg/99rfirWnrKzs5eV14cIFCwsLS0vL8ePHOzs7jx07Fh860DXCkAgpFMq+ffvc3d3/+OMPvKH5e3zxaLZfyP0/hs+nQCKReLqpEELS0tIjRozgqWa1v9zwaG1txYtxxq2019TUZGtrW1RU5OLiMnnyZHV1dWVlZQqFsnDhwvLy8k4lQiUlpZSUlGPHjj18+DA4ODg4OHjTpk2TJk26du0aXlnEW27bX4ZUVFQQQnhNtAfg4yTb9xKRSCRFRUUajdba2sqdCDs4ep3VtSPg5+eHtyp/044dO4YMGeLk5CQnJ0elUh89erR7927Ob6OionjGPPOppaWl/cm2f/9+zvmclJQ0ZswY/CPPGEKckpJSB5M0MRgMb29vnjMtMTHx48ePPJ1VAwcO/OIuSEpK7ty5EyHk5+dnbW2Nn28IocjISH19fe5bmS/uy86dOzljJZKSkkaOHImfANy3cXzui0A4ODhwzkA8PeC1WEVFRc5IqylTprR/xrxjEyZMCA4ODgsLGzRoUGxsrIeHB4lEmjBhwj///BMbG6uvr//x40dnZ2fu5opO+VoKwfh4JYOpqSneKdgpFhYWBw8e9PT0/Ouvv/BEiP+tc3Jy+Pn6uXPnbG1tfX19o6KikpOTvb29dXR0fHx8+JlM54uEIREihNzc3EaPHh0TExMSEtL+t3iDHpPJ5L4y1tfXC2rrSkpKJBJJTEwsISGha2O7ecjLy0tISNTU1DQ3N3+x9xshdOfOnaKiouXLl3NPN8Vms7u2X8rKygcOHDhw4EB5eXloaOjp06dfvHixYsUKfCgXngC4B47i8Jpu+/TAgR8Nni73Lh95PAVyxkRwtLW1ff78mUKhfPMOo8s4R4BnecdHYNKkSd+s/SOESCQSPlAZv2F68uRJc3Mzp6mquLi4qKiI8wxPp4iJiXE/ZYHjXKkxDIuKilq/fn0Ha2htbe3gRltcXLx9s97ly5enTZumoaHBT4RkMllOTq61tTUoKGjv3r2c5ZGRkfggkY73RUFBgfsreB/813S8LwLB3bIXERExaNAgvKX3O3HGy5ibm1OpVPzjhAkT8JZSfX19Tpk+BB9j//79e/yu3cHBYd++feHh4TQajfuprS8ikUjz58+fP38+lUqNiYnx8/O7du3azJkzc3JydHR0uhCMkCRChNChQ4ccHR3/+OOP9je22tra2dnZpaWl+BmDwwe8CISsrKyFhUV6evrbt287eHSGf2Qy2crKKiYmJjY2tv0jNbi8vDzU7uxPTk7+zvqZlpbWwoUL3d3d1dXVw8PD2Ww2mUweNmzYs2fPoqOjOQ9s4KKiohBC7edD4sBPytLSUu580P7I4+f9N7tnhgwZIiYmlpGRUVNTgzdZ4169esVms4cOHfr9lzk8kvazIVtaWgYEBERHR+Nd+hz4EfhaP4qKigpeZeyUiIgICwsLTo9gRESEuLi4vb09g8FobW3tVLJXVFSsra392m9TU1Orq6s5rVgsFqu5uZln/XV1dR3c6AhKQkJCc3Mz5yHaxsbGlJQUPENXV1fjzXSKiop1dXVfW0NOTk5paSlnEiI2m93Y2MjTeFBfX9/BvrDZ7FOnTvEzEbaxsXH7ntf2IiMjuVsIa2tru3Ay4IYMGaKhoREZGYknD7wvUEtLy9TUlM9E+LUTm0B4+5mkpCTeojBu3DgTE5OcnJzTp09/bV6nhoYGnvNTWlp64sSJEydOpNFot2/fjouL61rbSd9+fIKbg4PDxIkT09LSnj17xvMr/CqMT1uAKysr43744fvhvYMbN25s3/aC92x11ooVKxBCmzZt4qk/4W106H85Jjk5mfMrJpPZtYnB8vPzecaas1gsNpstJiaGn6Nz584lkUgXL17k7ux88eLFq1evVFVVO2iOwI/83bt3OUsaGxv37dvHUwwf9//NnjZ5efkZM2a0trYeOnSIs5Dz8IxApixXUVGRkZEpKiriycpz5syhUCjXr1/Pz8/nLIyLiwsODpaTk+Pnssi/Dx8+cI/MioqKGjp0qJycXGhoaPu2944ZGhryJA82m/306VM8O4aFhZFIJM4MDL6+voWFhTxrqKmp4X7mtZt8+PBBWlqa01UWExPDZDLt7OwYDMaVK1fwhYaGhjxJHcOwkJAQvB8LHy3CuSe7e/cuPl6aW8f7QiaTN27cuIUPHfy509PTk5KSEEKfP39OTU3lxNPY2Pg9FxwSieTo6NjQ0ODj42NiYsKpZTo5OSUnJ798+VJVVXXo0KEdrIHPf7EeU1JSgjf+Ozs740vIZPLff/9NJpN37Nhx+fJlnlbZ+vr6DRs2nD17FiFUXV3d/pYI/4ft8q2w8NQIEUIHDx58+fJl+1vgpUuX+vj47Ny5s6ysbPjw4fn5+RcuXDA2NsbnNxKIFStWPH/+/OHDhz/88MOCBQtMTU0ZDEZhYWFoaGh6enpjY2NnVzh//vxHjx7du3fP0tJy5cqVRkZGtbW17969i4qKev/+PULI2dlZXFz8+PHjbW1tDg4ONTU1Pj4+5eXlGhoaX3wpRwfwFolZs2aZmJhoaGiUlpZevHiRRqOtWrUKT4QWFhYbN248fvz4mDFjfv311wEDBiQlJZ08eZJEIp0+fZrTSdPe3Llz9+3bd+HCBSqV6uDgUFpaeuXKFU4/EIeVlZWsrKyvr29NTc3AgQMlJCSWLFnCmZ2A29GjRyMiIo4dO1ZaWurq6trS0uLj4/P69ethw4YJalZMR0fH4ODgMWPGWFtby8vLjxo1yt3d3cDAYMeOHXv27LG1tf3tt98GDhyYnp5+/PhxDMO8vb2/+HBbl2lqanKahgoKCoKCgvCHR9+8ecPdccgPS0vLtLQ07iW3b99esGDBuXPnli9fHhQUpKqqircuNjQ05OXlLVmyhGcNqampX2uTECANDQ1JSUl8GgoWi3X27Fl5eXl9ff2QkBBOpcrS0hJ/kIDj/v37Hh4e3t7eGzZsePTokZycHP5MXnNzc0pKCj5qhmdfvjh7g6DQaDQbG5vBgwcnJCT4+voqKChwmqD+/vvv7xxZPX78+Lt379bV1XG3SUyYMOHMmTMNDQ0eHh6cSTy+SEFBYcSIEe/evbO3tx82bJiMjMz48eMnTZr0PSF1ysmTJ/GR8AwGo7S0NC0tjU6n9+/f/8iRI5wyTk5Oly5dWrVq1fLly8+cOYPPNdra2pqSkhIcHFxfX4+/iCYpKWnGjBk//fSTlZXVwIEDm5ubnz9/HhAQoK+vzzNZcScI8FGMHsP9HCEPTs8Kz1yjV65c4YxYkZCQ2LhxIz6ZVvvnCENCQri/GBERgRBavHgxz4bwqxX3EiaTefDgQZ7HY/X19Tdv3swpg/dhREZG8qyNTqeTSCSeORoYDMaff/7J3RQgKSm5YMECTgF/f3/uFsLRo0fn5OTglTDOewNycnJIJNKUKVO+ciwxDMPOnj3Lc6esoKCwdetW7qcG2Wz2wYMHuYMZMGBAYGAg93qCgoJIJBLPg2tPnjxRV1fHv0Imk+fPn4/flnI/R4hhWFhY2Lhx47S1tfH/Z3xWC7zK4uXlxV0yMzOTexpuCoUyb948npeQ4DPL8Owm/kzx0qVLOzgUGIaVlZXNnz8fH4mKuGaWwTDs5MmT3K1b2traN27c6HhtXZCenm5ubu7n5+fj47Nv3z5/f/9Ro0b5+Ph0YVuxsbHa2trcSxITE4cNG3b//v3ff/89KyvL3d392LFjfn5+27dvb2ho4Pk6k8lUUFDgfn6RH12YWYbBYEyaNGnHjh3379/ftGnTmzdv9PT0bty4sWXLFs6Z/O7dO1VVVe6ngVNTU4cOHRoYGLh58+a0tLS5c+f+9ddf/v7+27Zta/9SGjab3a9fvzdv3nQqsE5hsViOjo6HDx8+e/bspUuX8KnmgoODd+3aFRMT850rx7tC0P/mTsPV19fj3fDtXyuhpqamrKzMswZ3d3d9fX282sQzswzP1/FZli5dutRBSF+bWYZH+5FBCgoKw4YN27FjR/s/E4ZheJ7j6SYcMWLE1atX8ZMhIyPDxsaG+1EuMpns7OzM87xyp5AwPsYFCYeWlhb8ofvBgwdz5w/BYrPZ2dnZ5eXlcnJyurq6Xeu55YZPUlNbW6uurm5gYMAzALWtrS0nJ6empmbAgAHtx0B3SkVFRUlJSXNzs7a2toGBwRcbGdra2tLS0urr67W1tc3MzPgcmNrW1pacnEyj0YyMjNpXB7umpKQkPz9fTEzM3Ny8BzqxOPBZdWprazU0NMzNzTu+De+ytra2rKysfv364edPeXl5S0tLF+Z+xDDM2Nj40aNH3A9+NDY2FhUVGRsb49eaDx8+yMrKfvFEjYuLW7NmDXfzOz86NViGW25uLoPBMDU1JZPJVCr1w4cP5ubm3Nc7CwsLX19f7j7ppqamwsJCIyMj/ImFvLw8KSmpL865mpycPHfu3KysrC5Mu8o//N9fRUUFn1elpqamoqLC1NRUIGPoRAp+Anz+/FlaWtrExATvJ+bW3NxcVFRUVVWlrKysr6//nRcBEUqEAIigS5cuJSUlXbhwoQvfnTVrlpubG8/0Bd8UFRU1atSoDhrMu8zX1zc8PPyff/7pwncXL17s6OjY8bBSILKEZ7AMAKC9ZcuW5ebmtn/m5Jtyc3MbGxs7mwURQg4ODt2RBRFCCxYsqKio6MLcQMXFxSUlJZxJTADgAYkQAGFGJpPPnDmzevXqTk2gTKfTN2zY8Pfff3dfYF1AIpHOnDmzdu1afiZC42AymWvXrj179my3NoqCPg2aRgEQflFRUZGRkfwPOt22bZurq2uXp+zqVvHx8U+ePGn/Upev2bNnz4QJEzp+5x8QcZAIARsFPrwAACAASURBVBAJPLMQCLBwz6uqquIMRf6mXr4voDeARAgAAECkQR8hAAAAkQaJEAAAgEiDRAgAAECkQSIEAAAg0np7IkxNTX348CH/5Tv1sBQQIDjyhMAwjOfNIaBnwJEnSncc+d6eCBMTE58+fcpnYQzDOG8pAj2JzWa3f2kq6AFw5InCYrF61ev9RAeLxerUjAr86O2JEAAAAOhWkAgBAACINEiEAAAARBokQgAAACKtE4mwpaUlJiamoKDgawWKiopiYmIyMzPZbDbPrzAMKygoqK+v516YkZERHx8v8G5PAAAAgH/8JsJ3794ZGhru3Llz7Nix69ata19gy5Yt48aN27lz5/Tp04cMGfLp0yfu3168eNHQ0NDHxwf/yGAwpkyZ4ubm9uuvvw4ePLi0tPQ7dwMAAADoGn4T4ZYtW3755ZeIiIjk5GQ/P793797xFNi1a1dhYWFEREReXp6BgYG3tzfnV8XFxWfPnh07dixnSWBg4MePH9PS0uLi4saOHXvgwIHv3xMAAACgC8T4KVRTUxMeHn7z5k2EUL9+/VxcXAICAkaMGMFdRlZWFv+BRCKpq6tzvwPTy8vr4MGDFy5c4CwJCAiYPXu2lJQUQmjx4sVubm7nzp37/p0BAAChxGKxMjIyJCUlAx890dXWuhv0YtQQs+dRcRPtR79OSp87Y1JWXtGy+bM/f/5sa2tLoVCIjreP4SsRfvr0SUJCQlNTE/9oYGCQl5fXvlh6evr169c/fPjAYDAOHjyIL7x69Wq/fv2mTp3KnQhLSkrc3Nzwn/X19aurq1tbW2VkZNqvk0qllpaWhoaGcpaMGTPmiyURQhiGsdns9j2UoLux/4foQEQOHHmi9MCRb2xsfPnyZXl1nff5q5i4ZLW0DiPvDWPyZtLRbdiii89OeqGFF96cXonmnorYvAGb/PtpRxeKsZ1a8042g7r7l2Xi4pKzPD2kpaW7L0JCdPbIk8nfbvjkKxG2trZKSkpyPkpKSra0tLQvJicnZ2RkRCKR7t69m5WVpa6uXl5efvDgwdjYWJ6SVCqVs0K8Xvi1RFhZWZmamspJq+Li4vv37zczM/tinBiGUalUfnYbCBabzabRaNzNAKBnsFgsmFmGEEwmk8FgdMf7XGtqalLSMzbvO9rKJJVrj2a/CcDm/02+uUbcwkROUZkdc26wrW15+L6xs9zT3p20nOURl3BKbfTonMRrbHkFWnN1RdknluPKVdsPkobN2HPOV5rM9j1zuH///vLy8gIPlRBMJpPJZPKfCGVkZL6ZFPhKhFpaWk1NTW1tbRISEgih6upqLS2t9sUMDAy8vLwQQrq6ujt27IiJiTly5Iient61a9cQQvn5+RiGDRo0yM3NTVNTs6amBv9WdXW1hITE114hra+v7+Liwhll0zEMw0gkkpycHD+FgQCx2WwKhcJpHgc9hsViiYuLf62NBHQfPBEKsL5Fo9GeP3/+KjHV52EYrbGG4XZI7O5GaUaC8WBT3ap7a275yEhJjBp1Wlxc/GuXdTqdXlxcnJOX/yIiJnSQwaeCV6U0Gst2ib3rfGkFpUv7NpuZmvzwww+CCpgoeCLEa1CCwlci1NXV1dTUjImJGT9+PEIoJiZm1apVHZSXkJDAb5QmT57MaVDlNmrUqOjo6JUrVyKEoqOjraysoDIBABBNtbW15y5fi07KCK+VZ6UGk4a7KuW+MC+4ufrvo+PtbdTV1flcj6SkpJGRkZGR0dQfnRFCVCr10j83X0QnRElINCkZzlm3VUrT4Nhqz1HDLUeOHNmdO9T38JUIxcTE1q9fv27dOryds6yszNPTEyH05s0bZ2fn2tpahND69evNzc3V1dU/fPhw+PDhEydOIIScnZ2dnZ3xlURHRzs4OOBdg15eXpaWlidOnNDR0dm1a9f58+e7a/8AAKC3+vz586GTf8dn5L0WH4xevZJW729oMnDXsjFjxvyqra39nSuXlpZev2r5+lXLMzIyMjIz12xLb6K1rv3TW0pD7/LWn+1txujq6gpkL4QAX4kQIbR582YVFZVbt26pq6tHR0fjTTEaGhorVqzAC4wbN+7Zs2c1NTWampoPHjzgflgC5+HhMXDgQPxnfX39iIiIc+fOJSYmXrhwwdXVVUC7AwAAfUBDQ8PJc5ci3ma8Ig8mxcfLK7y3GWd37dRhDQ0NgTePWVhYWFhYTJ40KT8/33nO0sbG6gUbdsoqKj2/dnL06NGC3VYfReqOzl4Bunz5ckJCAv99hC0tLdBH2PPYbDaVSoU+wp6HD5aBPsKe1+U+QgaD8eDBw/thsf4VCij+lpKyytiRQx/c8OmZ7qHGxsbc3NwJs5c1U+QlGc0qMmIxTwL09PR6YNOCQlgfIQAAgO9XVFR0I+DR3qfp7OxIeXnFkaNHhT3y68kAFBQURowYkfzifsLbpKX7L5ZKq3pu3DPaWPfYvp3i4uI9GUmvAokQAAB6wqPHQfO3HW39XCohr9JfRzM78RVRz3oZGBgYGBgwmcwX8Sm3EwqSaz6NfvjEZdJ4RUVFQuIhHCRCAADoXtnZ2T/OWVZDY1Kl+inJ1WXFPunXrx/hTzzPnzPLddqUPPf5ORWVC377U1Fqz9vnDwwMDIiNihDw7DkAAHQXOp2+77D3mr0nizTHMBHl3t6Vaa+eaWhoEJ4FcXJycnHPH4ZdPyktr1iP5NYf8Tnvc5XooAgANUIAAOgWra2tt/3v73mRj0rSXR1Hu8xaOWPaFKKD+gJLS8uI68fvRSQcfhD3MiHNeoTl8OHDiQ6qR0EiBAAAwUtLSxvnsaixtU1aRtZIW/XO6QOCHegoWCNHjtTR0bl75045Fbms+dPKQC3wnwv4VGKiABIhAAAIEpPJ3H/0RGTGxzp5PTlyRdozP319faKD+jYtLa2CpJg7gQ8WXH71/ENBelb2CMshRAfVQ3pFOzUAAAiN4Ocv/grJiskq2j5rwot/TvaJLMgxzXniT9p0TW0dh7lrRk+aISJTukONEAAABKO8vHyMs1tlC1NSQtKwn9yW1Uv63PwecnJyAVfP3fa7u/AyOansfV5RibmJIdFBdTtIhAAAIABZWVn3wuNLJHUlxGgvz++0trIiOqKuc3Odkfo+L7hM29pt6ShDrdB7t8TEhDlZCPO+AQBAz8jKyhrjvrSZRZk62Wna8EF9OgsihKSkpA7v2W587Z8VpaPjPsSXlVcM6C/MM3RDIgQAgO9y5oJPcMrH5laaggTpyNq5JiYmREckGLM9Z77LPPhca/rYn7d6jjY+sncn0RF1FxgsAwAAXZecnLzF59GL9GLvPVs/xL0QmiyIEJKVlT13bP8qS9lihcFn/UOqqqqIjqi7QCIEAICuYDKZC1dtcP3tMEZr1qcVek6y79evH9FBCd4s16nDG16TDUcf9o9ITU0lOpxuAYkQAAC64n3OB/83+Z8k+1/YuTb/XfT3v0q3d+rfv//b0Mdrx5mc8H/h6LagtbWV6IgEDxIhAAB0DoPBmOwxf9LPmxRV1MZJl81wdiI6om43fqihQmVaC0WugS1JdCyCB4kQAAA6J7+oOCa/pkl7+OGFTi/v31JSUiI6om430WnCp3cRczcfMrAc7fTTHKLDETAYNQoAAPxiMpnjXWcnFVT0t3SwU6bOdHMlOqKeIycnp18VTzf/MTb5AZ1Ol5QUnqoh1AgBAIBfZRVVCXkVVL0xv48beP7YAVlZWaIj6lEb13htHCqlMOvAX3fDi4uLiQ5HYCARAgDAt7HZ7Mnu88wneRrYTvllmOwcTw+iIyKAoqLi8b3b3KXz9vn4j3KaymaziY5IMCARAgDAt9XU1kWl5rYYjlsxVPHEX3+KWl2Qm7WhlkxtfhNJFiMJSQYRkt0AAIDus2Ttb2aTZmvYeywzYv+8aB7R4RBs0fw5H17eHbPz9uhFWx4/DSE6HAGARAgAAB1paWnxCw6vMXZebEK55H1AUVGR6IiIp62lpZNw4S1LZ9Ha34mORQAgEQIAwFedunB5zLz1cqPdXSjZqxaLel2Q22I3l34pt9hmE6hMokP5bpAIAQDgyxgMxvZDJ9Llhs4eiJ7e8tHU1CQ6ol5k/DiHqsyEsbNXDp6y8NotP6LD+S6QCAEA4AsSExPXHb7ENptkmPtg5ZyfiA6nl5J6daFI3/nXHfuIDuS7QCIEAABeGIZN8lx8Mb7M3Vw5NyHc3Nyc6Ih6qY0/z9V9e4liu4DOIjqU7wCJEAAA/oNGo0UkZ7eqGinmhHg4jiA6nF7NZszokneRoya7jZy97umzF0SH00UwxRoAAPyHxZhxRZiqi/Ok+/tXi4nBRfLbKM+9M5TsFq/7/XNuGtGxdAXUCAEA4P8x2KisiSEmIz9MgQZZkE+r5rurxp2VtHLDiI6kayARAgDAv7w2/qFqZm06Z0vokXW7/xCGJ+R6xhTnSVUZ8drTflnpfbOkpITocDoNEiEAAPzrzsOnTWOWGVW8srezIZPh8tgJZBJSCtp2KTLbccZsomPpNPhLAwAAepOYOOf3A5SpWz2ls4/u3ER0OH2SwzAz6bxIlooe0YF0GiRCAABA0xes8KtQGFYddfe894ABA4gOp0/a/vv6zKin1CXX06v72GQzkAgBAKKuuKqeOmisYsK1pTOciI6lbzNQVxydeXGklfWi1RuJjqUTYEwUAECkTXKbG5VVMsxpevzTJKJjEQb0DwltoxeERt8lOpBOgBohAECkvfvwkWk0tn9zPtGBCInrZ4/9bISkV1xj9p239kIiBACIqLCIKLvZq8Rd/zw30+LqmaNEhyMkNDQ0Lu/Z0F9Zeua2k/n5feP2AhIhAEBELVz3e2y/yTqJF70WzpGXlyc6HKHSfOPXh5/EpsxZSnQgfIFECAAQRaWf65nWC9TCDv6+HN4yKHjTJ9hLvb4qbzaG6ED4AoNlAAAix3PpqsevMw1H2mdkJRAdi3DatXnDvBVrrYOwKipSlyY6mm+BGiEAQOTEpecxTCeoN+YRHYgwG6QkNroixNzW6U7gA6Jj+QZIhAAAEZKTk7Ng0z7q1D9PTDe6988FosMRcoUPz1RP3L7joDfRgXwDJEIAgAj5afGqm/X6Ss/2/LJkrrKyMtHhCLkTe7fqvj5jNHc70YF8AyRCAICowBBiWfwoG+a92M2F6FhEwqQJ49Nf3k/WdCpo6tUvaIJECAAQCacuXlE2HkklSdblpezctIHocESFkgQaknpxyPBRpy5eITqWr4JECAAQCedvBDRM2iaRFCAOl72eVRwZ0OK849KN3jvpGpwRAAAhx2KxIuPfNU3b79Qa5+9zmuhwRM7tCycta16Zrei9c/fAc4QAACE3d/nah7kt/bCG0JhHRMciikYMHxZ2bZixP6O4GRsgRyI6nC+AGiEAQMjlNWJsRBogS3QcIkxFEg1Ju2QxcvT1O/5Ex/IFkAgBAEKruLh46ea9hZY/3927JuJxb7wEi47CFzebnHecvnyT6EC+ABIhAEBoLd2w9dpnbbHba9zGWUlKShIdjki7/vcxw8Jgy5/3EB3IF0AiBAAIrTZTJ6lX5zynTiY6EIDG2tuG3z73iGVR30Z0KO1AIgQACKGIyFcG1k4ZNYyyrKS/j+wjOhyAEEL9ZUlDKsNt3Jcmp6QQHct/QCIEAAih/X9fKRr9CyX8rDI0iPYm769uez/op192HCA6kP+ARAgAEDZsDKM7rtFLu3F83w6iYwH/8euq5TLPDg6dspDoQP4DniMEAAiV67fvrtt1SGyQdWUYTCLT62xat9LSzevXeBaGUO95ohBOEwCAUAkMjW4aMU/iYwJkwd7JSYeE6M0nH8awWCyiY/kXnCkAAOHxsaS0xGHr9P5YyJ2rRMcCvoyEUMPxn7Ycv+K1cQvRsfwLmkYBAELi2Olzu3zuiZOxlNRwomMBHeknJ1na0tSKtIgO5F+QCAEAQiK5qIqmOkj1cxLRgYBveP3i8YaH79la5kQH8i9oGgUA9HlsNvtlVEy82bI/Frkmhj0hOhzwDRISEn9OG+pfhNXQiQ4FIQQ1QgCAEDhw9MT+R2/Fq3L+yoPqYN+gLo1Ms/2sfnwceHz7MEtLYoOBGiEAoM8rZ8kwm+s15CSIDgR0woe7RwvM5+7xPkt0IJAIAQB9WVtb2w3/e4+kHS+eOZkZH0l0OKATtm5YI/PqnI3bYqIDgaZRAEBfduj46b0v8iTydy8tTCOT4c6+L/ltrZe68/LbeezNREcC5w0AoA9rVdJDxSn6mqqQBfsiTwNyUg2W14gRGwacOgCAPolGo53x+ed67YCAoJCM1xFEhwO6QpKClPxWDR023P/eQwLDgEQIAOiTTp69sNH/Xf2lZS4DZaA62HfV5SS2Dp8dHB5NYAxw9gAA+iSm5mCUF2ekoyEhAYNF+7Cgmz6GCsh6MZHTrcFgGQBAH0Oj0c5cvv73Z+PHEXEu+vC+wb7N2srqsPrIkxnsVcTFADVCAEAfc/Hy1a2PMht8N45VYxAdCxCA6QPI+Y0oqaKVqAAgEQIA+pr+Q1BejLFOP2lpaaJDAQIgRkakv13t7B1v3PEnJABIhACAPqOtre3k+cuH39EeRySkxoRSKBSiIwKC0Vb1kWro8C49k5CtQyIEAPQZN27e+t0/sfbWFmv5RqJjAYIUFeRvPths+LxNhGwdEiEAoM8g6VqgoreD1BUUFBSIjgUIkpmZ2c41i28VSxGydUiEAIA+gMViXbt5Z198853Q15mvI8TFxYmOCAjYT/rklFrsYzMBs8xAIgQA9AEBgYFeF0LKfLfYy9USHQvoFhJkJHV1oYXliBcvw3p405AIAQB9gYYh61OGjhxFUVGR6FBAt2Cz2Y0FaS0/zAiNiu3hTUMiBAD0ahiGhTx7vvct88LTuPzkOHhkQliRyeR7V8+qy0mOnb+uhzcNM8sAAHq1kJAQ910XWfXlk2fdI5H6Ex0O6Ebjxzn+0W/svQpsmkmPbhdqhACAXo2toNlWXaImxlCQlyc6FtDt5g4iP/rIbu7ZKYMgEQIAeq/U1NTDObJ77kQUpSUoKSkRHQ7odurSSCvhgp374pKSkh7bKCRCAEAvFRsbazd/ffyBuXM06+EVEyKira2tJPhSqsLIi9du9NhGIRECAHopOlmqtaVVgcyQloSnBkWFhITEPLfpYimPnKb91GMbhcEyAIDeqLq6+nKN3vLTgX/ZyKmoqBAdDug5F7wPoJ9Ybygkx57aItQIAQC9TlZWluGYif5rJmwYjCALiqB5huRbeewe21wnaoQVFRUPHz4kk8murq7q6uo8v21tbX316lVubq64uLidnZ2FhQW+PDc3NyEhoaqqSkdHZ9q0aTIyMnjhp0+fcr77ww8/mJqafve+AACERF1TS7OYgiyFRmJQiY4FEMBOk1Tf2BiU+HnaqJ54kILfGmFhYeEPP/zw7t27169fDx06tLS0lKfAq1evTpw4UVBQkJKSYm9vf+7cOXz57t27X7x4UVlZefHiRQsLi8+fPyOEqqur586dG/A/2dnZAtwlAECfxmKxQshDx2/wjg28YmLSsw+Ugd6B0dZWt8fB8+c1d/wDemBz/NYIT5w44erq6uPjgxCaP3/+mTNnDh06xF3A2dnZ2dkZ/9nOzm7//v2rV69GCN2+fRtfiGGYpaVlUFDQ0qVLEUJSUlL+/sS8gxEA0GuVlZWNcJr2mYrevgyyGKRNdDiAGBiGSVJQC0WKRmvrgc3xWyN89uzZ9OnT8Z+nT5/+7NmzDgp/+vSpf3/eCSCam5sbGhp0dXXxj2w2+/79+0FBQXgdEQAAEELFJSXV4hqSyhr06p57jAz0NpKSklmvwwfN2z1o0twe2By/NcKysjItLS38Zy0trbKysvZl6urqJk6c2NzcLC4uzp0pT58+7evrW1BQ8Ntvv02aNAlfaGho+OjRo7KysgULFty8eXPq1Klf3G51dXVKSsqBAwc4S+bNm8eJhAeGYXQ6Hd7P0vPYbDadThcTg0HIPY3FYtHpdGF6UXuy7DC9H3/eO4xtaWlJp9OJDuermEwmg8Egk2G8YXdRUFBYOG7orQ8Ma+X/vJiJyWQymUwSicTneiQkJL5ZmN8rF4lEwjCM+2P7MvLy8hcvXqytrT148OC2bduuX7+OL3d3dx89evS7d+927949ceJEKyurAQMGpKam4r+9ePGil5dX+05HHIZhbW1t9fX1/4YrJsYdBgBAaDQ1NU30WJhR0/bEz9fRUJXocADxZhtgNsEk71FIvJvvN/hNhJqampWVlfjPFRUVX6yTiYmJjRgxAiFkbGysr69/+vRp/IUpOjo6Ojo6VlZWGRkZV65csbKy4v7W1KlTV65cWVdXp6ys3H6d/fr1s7KyOnr0KD9BYhjGZDIlJSX53CkgKGw2m81mw5HveSwWCyEkHEc+JSUlq1WarKJBKs2QNJ9EdDjfQKFQyGSycBz5XstQEiml3VmSlHdt70Y5OTl8IYVCoVAogj3y/ObZSZMmcR54ePr06cSJE/Gf8/PzaTQaQoi7olZYWCgjIyMrK4thGPfyoqIiNTU1nsKxsbEqKiowiyAAIq5Fd6S0kdUqO4Px48cTHQvoFWpqasoenriX0/LPjZvduiF+a4QbN260trYmk8ksFiskJCQxMRFfbmZmFhoa6uDgsGXLlqKiIiMjo+rq6sDAwP3794uJiVGp1OHDh0+cOFFOTi4+Pj4/P//y5csIIW9v79jYWFNT04qKisDAwDNnzvDf4AsAEDIsFmv1ph1+GXU+p456msErJsC/lJSUhhgOSPgQNcZ2frduiN8aoaGhYWpqqqmpqYWFRWpq6oABA/Dlfn5+5ubmCKFNmza5u7srKSkNHz48JiZm48aNCCEpKanr16+bmJioqqquXr06JycHb1NduHChp6ensrKyjY1NUlLS4sWLu2XnAAB9QUZGxo2Y7FYkKfkhnOhYQC9CoVDiggPHnowuVRzcrRsi9fKxJ5cvX05ISMCfX/wmDMNaWlo4Tcmgx7DZbCqVKisrS3QgIgcfNYpP2NR3ZVfTh7p5Wco2hdzy6SsTquGjRqWlpYkORPidf8+OrcRuOv47NBofNSolJSXATcB4dwAAke49eLg7smL/6UubLOFFS+ALPAzI298yqEyKdLflK3gIBgBAmIKCgsU7jr9PStDLCyI6FtBLqUlicsF77d0Xdd/sK5AIAQCEkVXRoLNJatXplkN+IDoW0Et9+vSpNiMmiWzw+HF33S1BIgQAEKOkpOTXx9kLL4VXZCUaGxsTHQ7opXR0dKbaWFKKk8dP+rGbNgF9hAAAAtTU1FiOm1IvrXVz33oSyYXocEDvRSaT/S6eGPeUmc4mG3TTJrpntQAA0BGKmFgzA8mwWjQUBDn8DwgrdwPyvaLuesYBaoQAAAL4lUhaHAwPGNM80ECf6FhAH+BhQN79jtHGpnRH7Q1qhACAHtXW1mY8wm7NVJsNqrmQBQGfNKQw6ecHJy1Y09TUJPCVQyIEAPSo5ubmkvpWst4wRlkO0bGAPiMvL68uIzamTpYz67UAQSIEAPSofJay3IKTpxaPWzB/HtGxgD7DwMDAxkSbXJY51sFR4CuHPkIAQM9Zuu53/5jMvd6nV48fS3QsoC8RExMLvXtt2ANmIRmpI6ZgVw41QgBAD2loaPB/FtliPJ6ZIvjWLSAKftInPywW/NhRSIQAgB5Cl1DErGY7sLMWzp1NdCygT3LVIz34KPjVClUifPv27Y1bt/EXBQMAepXgZ88nrT+4bNmyyIBrmpqaRIcD+qQhKiRGZthx3weCfW+S8CRCOp0+edaSdYGpp85dJDoWAMB/0On0Oat/Ty1vUXh9mehYQB+Wm5tbG3hg75V7z58/F+BqhScRSkhIaGmqs/MSBpt17yscAQCdJSYhyVY3Ucl9NtnBhuhYQB+mpqbWj0JTbCgwMBDkbGvCM2qURCJlxIbZBDZSLBSIjgUA8P/q6+v/epk/fIdf5FQxEtHBgD5NWVk5LymWTqfLy8sLcLXCUyPEuQ6UeFDEJjoKAMC/MAwzH+3o/dde25zLkAXB9yOTyeLi4gJep2BXR7hpuqzHxWxWd03NCgDotAYGSQKxtGUpRAcCwJcJT9MoboAs0qSVXXhaunrKGBIJbkABINjrKqSwPfzesE9jhsGrd0EvJWw1QiqVWrB/6m97jly/eZvoWAAQddPmLHEcM3IR5Q1kQdCbCVsipFAoshJkRlOdvJws0bEAIOoiX79hDndvyYomOhAAOiJsTaMSEhL5SXHDfcu1bQcRHQsAIq2KiigL/l7Ojv/zt1+JjgWAjghbjRAhJCsrO9vK4F4hjB0FgDD/3L478qel7kO1Lh7YqqKiQnQ4AHRECBMhQsjdgHyvSKAz8AAAOmPj9r0lA8a3vjhPdCAAfJtwJsIhKiSxpqqA6HSiAwFAFLExJG07R+vNpbWLYXJt0AcIWx8hrqmpqfzA5EUKWuJ71vzkOoPocAAQIaWlpfsfJAxyW//q5g54gAn0CcJZIySTyVIUxGAwKWTh3EEAei3bH90uBkVrhO6FLAj6CuGsEcrKyn54EzXMt0LPzozoWAAQLTSZflKlKSOd3YgOBAB+CWciRAipqKjMsVUMKGAPVYGJnQDoIW/K6eR191Mc6ox1NYiOBQB+CXPL4UwDcmAhDB0FoIfMmPeznZ3d2NyrkAVB3yLMiXCEGqm1LM/3WSzRgQAgEl6nZDCHzmAWJhMdCACdI7RNowih6urqmlOzvZR0lBgrp0+bRnQ4AAizRgYi/Xx9Mxa1eclqomMBoHOEuUYoISEhK4YxGquVlZSIjgUAYXbzbqCpo+tomZpDv3nBPDKgzxHmRKigoFCYHDdg8z0FczuiYwFAmG3686/y0Surn5wmOhAAukKYEyFCSFZWdtYwLf8CmHcUgG6kOnW9bsyJvb9Boyjok4Q8ESKEZhqQ/fJh7CgA3aKsrGytty/b4sfC188njHMgOhwA0k0MQwAAIABJREFUukL4E+FwNRItJXjHqassFovoWAAQNhPd5p6LK5YLWC8m/NcSILSE/+QtKCiofXjkSGDU48ePiY4FAGGDaZpKZj1zGD6Y6EAA6DphfnwCp66uriGFfSrNNDXbSnQsAAiV9PKmao8zKRNbTTQUiY4FgK4T/hqhnJxcUXKsyaHYBhUTomMBQHgsXfvbSMfJxlF7IQuCvk74EyHO01D8LowdBUBwopPSGUOmo0+ZRAcCwPcSlUQ4exDpbgGbBaNHARCENjZiL7m6frTGA9+LRMcCwPcSlUQ4UIZJv/WrvduixsZGomMBoG979CTEwMalX0XSid+W9OvXj+hwAPheopIIMzMzqdWlia3KERERRMcCQN/2x4FjZY5bPj88SnQgAAiGqCRCc3PziWaapKYqm7GORMcCQN+mNuM37agjx3b/QXQgAAiG8D8+gRMXF3/0zzm7IGZiM8VFmehoAOibPn/+fPxWUKGiU270UxlRuXgA4ScqNULc7EHkO/kwdhSALvJcuuZwbJXY1cWQBYEwEa1E6GlAflJAq22mER0IAH0Sqb+5eOoj+xE/EB0IAIIkWomQVV9O3Wc7aNiY/Px8omMBoI/Jr2zMsNkWHxl648IpomMBQJBEKxHW1NSIySq2yGiWl5cTHQsAfcnWPX9ZOLoo+a0api1HdCwACJhoJUILCwv/U3slJq8zHwWv6gWgE2JSc+imE8lVuUQHAoDgiVYiRAi5THBwmTw5sBCGzADALzaGaDO9V9gYhD+8Q3QsAAieyCVChNBcQ9ItGDsKAH8io6INbJxrE4LOb1qora1NdDgACJ4oJsIfdUlv/c/vOn4Ow2DuUQC+Ye+pi8Wj1zWFnCIRHQkA3UQUE+GbuBhGdtRR/4jY2FiiYwGgt1NzWaORcOnIbnidJxBaovhYrJmZmRbzcwUVmZqaEh0LAL1XfX39iX8Cwll22eGP1KSIjgaAbiOKiVBNTa3wXdTAu8xSMkWN6GAA6LXWbN5557OayvuFar8lEh0LAN1IFJtGEUIkhOYNIt3MhSEzAHyV2IAhlPSnVkPMiA4EgO4lookQIbTQiHzrfUsLFaZbA+ALiitqwgcufhYWFezvS3QsAHQv0U2EpKq86j/H6g0dDbPMAMDj+N8XTMe70k+5TtCTJjoWALqd6CbCqqoqcWXNZopCbW0t0bEA0Lu8+/CRpjtcsqWK6EAA6AmimwhtbW19/9ok7r5Xz9ic6FgA6F1qJm5f4jI2LuQ+0YEA0BNEcdQoh/uPE3zFWPeL2AuNRPeGAABubxITPVZtphpPCLq5Qxz+LYBoEPUzfaEhyRfGjgLwP+d9/UuGLWmLuwNZEIgOUT/Z7RUbo/fNn75oFYPBIDoWAIgnO2GZVkHokd1/EB0IAD1HpJtGEUJxryLZilph74uysrKGDh1KdDgAEKahoeHA2as3acMzH/vqyMLEokCEiHqNcMKECaPl6sX66Q02hyEzQKQdOHbqWEId03e1ljRMRg9Ei6gnQnl5+ejAazpLjiXWUIiOBQAiKZhak9OeWA21IJNF/bIARI2oN43iFhqRr+eybTQgFwIRlZWT54uNeRz11mUAZEEgcuCkRwihuQasm5fP3bh7j+hAACDA1es3R85cWbrbwVmX6FAAIAIkQoQQSgx9TC9IXr3vdEFBAdGxANDTPtW30iQU5chMeFU1EE2QCBFCaOjQocqf0xGZoqmpSXQsAPS0JLOlK375NT02jEKB3gEgiqCPECGEDA0NSzMTdW8zKlhiA4kOBoAe8y4paeqiNY2KA2uibklBEgSiCmqE/5Igo9mDyDDLDBApIWFRlQMnUcozWbQWomMBgDCQCP/fYiPS+YDnKalpRAcCQA+hWS8aJM++dGinrKws0bEAQBhoGv1/ha8e1ry85RhY+CH2mbq6OtHhANCNamtrvf7YF8I0zrywV08O5pEBIg1qhP9PS1NTuqGYwcZkZGSIjgWA7nXj1p17n5WwmH/k6PA+TiDqIBH+Pxsbm4zoFxLbY2hickTHAkD3krd0Ev8QOWGUhYqKCtGxAEAwaBr9Dz0N5akDWXfy2evM4RYBCCc2m/0k5PnebK0Hz8Jd+kOjKABQI2xnqTHp1KOYwsJCogMBoFsEBATO3Het4u8l1jLQKAoAQpAI26uIDiwMPDZi4ozGxkaiYwFA8MjKWsyqQjVpkrS0NNGxANArQCLkpaykKEWtaWOTxMSg3RgIm8bGxn+YI3bdeF6YmgCDwgDAwbWel/PkybFPLMZFyrHE4TIBhEpiYuLEuV6tTMw3OUpcXJzocADoLaBG+AWWg3Qm6MvdyYdZZoBQKS6rbFIYICNBobc2Ex0LAL0IJMIvW2qIHbv1NCcnh+hAABCYJG3ncXNWRt69rK2tTXQsAPQi0DT6ZcWhvnlhYTZ30kvS4qErBfR1ZWVljm7zi5jyuaF39ZSliA4HgN6F30SIYZivr++LFy80NTU3bNjQv39/ngKJiYmBgYGfPn1SUlLy8PAYN24cQojNZp87dy49Pb2lpcXMzGzlypWqqqp4+fT09PPnz7e0tHh6ek6ZMkWAuyQQA3S0pWvyGWQJ6EoBQiAhIaFA1lSyrqilvAApDyY6HAB6F36bRk+cOPHXX3/NmDGDxWLZ29vTaDSeAllZWSoqKtOnTzcwMHB1dX348CFCiM1mp6enW1tbT5kyJTEx0d7enslkIoRKSkrGjh3bv39/JyenxYsXBwUFCXavvp/Lj85vw5+KbYumIkiEoM9rNHVWVlHdPv/HwYMhCwLQDsYHJpOpq6v74sUL/OPw4cN9fX07KL9hw4alS5fyLKTRaGQyOTs7G8Ow7du3z5o1C19+7ty5sWPHfm1VPj4+y5Yt4ydIDMPYbHZTUxOfhfnh/pJ5PoslwBUKKxaL1dzcTHQUoojJZLa0tHRQgEqlLlq3Sd59Z0w5s8eiEgUMBqO1tZXoKEQRg8GgUqmCXSdfNcKSkpLS0lIHBwf8o6OjY3x8fAeFX716ZWtry72QyWQGBATo6OgMGDAAIRQfH89Zm4ODQ0JCAoZhXczk3WmhbsveIyfDwyOIDgSALgoPD7+d2dhWkiX5KZnoWADopfjqI6yoqJCXl5eQkMA/qqmpJSYmti8WGBg4c+ZMhNCsWbMWL17MWW5tbf3mzRsFBYUHDx7gk1lUVFRwOgvV1NTodHptbS1nCbeioqKXL1+6ublxlmzbts3U1PSLcWIY1trays8e8eltwNnyqoaffv6l4G2kpKSkANcsZNhsNo1G6513M8KNxWLR6XQ2+6uP+mC6Q1DtuRFaUnp6es3N8NSEwDCZTAaDwWKxiA5E5DD/h8/yMjIyZPI3qnx8JUJpaWk6nc75SKPRvjiQ0sPDA8OwgoKCRYsWbdu27dChQ/jyhIQEOp3+4MEDNze31NRUPT09aWnptrY2/Lf4mr82MlNDQ2PQoEGzZ8/GP1IoFBMTk68Vxiu5AhzkOd7e5tQ/G6W0BikrKwtqnUKJzWaTSCQYXtvzWCwWhUL52pEPuPdgVyrZ5+6jBYYwubaA4YkQpqnreXgWlJLid/DzN7Mg4jMR6ujotLW1VVZWamhoIISKi4t1dHS+VnjgwIGrV68+ePAgJxEihCQlJWfPnn348OHo6Gg9PT1dXd3i4mL8V8XFxaqqql87n6SlpQ0MDDw9PfmJE8MwMpnMz27zydHBoSDjnVEAs7aNrAZjzjsk2CMP+NTBOR8bG7vk0DUGgznYtT+ZPLLnYxNu5P8hOhCR0x1Hnq91qampOTg43LhxAyFUU1MTHBzs7u6OEKqurvbz88PLVFZW4j9gGPby5UsTExOEUENDA6cqWVBQkJeXZ2xsjBByd3e/e/cuXim8ceMGvrbeSVWKNEOPfO0DzDID+hgZDT16bXk/VnV/XV2iYwGgV+P3OcLDhw9PmzYtNDQ0JyfH1dXVysoKIfThw4c5c+bg7ZYzZsyg0Wjq6uq5ubny8vKPHz9GCCUlJXl6euIjtlNSUjZs2IB/0dPT09fXd9iwYf369SsuLo6MjOym3ROI6ZK5c5ZuaXSz27dtE9GxAMCXgoKCXa/bttx9vW8EhUKhEB0OAL0aif8BDo2NjUlJSRoaGmZmZvgSOp1eWVmJDwRlsViZmZmfP3/W1tY2NTUlkf7tk6iurs7OziaRSKamptzDYTAMS05Obm1tHTVqVAfjUC5fvpyQkODj48NPhBiGtbS0yMkJ+P3yv2//07tURy7iZH1BGlxTvoj9f+3dd0ATd/8H8EvCDJFhIDIUZE8BGW4RBAQHoqg4oYJFWqs+2opVa6s+1r1aWwSV1m0diCgUKiggKggoKsoSRUSQPRNC5t3vj3t+PD64QIEvST6vv+6Ob5J3Lsd98v3e5Q7H29vbVVRUUAeROeTJMp2OET59+nSEzwKOWC7j4pERw21RZZNucIwQle4eI+yKblxiTVVV1dXV9c0lioqKZBXEMIxGo9navuNfTlNTc9y4cW8vp1AoDg4O3UiKTujihVeXrsW8gqEKAokgwmgcAa4iL1JRgC0WgI+Da41+nKmp6ePrlw3OCZ+2EGZqcPYd6NcEAsHp+kHeu65EjqEMhqODAHQBnPLUJYo0LMiMevBeM+ogAHxIbW2t/rARuxaMX2vWBlUQgC6CQthVmnnnI0InO7pNRh0EgPdqbGpqxBiKTD28uRp1FgAkBgyNdlVdeSlmPOpZyQ3UQQB4r0yKqUnQtj32vE6H8wEAHwCFsKs2r18jf+zyaWEogWFwnBD0N21tbaMm+RbWtsVeODt1uCHqOABIEhga7SolJaWtX89XHWx6vRKuqAn6nYqKipJWjKpnLSx7iDoLABIGCmH3LNau/2bj9ozMTNRBAPgf1apmDJeADZNtfHx8UGcBQMLA0Gj3PDyzp4Sr4xsYWleShzoLABiGYSKRaNr8oOtF1VFHoxY5vfciwACA94EeYfdMn+TKuH9W3XI06iAA/EdFRUVKQYVQzx4rghtnAvApoBB2zyxfn+K8+w1zfm/kf7wxAH2gTd2QYuHmO4gzffp01FkAkEgwNNptunTK1CHUw4/b1jvBpTUBSgRBfPvDlpN3X/y8fc/qMdqo4wAgqaBH+Cl0co/9ONdt0qyFqIMAmVZVVRUVl9Yy0Iz+JA51FgAkGBTCT1FT8gS38X6QX4Q6CJBpIlVtXN/RriVn2tQpqLMAIMFgaPRTHNy5hRF5Md37DOogQHb9cujI/qT81d//sHEEo9NtmAAA3QI9wk+hpqb229ovhZom6dXw43qAQFNT04+/RFVQNdVz/0KdBQCJB4XwE1EpmB/t0ex5C/+Kvow6C5A5YmV1ka6NwauUad6eqLMAIPGgEH66+38drLOdu2rDZtRBgGw5ez7aJeTHoB92v8hJNTc3Rx0HAIkHhfDThX0dpJV+wNAnFHUQIEO4XO5XG7YWceRZ90+gzgKAlIBC+Ok8JroWZaU9sw+p4qKOAmQGoUAX69oOevaPzyQ31FkAkBJQCD/LQEXMj9UcsvvP8vJy1FmA9ItLSHT5+t9+6w9UPc50dHBAHQcAKQGF8HOV/hGWkFfp4Qc/rge9SyQSLfpmTW4Dxso6gjoLAFIFCuHnGudgq/AkQVEPzlkAvUtEkcOHDGcWxs3wnIA6CwBSBX5Q/7k2r/t28qKlM9MVeWJMiYY6DZBSyddv/Hj+tue3+2NmsFBnAUDaQI+wB4wczLDXwLfG5ggEAtRZgBQiCGLOlyuzmhS07vyOOgsAUggKYc/gn/rXrgO/TZzujzoIkEJCgoIbjdB4fGnOZDhTFICeB0OjPUNTiUKhUOracdRBgLRJTEr+6WzK2NCfE+fC3ecB6BXQI+wZZ4/+vn/dNwrL/oJrj4KetfCr1feEOsz0X1AHAUBqQSHsGTQabcXUEcpKilF3SlFnAdJDiGMUS1dm7ulFM7xQZwFAasHQaE9SvrBqWUHZ31as2DN/os4CJN7luL9/OpHoOG9VUoAJ6iwASDPoEfYk+bZ6XMeyuKIWdRAgDb5cve6JujMjFQZFAehdUAh7UuzpqLBZEwZ8fRJ1ECDx+GKM6uindTfyq4V+qLMAIOWgEPYkBoOxPXhau9yAE5nPUWcBEuzYmXOmHv42rtNrn2ROcp+IOg4AUg4KYQ+jUjClM6Ehy1fPDf4adRYgqcI2bX81bIFCWjjqIADIBCiEPU+Z3yzWNC553YA6CJBIXBFGnbjM4O7vG1csQZ0FAJkAhbDnJVw4uXaBt0LwUdRBgOTZ+Uu4nv04GxODsuzr48eNRR0HAJkAhbDnMRiMbYFerSJs98VUPp+POg6QJPsioprd1+IZp1AHAUCGQCHsFVQKJo5cuCHyoj8cKQRd1izAKDO3Or+89Ou/16HOAoAMgULYWwZrqFA4dS2EMuogQDIsX/ujjqWDOV6ZHXvCztYWdRwAZAhcWaa3JMWcjUwtPNJijhMYlYI6Dej3LiVe57ksoxYmYVgo6iwAyBboEfYWGo32jYcNpaZ4/ro9lZWVqOOAfu3e03LhosjFrJpT4XtRZwFA5kCPsHfVRYY8sgtoWB52/fJZ1FlAP+Uzf3FSYa21lcWxs/tRZwFAFkGPsHdNmeiicPvowGHjUQcB/VdeWY1Ib5gWvwZ1EABkFPQIe9fRX3aFbNjud4PgijA6rGzwvxobG0/G3eD4HwzXLVw0FS6lBgAa0CPsdSNYNGZWlKGTa3JKKuosoH+ZOm/xdzGPmVfCvpo7ncFgoI4DgIyCQtgXqhMja13X7Pj9D9RBQP/SMkBf/vmd0ZZDUQcBQKbBaF1f+GXbprDw00OC1qAOAvqL3AcPv9sX1Wi/4NGvJuaDWajjACDToEfYF+bPnpnz97kEwjo97ynqLKBfCF69Po05Sf5SGFRBAJCDQthHdOiY3G/TJwV9u37LdtRZAGIEholHzFe99vPqkC9QZwEAQCHsQ6oYT8jQflbTgjoIQOnI8dMa5k7NddUNBdnffrMUdRwAABTCPnQ3+eq/ghfUe21GHQSgFHEmpsX9e+r9S3LwzwdA/wD/i31HQ0NjT+DEZzdjHKcHFhUVoY4D+hqO4xn3HjZ7/+TKvXv+6G+o4wAA/gPOGu1TNArGid1RMemHTXt/Px/1O+o4oE8tXbX29L0KNVFLanYi6iwAgP+CHmFfW/7lF4w7h1mu81EHAX2tsFEkoinpKolRBwEA/A8ohH1t64awB3dS/+KaJ97KJggCdRzQF0pLS2eFrikwmHbyx5CMa1dQxwEA/A8ohAiYqFKEe71nrN37826424BMWLlxW4x4GOVi2IJJY5WV4V7NAPQvUAjR0KTLifi8OgEco5V+BEGwrXxU0sMDZ/uizgIAeAfYEaPxOCNlf+qzRNyKwDC4fb0UO3n2/IpNuyhDneoKspXhvw2Afgl6hGjQ6fQNU2wLDyxhWjinpN5EHQf0lgtJt1sd5iuWQxUEoP+CQogMBSPw0ntNNn6Xk+D2TFKIzWb/evjPh8b+cw2xf84dQx0HAPBe8DUVGQqFcv7or+v+ut3usgx1FtDz1m3ZEVGM0QsOnnv+EHUWAMCHQI8QJS9Pj5SIzRfOnhrps6CsrAx1HNCTqgfayBUkO1gYow4CAPgIKISIMahC/NaxbPWxh4+fQZ0F9IwbKWksqxGJKelF2TfT/76EOg4A4COgECImLy8fsnD2gIK/ac5+qLOAnnHqalKd7VyF57eMmHTUWQAAHweFEL0D2zbfT/n74LW8kLCfmpubUccBn47L5R46eixjoKsLo/H0oX2o4wAAugROlukXdOW4ovgdf9j6GkcdW7dmNeo44BPt/uX3rTerFIr2FZY+otFoqOMAALoEeoT9Ap1OHz/cRr4giWrhgjoL+ERisfiVug21JN3aaAhUQQAkCPQI+wUKhXLt0umsWmLy7tjGR6k/rV5Gp8PhJUmSkpo2O/RbjjKr6N4d44FKqOMAALoBeoT9iKlcU3vstj23q4+dOIU6C+iexMyHTYaudHalHnyBAUDSQCHsR9TU1OyN9KhPbwqGOqPOArqqoaFh0TffnS5T9B9lGn8qUkkJuoMASBgYGu1HaDRaZtKV9GpixpZjz29e3fVjmIqKCupQ4CP+PHnmbJW6Yt6fh+8lq6uro44DAOg26BH2OyZ4VXvanxGP2s6fv4A6C/iInJycTHlr+tPrXqPs1NTUUMcBAHwK6BH2OywWy8lIO+t5JtdwCeos4ENiYi4H/BwlqH/1IvfWYE2oggBIKugR9jtycnK3/o6+mXbz+/1/esz5orGxEXUi8A4EQRRwlXht7IFKmLoS/FgCAAkGhbCf0uK8IF49SuUNTkhMRJ0FdJafn69t5bRl++6/jh8puXebwWCgTgQA+HRQCPspQ0NDv5EWjPri1wYeqLOAztJyHtWxHJQI/pghDFVVVdRxAACfBY4R9lM0Gu304V8fV7EdPHwP0rh34s4bGBigDgWwxsbGOUtWPBQyfcbaB4/2HTx4MOpEAIDPBT3Cfk2+pVJRjlqlbpWdCzd37RcSk5LTBEO4r5/vDJjo6zMNdRwAQA+AQtivWVhY/PZdkJ2dXSzfrKCgAHUcmSYWi3ceOPhrVh2T83LOKDMzMzPUiQAAPQOGRvu7oMCFDiOfO3r7x9GVbxw/4OwMF51B49q1a5su38ebazJP7nGyH4Y6DgCgx0CPUAJoKMuryonb2vk1QnnUWWRUUnLyueI2rKbETEVgbjwUdRwAQE+CHqEE0NfXL85IjnjI/mLHwcl62B+/7lJUVEQdSoZkZGTMDNvL4/GST0dOdLZFHQcA0MOgRygZtLS03JQqOU0NF4pa79y5gzqODCksLPynpJXXXKuFtVoO1kIdBwDQ86BHKDEcHR3Hs/CsGm6GeMiwujotLdgp97r8/Pwxs5dwRNSTUYdnj7eDjjgAUgkKocSg0+nXY84kZuVPmzt3vzx+P+myoaEh6lDSrLKy8npxHaedryqPOQ1RgyoIgLSCQihhdBSFDGWlFoHoVQsPymDvefHihcMkv1ZcYc/evYtchrFYLNSJAAC9BQqhhLG3t8+8FBX+gD3tqx/G6DOunDoKPZUeV1VVdS3vZauQwlCmuBkPhCoIgHSDQih5rKysvEvjjwzQTS2revzkiZOjI+pEUqWysnLYhCktFJUNP24KchtmZGSEOhEAoHd146zRgwcP2tra2tvbHz58+O2/PnnyxNfX19ra2tbWdvny5eTNgwiC2LFjh5ubm4WFxdSpUzMyMsjGtbW1nm84d+5cj7wZ2eE+ceKsoRSNISabY+/HXo1DHUd6lJaWXrqd1yKkqMhhU60HQRUEQBZ0tUd4+fLlPXv2xMbGikQiX1/foUOHenl5vdlAXl4+KCjI0tKSx+OtW7du6dKl0dHROI7n5+dv2LBh6NChV69e9fLyKiwsHDx4MI/Hu3PnztWrV8nHmpiY9PDbknZ0Ov3c0d9ORl9dvP/izcvluZYWpqamqENJvOrqaicvvxZFre/Wfr/Uyxk2SwBkRFcLYWRk5OrVqx0dHTEMW7FiRWRkZKdCaG5ubm5uTk4vW7Zs1apVGHkLhdOnyYXfffddREREdnY2ecF+Go3m4QE3GPoso+0sWc2F9SLFpKeNKiqvdXV1USeSYDk5ObF3C1v4OEOO7WuvD1UQANnR1UL45MmT9evXk9OOjo4nTpx4u41IJHr06FFLS8u+ffuCg4M7/bWqqqqiosLKyoqc5fP57u7uCgoKkyZNWr58ubw8XDys20xNTSsfZ/0afWPlmnU/iRqLM29oamqiDiWRamtrPRZ+xR5otmb16pVzPOHmSgDIlK4Wwvr6ejU1NXJaXV29trb27TZsNjs0NLSlpUVRUXHGjBlv/kkgECxatGjp0qUWFhYYhjEYjIMHD9rZ2VVXV2/cuLGgoODo0aPvfN2ioqJTp05FR0f/J66cXHR0tIODwzsbEwTB5XIJgujim5IO5gwhnSJq5uE5pdUjqVQkXylwHOfxeDiO9/1Lf74Lly7HZBW3CTG1xuIpw0PV1NTYbDbqUF0lFosFAoFYLEYdROaIRCKhUCgSiVAHkTkikYhc+V1sT6fTaTTaRxoRXaOlpXXr1i1yOikpaejQoR9oHBERoa+vLxaLyVmhUOjn5+fr6ysUCt9unJ2draCgIBAI3vlUR44cCQgIaHwDjuPve10cx9lsdhffkTQpKSlZdiiWZuikZ+lAfhXoY2KxmMPh9P3rfr76+noVUydsctj63Yck8S2IRKK2tjbUKWSRUChE8r8GhEJhe3t7zz5nV88aNTIyKi4uJqeLi4s/fDbd9OnTy8vLW1paMAwTi8WBgYE8Hu/8+fNycu/ogGpqagqFQj6f/86nolAoioqKGm+gUChdzCw7TExMXAbyaSzDmjZxSVUjIWN94k+2bM0Ge59AkYqWXkX6gsnjVVRUUCcCACDQ1UIYGBh46NChtrY2Npt9+PDhgIAAcvnatWtLS0sxDMvJyamvr8cwjMfj7d+/38rKSkNDA8fxoKCgsrKyqKgoLpfb1NREFrzHjx9XVVVhGMbhcDZu3DhmzBgGg9Er709m+PvPubI52HvxSkeP6aaOY2HE5qNq6+qPxfxTYeD+7RyPiry7NjY2qBMBANDoaiEMCQmxsbHR09MbMmTI6NGjOwrh2bNna2pqMAzLzs42MzNTU1NjMpkPHjwgj+q1trbGx8cXFRVZW1sbGxsbGxuTJ5E+ePDA0tJSVVV10KBBzc3NHWeWgk9GoVC8vb19hmDYELuXDW1FZVVw3Oh9CIKwG+ehP8pL3dzZXfBgWaA/6kQAAJQo3RpG4/F45Fjl+xqw2WwVFRUqtUv1tbW1lcFgfLhxVFRUVlbW+06l6YQgiLa2NlnuXIpEouhLMYfvvr4d95c+g1KSm9HFz+Iz4Tje3t4uEUOLJSUlqdkPv1m3hbCYuHmS4caw1agTfRaxWMzn8+l0OuogMoc8X0NZWRl1EJlDniworCWdAAASIklEQVSjpKTUg8/ZvUusffS1BwwY0PVnU1VV7darg4+Sk5ObN9e/rvHwHV3L8rL7yRn3Rw0z6zjdF+A4Ptp7ZqP+WOdx7ku9HBbOm4s6EQAAPbjWqBT6OmTJ0MG6+xJNpq7YwuK+qih80Df9wn5u8879p+KuN1NUBlTnrQxdvXAejIgCADAM7lAvleTk5Hx8fPyH6yjSKLVc/If9kffv30cdCiU+n19QWLT7UFSp1bzp7uPLMhKgCgIAOkCPUGotC1kybqTz2kPndqa9igwPbyx9IrO/PLEd6/ECVxsw1M6uMm7bvn9raGigTgQA6EegRyjNbG1tl04erf40iUNnGY+dsm3fQdSJ+tofp846TvZ/Xs+hqgxc4madmXDR0tISdSgAQP8ChVDK+fn61ObfnTZ2+Ithi7b/GvmqokJCL4TWXRwO5/adO6t+2pFrPG+4ve31HUt3btqAOhQAoD+CoVHpJy8vv23Vl1Xf/VRi5GAyebGVBu1B+jXUoXrdCE+fZ/IGikx9i9xDv+zdOnbMaNSJAAD9FPQIZYKVldXdxOh5zgbEINP8l1VTApfFJ0ptLfzh511als5Pa9vkBJwv3IYXZlyHKggA+ADoEcqQA9s2Tb1+fV04O5E+PnPF2qwEQ2Nj449fl11yZGXn/H7iQkxiMtd7w/DnF07s3TRs2DDUoQAA/R30CGWIgoLClClTti8P0Mv4lc/Qtl0QNtZ7JupQPaO1tTUhIcE3eMVpbKQyQzWAknXut+1QBQEAXQGFUOZMm+xVkXfXd4Q5oa57/1mFnu2YfeFHUIf6LJWVlRP9FvqGpzUJKFopO39e9eXJ8H1mZmaocwEAJAMMjcqo47/t+fL27W+2VhYP/2rTnu+LnpctC/Qfbm+POlf3iMXigK9XX8152t5Up8Qkpk9w+uvob6hDAQAkDPQIZZSioqK7u/uZ/Vu8qy/JMdSjePZTFn114cKFhoYG1NG6atPOvUzLEReTb/OGONoZ6RTEhEMVBAB8AiiEMs3RYXjiX39sXv4F88bOhnbR/L8KJs4OvH37tkAgQB3tQ8KjjjPNHHYePdcyMlhbVTk62Dkl5rSBgQHqXAAAiQSFEGCrvvqyvjh3vpcL/Vnak2evPDadmjw/ODs7u7/d0VAoFEadPKNp7rBm5++N474ZICfeYS+8efnUjBkz1NXVUacDAEgqOEYI/uPEof37GxpmBX+TXd+aevfRhLBD7nq0fwXNc3Fx+cAdKPtGeXl5/LUbP+wN5wrEApdQ9duHvtap+OLCcTs7O7TBAABSAAoh+C8mk3kj5kx+fv6Sbzc+4YkSUjOT65XHHjrm6z5u0fy5TCazj/Pw+fyEhIRb9/OOxt9ur6vARy5Ue3RhpX7d/IvHHRwc+jgMAEBaQSEE/4NGo9na2t5JiH7y5Mmy9f9+XF+WllecJjfs0LkAA23mL/9eb2lp2dt3sSAIorKy8rtN21/VNuW0DxTnJdKs3Zm8hg0TB07adh6umg0A6FlQCME7KCgoODg4ZCTGlJWVrd60IzXr0vPWlqe2i0dMD1CkEltWh+oO0po6dWrPDpmKRKLLl2Nb2K3rd/3OF+HsMSHUrH1KTJ2hJoYH1s+xt9/GYrF68OUAAIAEhRC8F5VKNTIyunLqqFAo3HEg/GR0xGs+t2nkopU/7aDZTtHbdUiBEO1YtxLDsGlTp37CmTV8Pp9KpZ6/GE0h8E0HDvOF4mo1c7wwFR8xT/FhtH7+ea8Fc3b+9L2GhobM3kkRANAHKARBoM7wIVFRUVlZWUePHu1KY4Ig2traGAxGb6eSWTk5OddS009GX33NFvCaG8SuX1Fv/Ea1cFN5nUujUlzHjuS38/ynejazOeNGObe2turr63O5XAUFBaFQ2NjYiGHY/bx8KobF30gX4eKHBSVCbmv70DHE09v4yAVyDy8rKiprqtBCA+aNcnJwc52A+u1KALFYzOfz6XQ66iAyRyQSCYVCZWVl1EFkjkgkEolESkpKPficUAhBtxEE0dzcfPFK/PX0zMyce03tYgGfJxo6Eqt8TJi7UgtvYDoWtMrHmCqLyq6lyCsRYhFG4ISyGs5tIXQsiZoS3Gw85VkmRdNArqpQQZ6qrqzg5e7qaGsdMHc2nU6nUuFXPV0FhRAVKISo9EYhhKFR0G0UCkVDQ2Pp4oCliwNwHOdwODn37uU+evykRDX/aVGDhjK7vkBIxXFuPY7hmFhA4GIKhUJtb6RSCOXW58rqSsYqdYO9xniOG6Gn983oUaMUFBRg8BMAgAoUQvBZqFSqqqqqm6vrqJEjVVRUyIU4jvP5/Pb2dhzHeTyeoqKiWCweMGAAhUKBvgsAoL+BQgh6HpVKVVZWhlEjAIBEgIMxAAAAZJpUFcKKioozZ86gTiGLXrx4cf78edQpZFFxcXFMTAzqFLLoyZMn8fHxqFPIotzc3GvXrvXsc0pVIczPz4+OjkadQhbl5eXFxsaiTiGLcnNz4+LiUKeQRdnZ2QkJCahTyKLMzMykpKSefU6pKoQAAABAd0EhBAAAINOgEAIAAJBp/f3KMhs2bAgPD9fU1OxK4/b29paWFm1t7d5OBTrhcrlsNnvQoEGog8ictrY2LperpaWFOojM4XA4PB6vi7sm0IPYbLZAIOj6XeEWLFiwdevWD7fp74VQJBI9f/5cXl6+i+35fD7yu8jKIIIghEKhgoIC6iAyB9Y8KjiOi8Xiru+aQE8Ri8UEQcjJdfVH8Do6Oh/9TXN/L4QAAABAr4JjhAAAAGQaFEIAAAAyDQohAAAAmQaFEAAAgEyTsLtPpKWlxcXFaWpqBgcHv/Nk/ZcvXx47dozL5c6ZM8fZ2bljeXR0dEZGhr6+fkhISMfdgkDXZWdnR0dHq6ioBAUF6evrd/prQ0NDXFxcfn6+mprarFmzLC0tyeVxcXFVVVXktIaGxpw5c/o0tFRITExMTk7W0dEJCQlRV1fv9NdHjx5lZWV1zM6bN09VVRXDMBzHz5w58+DBAxMTkyVLlsCp1J/g5s2bV69eZTKZwcHBb/8oKzk5+cWLFx2zDAZjwYIFGIZFR0c3NjaSC7W1tadPn95ngaXDixcv7t+/39jY6O/v//YGT4qLi0tNTdXT0wsJCSE3eAzD6uvro6Kiamtrp0yZ4uHh0a0XlaQeYXR09Jw5c/T19UtLS0ePHs3hcDo1qK6udnJyam1tZbFYnp6e6enp5PIdO3b88MMPJiYmKSkp3t7efR5c4qWlpU2aNInFYjU3Nzs7O9fU1HRq8O2338bHxw8aNKipqcnJySklJYVcvm/fvn/++ae0tLS0tPTVq1d9HlziHTlyJDQ01MjIKDc3d/z48SKRqFODpKSk8PDw0v/X0SAsLGzv3r2mpqYxMTFz587t8+AS7/Lly7NmzdLX1y8rKxs1ahSbze7UoKampmO1Hzp06PTp0+TyLVu23Lx5k1z++vXrPg8u2err6x0cHCIjI0NDQzu+Q3dy8ODBlStXGhsb3717183NTSwWYxjG4/HGjBlTWFhoaGgYEBDQ7bsvEJLDwcHh+PHj5PSYMWOOHDnSqcHmzZtnzpxJTu/atWvq1KkEQbS3tzOZzMzMTIIgBAKBjo7OzZs3+zC1NPD29t69ezc5PWPGjK1bt3Zq0N7e3jG9atWqefPmkdMTJkyIiYnpm5DSRywWGxoaxsXFkdOWlpaXLl3q1Gb37t1LlizptLCxsZFOpxcXFxMEwWazGQzG48eP+yaz1HB2dv7jjz/I6fHjx0dERLyvpVgsNjAwiI6OJmdtbGxgD/PJcBwnJzAMKygoeLuBUCjU09NLSkoiCEIkEhkbG8fHxxMEceLEieHDh5MPv3DhgrW1dbdeV2J6hBwOJzc319PTk5z19PS8efNmpzbp6emTJk3q1CA/P18oFI4cORLDMHl5+QkTJrz9QPBh71yxb1JSUuqY5vP5DAajY/batWv79u37+++/CfjFajeVl5eXlZWRgzxUKtXd3f2dm25JScnu3btPnjzZ0tJCLsnOzmaxWGZmZhiGMRiM0aNHd4yOgK7gcrk5OTkf3tt0SE5O5nA406ZN61gSGxu7f//+5OTkXg8qdSgUyocbPH/+vLa21s3NDcMwGo02ceJE8qNJT0/39PQkH+7p6Zmfn19XV9f115WYQlhVVUWhUDouJTVo0KC3hx2qqqo6GrBYLA6Hw2azq6urtbS0OtbvOx8IPqClpeXNi3ixWKz3DVlgGPbgwYMzZ86sXr2anLW0tFRSUqqtrV21atWUKVNwHO+LxNKiqqpKVVW140vGOzddJpNpYmLS2tp64sQJS0vLly9fYhhGbvMdbWCb7y5yC2exWOTsh1fgn3/+GRgY2HEU1t7enkqlvn79Ojg4mDxqCHpQdXW1hoZGx2VlOj6aN3f+6urqCgoKH9hNvU1iTpaRk5MjCKLjmkZCofDt4/9ycnIdx0jICTk5uTcXvu+B4APIFf7min3fBb2ePXvm4+Nz4MABKysrcklERAQ5sWHDBjMzs4SEhDe/OIMPk5eX/+imGxwcHBwcTE7Pnj17586dERERcnJy5IGTjgfCNdi6pWObJ1f4B3YaDQ0NV65cyc7O7lhy6tQpcmLNmjWmpqZ3794dNWpU70eWFe/bn7+5nBxi7dY2LzE9Qh0dHSqVWllZSc5WVlbq6Oh0aqOnp9fxxa2yspLJZCorK+vq6tbW1gqFwg88EHwAnU5XV1d/c83r6uq+3YwcxNu0aVPHfvlNampqNjY2b55lBz5KV1e3ra2tubmZnP3opjt27NjS0lLyga9fv+4Yi37fRwbeR1tb+6N7G9Lp06ft7e1tbW3f+SRGRkbkJwJ6iq6ubnNzc1tbGznb8dG8ufOvqqoSi8Xd2s9LTCFUUlLy9PS8ePEihmECgeDq1avkeclcLjclJYX8LuDj43Pp0iVy/O3ixYs+Pj4YhllbW+vo6MTHx2MY1tDQkJqaSi4HXefj40OueRzHL126RK5AkUiUkpLS3t6OYdirV688PDzCwsJCQkI6HiUSiTq+o1VWVubm5lpbW6OIL6l0dXWdnJyio6MxDONwOImJieQ239LSkpqaSrYh1z+GYWKxOCEhwcbGBsOwMWPGiESiW7duYRj26tWr+/fvT548Gc17kEwKCgre3t4de5srV66Qa769vb1jb0M6fvz4m9/8BAJBx/h/cXHx06dPYZvvESUlJfn5+RiGGRoaWllZxcTEYBjW2tqalJREfjQ+Pj7x8fHkv0N0dLSLi4uamlo3XuBzzvDpY3fv3mUymQsXLhwxYoS7u7tQKCQIoqCgAMOwpqYmgiA4HM7w4cMnTJjg7++vra1NnjVHEMTFixc1NTUXL15sbm6+dOlSlO9BMhUWFrJYLH9/fxcXFycnp7a2NoIg6uvrMQwjV7Kfn5+ysrLj/1u8eDFBEC9fvtTW1p45c6a/v//AgQNDQ0MRvw0JlJSUxGQyAwMDbW1t/fz8yIW3bt2iUqnk9IgRI7y8vBYtWmRmZmZnZ1dfX08uP3LkyKBBg4KCggwNDdevX48mvSTLyckh9zYjR450c3MTCAQEQRQXF2MY1rGSc3JylJWVyZ0PKTc3d8iQIbNnz549e7aqqiqs+U/g6enp6OiIYZi1tbWjo2NLSwtBECtWrOg4Fz0+Pl5TU/OLL76wtraeP38+uVAsFk+ePNnBwSEgIIDJZN66datbLyphd5+oqalJS0tjMpmurq7k8VIej5eXl+fo6Eij0TAM4/P5KSkpHA7Hw8NDQ0Oj44GlpaVZWVn6+vpjx45Fll6SNTU1Xb9+ncFguLu7k4PvIpEoNzfX1tZWSUnp2bNnHacsYhimoqJiYWGBYVhRUVFBQQGO48OGDTM3N0eWXpJVVFTcvn1bW1t7woQJ5DlfHA6nsLCQvF5EbW1tTk5Oa2urgYHBqFGjqNT/jvEUFhY+fPjQxMTkzStLgK6rra1NTU0dOHCgm5vbm3sbBwcHcra2traxsZHc1EkEQeTn5xcVFVGpVHt7eyMjI2TpJVZeXl7HkSwMw+zt7Wk0Wnl5uVAoNDY2JheWl5ffuXNHT09v/PjxHSdCisXitLS0+vp6FxeX7h7/krBCCAAAAPQsiTlGCAAAAPQGKIQAAABkGhRCAAAAMg0KIQAAAJkGhRAAAIBMg0IIAABApkEhBAAAINOgEAIAAJBpUAgBAADINCiEAAAAZBoUQgAAADLt/wCDdu9F691MIAAAAABJRU5ErkJggg==", "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/html": [ "" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "q = 1\n", "f(x) = x^2 \n", "n = 2^8\n", "h = 1/n\n", "X = 0:h:1-h\n", "\n", "# write down the matrix ℱ (exercise 1)\n", "ℱ = # your code here\n", "Fhat = ℱ*f.(X)\n", "\n", "# write down the Fourier coeficients \n", "# of the numerical solution (exercise 2)\n", "Uhat = Fhat ./ # your code here\n", "\n", "# leave the following as it is \n", "# (in Julia A' is the conjugate transpose of A)\n", "U = ℱ' * Uhat \n", "U = [U..., U[1]]\n", "plot( 0:h:1, real.(U); \n", " m=1, legend=false,\n", " title=L\"Numerical solution to $-u''(x) + u(x) = x^2$ with PBCs\" )" ] }, { "cell_type": "markdown", "id": "a32395c5", "metadata": {}, "source": [ "**Exercise 4.** What is the computational cost involved in solving $\\mathcal L_h U = F$ by using the code above? " ] }, { "cell_type": "markdown", "id": "cb37bc49", "metadata": {}, "source": [ "::: {.box} \n", "\n", "***Your answer here***\n", "\n", "\n", "\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "1a124775", "metadata": {}, "source": [ "A Fast Fourier Transform (FFT) refers to any algorithm that improves on this computational complexity. In order to get the main idea, it is useful to first consider the case where $n = 2^m$ is a power of two. Let \n", "\n", "\\begin{align}\n", " E_k &:= \\frac{1}{\\sqrt{n}} \\sum_{j=0}^{\\frac{n}{2}-1} U_{2j} e^{-\\frac{2\\pi ijk}{n/2}}, \\nonumber\\\\\n", " %\n", " O_k &:= \\frac{1}{\\sqrt{n}} \\sum_{j=0}^{\\frac{n}{2}-1} U_{2j+1} e^{-\\frac{2\\pi ijk}{n/2}} \\nonumber\n", "\\end{align}\n", "\n", "be (up to a constant) the Discrete Fourier Transform of the even and odd terms, respectively. \n", "\n", "**Exercise 5.** Show that \n", "\n", "\\begin{align}\n", " \\hat{U}_k &= E_k + e^{-\\frac{2\\pi i k}{n}} O_k \\nonumber\\\\\n", " \\hat{U}_{k+\\frac{n}2} &= E_k - e^{-\\frac{2\\pi i k}{n}} O_k\n", "\\end{align}\n", "\n", "for all $k =0,\\dots,\\frac{n}{2}-1$." ] }, { "cell_type": "markdown", "id": "f03acf6f", "metadata": {}, "source": [ "::: {.box} \n", "\n", "***Your answer here***\n", "\n", "\n", "\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "f4ad98b5", "metadata": {}, "source": [ "This is the basis for the (radix-2) Cooley-Tukey algorithm where one computes the Fourier transform on $\\{0,1,\\dots,n-1\\}$ by iteratively subdividing the domain into two parts and computing (rescaled) Fourier transforms on each smaller piece (this is an example of a *divide-and-conquer* algorithm).\n", "\n", "**Exercise 6.** When $n=1$, what is the discrete Fourier transform of $\\{U_0\\}$? Hence, complete the following code sample:" ] }, { "cell_type": "code", "execution_count": null, "id": "1584ad4f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "FFT (generic function with 1 method)" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function FFT(u; normalise=true)\n", " n = length(u)\n", " if !ispow2(n)\n", " @warn \"n = $n is not a power of two\"\n", " else\n", " U = zeros(ComplexF64,n)\n", " if n == 1 \n", " # your answer here:\n", "\n", " \n", "\n", " else \n", " E = FFT(u[1:2:end]; normalise=false)\n", " O = FFT(u[2:2:end]; normalise=false)\n", " for k ∈ 1:Int(n/2)\n", " t = exp(-2π*1im*(k-1)/n) # Twiddle factor\n", " U[k] = E[k] + t * O[k]\n", " U[Int(k+n/2)] = E[k] - t * O[k]\n", " end\n", " end\n", " return normalise ? U / √n : U\n", " end \n", "end" ] }, { "cell_type": "markdown", "id": "930c035e", "metadata": {}, "source": [ "**Exercise 7.** By using the function ```FFT()``` and the complex conjugate ```conj()```, implement a one-line function that computes the inverse discrete Fourier transform:\n", "\n", "\\begin{align}\n", " U_j = \\frac{1}{\\sqrt{n}} \\sum_{k=0}^{n-1} \\hat{U}_k e^{\\frac{2\\pi i jk}{n}}.\n", "\\end{align}" ] }, { "cell_type": "code", "execution_count": null, "id": "95535500", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.466075386898995e-15" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "iFFT(u) = # your code here\n", "\n", "# this should be (close to) zero:\n", "norm( iFFT( FFT( f.(X) ) ) - f.(X) )" ] }, { "cell_type": "markdown", "id": "4095de7c", "metadata": {}, "source": [ "**Exercise 8.** Look up another use for discrete Fourier transform in the literature. Summarise the key ideas in a few sentences. " ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.11", "language": "julia", "name": "julia-1.11" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.11.6" } }, "nbformat": 4, "nbformat_minor": 5 }