{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 3.21 (Diagonals and Multiple Laurent Expansions)\n", "Compute and plot points on the contour of an amoeba, and examine terms in two Laurent expansions. \n", "*Requirements: None*" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[-x*y^3 - x - y + 1, (3*x*y^2 + 1)*t*y - (y^3 + 1)*x]" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Define the polynomial under consideration and the system defining contour points\n", "var('x,y,t')\n", "R = 1-x-y-x*y^3\n", "sys = [R,diff(R,x)*x-t*diff(R,y)*y]\n", "sys" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "At a point of the contour y and t are roots of\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}y^{4} t + \\left(-\\frac{1}{2}\\right) y^{4} + \\left(-\\frac{3}{2}\\right) y^{3} t + \\frac{1}{2} y^{3} + \\left(-\\frac{1}{2}\\right) y t + \\left(-\\frac{1}{2}\\right) y + \\frac{1}{2}\n", "\\end{math}" ], "text/plain": [ "y^4*t + (-1/2)*y^4 + (-3/2)*y^3*t + 1/2*y^3 + (-1/2)*y*t + (-1/2)*y + 1/2" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "and x = -1/3*((2*t - 1)*y - 3*t + 1)/t\n" ] } ], "source": [ "# Parametrize the contour by Pyt(y,t) = 0 and x = R(y,t), \n", "# where Pyt is a polynomial and R is a rational function\n", "J = PolynomialRing(QQbar,3,'x,y,t',order='lex')\n", "GB = J.ideal([J(k) for k in sys]).groebner_basis()\n", "Pyt = GB[-1]\n", "X = SR(GB[-2]).solve(x)[0].rhs()\n", "\n", "print(\"At a point of the contour y and t are roots of\")\n", "show(Pyt)\n", "print(\"and x = {}\".format(X))" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Determine the \"bad\" values of t, which cause the denominator of x to be zero\n", "bad = [k.rhs() for k in X.denominator().solve(t,multiplicities=false)]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Find points on the contour when t lies between -5 and 5\n", "pt = plot([])\n", "for tt in range(-500,500):\n", " T = tt/100\n", " if T not in bad:\n", " ys = CC['y'](Pyt.subs(t=T)).roots(multiplicities=false)\n", " for k in ys:\n", " xs = X.subs(y=k,t=T)\n", " pt += point([log(abs(xs)),log(abs(k))])" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAGECAYAAADEN3+HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAA9hAAAPYQGoP6dpAABgUklEQVR4nO3dd3hTZfsH8G/aQlukLZTSUqAMkVVARGZxgSgKyI+hKA4EB8pSEVHGq4Czr6AoSh0ogoi4XqYyBJWCCMiQsilDkDLKpqWMzuf3x21ITnMy2iYnSfv9XFcumnNOkjtpGu48z33ux6SUAhERERHZF+DtAIiIiIh8HRMmIiIiIieYMBERERE5wYSJiIiIyAkmTEREREROMGEiIiIicoIJExEREZETTJiIyKeYRLjJZDJ5OxYiIrOgIhzLDpdE5HEZGRmIiIhARkaGt0MhorLD6Rc0jjAREREROcGEiYiIiMgJJkxERERETjBhIiK7EhMT0bp1a4SFhSE6Oho9e/ZEamqq09utWrUKLVu2REhICK699lp88sknBkRLROQ5TJiIyK5Vq1Zh6NChWL9+PVasWIG8vDx07twZFy9etHubgwcPomvXrrjllluwZcsWjB07Fs8++yzmzp1rYORERO5lUsrlk994lhxRGXfq1ClER0dj1apVuPXWW3WPGTVqFBYtWoTdu3df3TZo0CBs3boV69atszk+Ozsb2dnZV69nZmYiLi4OGRkZCA8Pd/+TICKyxbPkiMh9zKf6R0ZG2j1m3bp16Ny5s2bbXXfdhU2bNiE3N9fm+MTERERERFy9xMXFuTdoIiI3YMJERC5RSmHEiBG4+eab0bRpU7vHpaenIyYmRrMtJiYGeXl5OH36tM3xY8aMQUZGxtVLWlqa22MnIiqpojSuJKIybNiwYdi2bRvWrFnj9NjCTbrNU/96zbuDg4MRHBzsniCJiDyECRMROfXMM89g0aJFWL16NWrWrOnw2GrVqiE9PV2z7eTJkwgKCkKVKlU8GSYRkcdwSo6I7FJKYdiwYZg3bx5+++031K1b1+ltEhISsGLFCs225cuXo1WrVihXrpynQiUi8igmTERk19ChQzF79mzMmTMHYWFhSE9PR3p6Oi5fvnz1mDFjxuDRRx+9en3QoEH4559/MGLECOzevRtffPEFpk+fjpEjR3rjKRARuQXbChCRXXo1RwAwY8YMDBgwAAAwYMAAHDp0CMnJyVf3r1q1Cs8//zx27tyJ6tWrY9SoURg0aJBLj5mZmXl18V22FSAigzhtK8CEiYh8ChMmIvIC9mEiIiIiKikmTETkE5KSkhAfH4/WrVt7OxQiIhuckiMin8IpOSLyAk7JEREREZUUEyYiIiIiJ5gwERERETnBhImIiIjICSZMRERERE4wYSIiIiJyggkTERERkRNMmIjIJ7BxJRH5MjauJCKfwsaVROQFbFxJREREVFJMmIiIiIicYMJERERE5AQTJiIiIiInmDAREREROcGEiYiIiMgJJkxERERETjBhIiKfwMaVROTL2LiSiHwKG1cSkRewcSURERFRSTFhIiIiInKCCRMRERGRE0yYiIiIiJxgwkRERETkBBMmIiIiIieYMBGRT2AfJiLyZezDREQ+hX2YiMgL2IeJiIiIqKSYMBERERE5wYSJiIiIyAkmTEREREROMGEiIrtWr16N7t27o3r16jCZTFiwYIHD45OTk2EymWwue/bsMSZgIiIPCfJ2AETkuy5evIjmzZvjsccew7333uvy7VJTUzVnuFWtWtUT4RERGYYJExHZ1aVLF3Tp0qXIt4uOjkalSpVcOjY7OxvZ2dlXr2dmZhb58YiIPM2vpuQOHgTuvx/o0gX49VdvR0NE9rRo0QKxsbHo1KkTVq5c6fDYxMREREREXL3ExcUZFCURkev8qnFlw4bA3r3yc2gosGsXUKeOV0MiKjNMJhPmz5+Pnj172j0mNTUVq1evRsuWLZGdnY2vvvoKn3zyCZKTk3Hrrbfq3kZvhCkuLo6NK4nISE4bV/rNlNylS5ZkCQAuXwb27GHCRORLGjZsiIYNG169npCQgLS0NLzzzjt2E6bg4GAEBwcbFSIRUbH4zZRchQpAQoLlepUqwI03ei8eInJNu3btsG/fPm+HQURUIn4zwgQAS5YA77wDZGYCQ4YA0dHejoiInNmyZQtiY2O9HQYRUYm4nDCtWQPcfLMnQ3GuUiXgjTe8GwNRWZKVlYX9+/dfvX7w4EGkpKQgMjIStWrVwpgxY3D06FHMmjULAPD++++jTp06aNKkCXJycjB79mzMnTsXc+fO9dZTICJyC5cTpq5dgf37OapDVJZs2rQJHTt2vHp9xIgRAID+/ftj5syZOH78OA4fPnx1f05ODkaOHImjR48iNDQUTZo0weLFi9G1a1fDYycicieXz5IzmaDWrwfatvVwRERUpmVmZiIiIoJnyRGRkZyeJedy0XdcHBAfX7JoiIiIiPyRywnTmjVAWJjl+s6dwGOPAYMGAUeOeCI0IipLkpKSEB8fj9atW3s7FCIiG8VqXHn2rDSRPH1arjdoIE0kAwPdHyARlS2ckiMiL3DflJy11FRLsgRIQ8mTJ4tzT0RERES+r1gJU4MGQGSk5fp112nPnlNKLkRERESlQbESpipVgN9+Ax58UOqYfvnFMh03ZYp05Q4PB77+2p2hEhEREXmHWxff/ftvGW0y32X58sCpU5I8lUa5uUC5ct6Ogqh0YQ0TEXmBZ2qY7MnI0E7F5eTIormlze+/AzExQEgIMHAgpx+JiIhKO7cmTNdfD3TubLn+yCNAtWpSID5gAHDHHaVjmu6xx6TIvaAA+Pxz4McfvR0REREReZJbF98NDAQWLwaWL5fpuE6dZPsjjwA//yw///YbUKcOcNNN7nxkY50/7/g6ERERlS5uHWECgKAgWXfujjsA078zgikplv1KAVu3uvtRjTV6tOXnxo2BHj28FwtRacHGlUTky9xa9G1P//7Av4uZIzgY2LRJisHHjZNRqYkTgTZtinvv3rF5M5CeDtx6q7YDOhGVDIu+icgLnBZ9G5Iw5eQA770nS6g8/DBQrx5w7bVAVpbsr1IFOHxY2hEQUdnGhImIvMA3EqbCtmwBbrxRu+3QISAiQpKomjXd9UhE5G+YMBGRFxjbVsBVjRpJt3CzG2+UxX1jYoC4OOD+++UMNCIiIiJf4JURJkBOy//kEykSHzIEqFtXe7bZ4sVSPE5EZQtHmIjIC5yOMLm1rUBRREdL0bdZXp52//79svRKuXJy3HXXGRsfERERkZnXRpgK++QTYOhQmYq77TZg505peAkAtWoBe/fKGXZEVLpxhImIvMA3i77tOXwYOHdOzqor3GZg1ixJoDp1ko7iZcmJE8BPP0nX9G7dvB0NkWcxYSIiL/CvhMksMxNo2FD6HAFAVJRltCk4GFi92v/6NhXXyZNAy5bSkgEAXnoJePtt78ZE5AlJSUlISkpCfn4+9u7dy4SJiIzknwkTIFNwkyZJUfjGjdIo0mzAAGkY2bo10LSpkVEZb+ZMWbvOrHJl4OxZr4VD5HEcYSIiL/DNtgKuaNAA+Owz4OOPtS0IAOCrr4DHH5eRl19/9U58RqleXXs9NtY7cRAREZVlPpswWZsyRWp3ateWs+Xy82V7Tg4wYQLw3/8C27d7NUSP6dwZePlloGpVoFkzYM4cb0dERERU9vjslJw9zz0HfPCB7fbQUGD9+rJXEE5U2nBKjoi8wH+n5Ox59VXg7ruB8HCgYkXL9suXpW/TI48Af//tvfiIiIio9PG7ESZrPXoAixbZbq9bF9i6FQgLMz4mIioZjjARkReUvhEma59+KrVNcXHa7QcPyghU9+5S50RERERUEn6dMFWrJg0dU1KkKLqwn34CevcGfv7Z8NCIiIioFPHrhMksMlKaWT71lPQpsrZ4sdQ8ffEF4PrsY9mgFPDMM0ClSkDz5sDu3d6OiMqypKQkxMfHo3Xr1t4OhYjIhl/XMOn5+mugf39L6wFrDRsCK1bYTuGVVd99B/Tta7nerh2wbp334iECWMNERF5RumuY9Dz8MLBrl/xbWGoq0L49MHWq8XH5ohMnHF8nIiIiUeoSJkA6g0+fDgwebDtFd+SITEMNGACcOuWV8HxGr15ATIzl+qBB3ouFiIjIl5W6KbnCtm0DOnQAzp2z3RcbCyxZAtxwg9FR+Y7jx4Hly6WLeocO3o6GiFNyROQVZW9KrrDrrwf27AH69bPdd/w4cOONwPvvGx6Wz4iNlZovJkukZ/Xq1ejevTuqV68Ok8mEBQsWOL3NqlWr0LJlS4SEhODaa6/FJ5984vlAiYg8rNQnTAAQHQ3MnCmJUYUK2n1KAc8/L2u2ZWd7Izoi33Xx4kU0b94cU10s/Dt48CC6du2KW265BVu2bMHYsWPx7LPPYu7cuR6OlIjIs0r9lFxhGzcC998PHDpku69KFeCXX8r2FF1JFRQAWVnSOJRKF5PJhPnz56Nnz552jxk1ahQWLVqE3VY9KgYNGoStW7dinZ1TMLOzs5Ft9W0lMzMTcXFxnJIjIiNxSq6w1q2BffsAvc/8M2eAli2Bjz4yPKxSYcMGmeKLiAC6dOGIXVm0bt06dO7cWbPtrrvuwqZNm5Cbm6t7m8TERERERFy9xLHvBxH5oDKXMAFAUBAwf74srRJQ6BUoKACGDpX/8LmsStE88wxw8qT8vGwZMGOGd+Mh46WnpyPG+tRLADExMcjLy8Pp06d1bzNmzBhkZGRcvaSlpRkRKhFRkZTJhMnsqaekUWNUlO2+ZcuAmjXZ/booLl50fJ3KBpNJO7JtnvYvvN0sODgY4eHhmgsRka8p0wkTALRpY2loWdipU3KW3ezZxsflj155RUbvAKBePeDRR70bDxmvWrVqSE9P12w7efIkgoKCUKVKFS9FRURUcmU+YQJkLbrff5epuMLy8qQlgd4+0nrgARmR+/VXYMsW/QWRqXRLSEjAihUrNNuWL1+OVq1aoVy5cl6Kioio5Jgw/SsgQJZMmTXLMkpi7aOPgIQEwE7dKv3ruuuA228HwsK8HQm5Q1ZWFlJSUpCSkgJA2gakpKTg8OHDAKT+6FGrocRBgwbhn3/+wYgRI7B792588cUXmD59OkaOHOmN8ImI3KbMtRVwxfbtwK23AufP2+6rUgVISZH6JvKMbduAnTtlmrR2bW9HU7YlJyejY8eONtv79++PmTNnYsCAATh06BCSk5Ov7lu1ahWef/557Ny5E9WrV8eoUaMwqAjr7rDTNxF5gdO2AkyY7MjKAtq2lYV8CytXDli9GmjXzvi4Srt586RPVn6+9HJaswZo1szbUZGRmDARkRewD1NxVawI7NgB9O5tuy83V6bnpk83Pq7SbupUSZYAIDOTrQmIiMg3MGFywGQC5s4Fxo3T3//kk8BzzxkbU2lXuFBcr+UDlU5JSUmIj49H69atvR0KEZENTsm56NtvgYcekrXnCuvSBViyxPiYSqMjR4BevaSO6e675XUPDfV2VGQkTskRkRewhsmdNm2Sqbi8PNt9zZrJf/JEVDJMmIjIC1jD5E6tWsk6dHojHtu3AzVqyNIq5D2rVgFjxgBz5ng7EiIiKk10Og6RI3XqyHppcXG2bQeOHQMqV5ZFfPV6OZFn/fYbcOedlqT12DGA7X+IiMgdOMJUDBUryrIp1arZ7svMlKaNbHBpvJ9+0o7wLVzovViIiKh0YcJUTEFBMoKh11jxyhXgmmuAnBzj4yrLGjd2fJ2IiKi4OHFUAiYTcOiQLAdy4IB2X26ujERdusTpOaM8+SSQlgYsWwY0bQq8+663IyIiotKCZ8m5SYMGUhBeWLlyQHa2JFfku375RZLezp2BunW9HU3ZxrPkiMgLeJacUfbu1Z+ey80FKlQwPh5y3bvvSrH4oEHAjTfK75KMx8aVROTLOMLkZjExchZdYRER+ov5kvc1aaJdMzAxERg92nvxlHUcYSIiL+AIk9FOnJCC78IyMqQlAfmeGjUcXyciImLC5AEXLgABOq/sP//Ich/kW6ZNA9q3l9HBZ54BHnnE2xEREZGvYcLkASYTcPmy/r6ffwY++MDYeMixOnWAP/4A0tPld+NKgf6XX8pUXvv2QEqKpyMkIiJvYw2TBx04IC0H9Jw4AURHGxsPuceOHUDz5pYmmXFxwOHD3o2pNGENExF5AWuYvKlePeDzz/X3xcQYGwu5z6FD2o7iR46wSSkRUWnHhMnDnngCuOUW/X1sN+Cf2reXUSWzXr2A8uW9Fw8REXkeEyYDrF4NBAbabr98GXj1VePjoZKJjATWrwcmTgQ+/RT49lvHx+fkAC+9BHTsKL9v69EpIiLyD6xhMpC9YmLXfwXkj15+GXjzTcv1998HnnvOa+H4rKSkJCQlJSE/Px979+5lDRMRGYk1TL7ku+/0t3PZFM8YNEjOgBs+3LtxbN3q+DqJoUOHYteuXdi4caO3QyEissGEyUD33w8EB+vvW77c2FhKu1tvlemyf/4BpkyRpHTYMO/E0rmz4+t68vI8EwsRERUPp+S8gFNznhcYqF8rFBwMPPggMGOGsfHMmgVs3Ah06ADce6/94/76C+jRAzh6VBLs2bOBoCDDwvQJbCtARF7gdK6HCZMXDBokox+FPfyw/AdJJRcaCly5Yn9/QABw003A0KHAAw8YF5czrVsDmzZZrs+cCfTv77VwvIIJExF5AWuYfNEnn+hv//prY+MozZKTHdeGFRQAv/8O9O0L/Oc/hoXlVOEFmrlgMxGRb2DC5CXz5ulvf+cdY+Mordq2laTo++9ltMmRxEQZcQoMBG6/3Zj47Bk1ypLo1akjCZ09M2YAAwfKMi1ERORZnJLzItYyGefKFaB7d+CXX5wfO2QIMH68TI21ayd9l4y0bZsstXLTTUDlyvrHfPIJMHiw5fq0aZI8lQackiMiL+CUnC976CFvR1B2hIQAK1bIOnB9+jjusr5ypSxd060bEBUFrFljXJwAcP31wD332E+WAODXX7XXf/vNszEREZV1TJi8yF7NEvsyeU6TJjJNd/GijOTVq2d7TGqq5WelZHmbdetkhGrQIN845b9lS8fXzfLzpXD83XflzDtflpSUhPj4eLRu3drboRAR2eCUnJdxWs67rlwBHn1UCsAjIyUhevZZ7TGRkVJ8bW5T0KgRsHu34aFq5OcDb70FrF0rU3djx0odVmGPPgp89ZX8XL06kJICVK1qaKhFxik5IvICthXwddOmAU8/bbudCZN35OUB5cppt7VpA2zYoN22a5dMm509K1N806YZF6OrlJJFga1HxH74AbjvPu/F5AomTETkBaxh8nVPPeXtCMhaUJCs9VaunIz+9e6tn2C0bQv8/beMPH32GTB1qmw/eVIW2/UFJpOcaWfvurXZs4GRI4Fly4yIjIjI/3CEyQfoTctxhMl37NwJNG3q+JhevaSgfN8+uf7667Lorrft2CFJ+dmzsqbeoEG2x7z3HjBihPxsMgGLFsnombdwhImIvIAjTP6gcHLEZMm3NGkCxMZarl9zDVCzpvaYpUstyRIgbQkA6ZH0+ONyhp43NG0qdU579ugnSwDw00+Wn5WyHWX66KOPULduXYSEhKBly5b4/fff7T5ecnIyTCaTzWXPnj1ueDZERN7DhMlHKGW5kO85dAh4/nmpNzt0yLYlROFlWAoKgOeeAwYMkAaTnTv77rI3hUfPmjSx/Pzdd99h+PDhGDbsP0hI2IJDh27BHXd0weHDhx3eZ2pqKo4fP371Ur9+fQ9ETkRkHE7JERXDwoVAz56W60FBMp2VmyvXAwJkJOrCBcsxN9wgU3VPPSWF2K+/rl/wb7TLl6V+KSUFuPNOGR0zTxO3bdsWN954I9LSPsbixeZbNEbPnj0xf36izX0lJyejY8eOOHfuHCpVqlSseDglR0RewCk5f2EyWS67dnk7GnKmRw/pCB4aCkREAP/7n7YZZkGBbfF3VJQkWcePA6dOyRTZzp1y7J9/AseOGfoUrgoNBZKSgD/+ACZMsCRLOTk52Lx5Mzp37oz9+61v0Rlbt661uZ/cXMvCwS1atEBsbCw6deqElStXOnz87OxsZGZmai5ERL6GCZMPKFz0bT0lQr4rKQm4dEnOlOvRw/b32KiRZZvJBDz8sPRPsvbbb0B0tCzBUqMG8OqrhoTuktOnTyM/Px8xMTHo08eyPSgoBgUF6Zpj8/KAu+4CXnwxFsA0tG07F/PmzUPDhg3RqVMnrF692u7jJCYmIiIi4uolLi7OQ8+IiKj4OCXnA3iWXOnw8ccy6gRI/6M2bbTLqtSuDZw4Yal3CgiQxXXnzLEcExQkIzUvvACsXg3cfbdM3XnDsWPHUKNGDaxduxYJCQn45hvgwAHgxIk3sWLFV5pC7t9/B269VXv7zEwgLAzo3r07TCYTFi1apPs42dnZyM7OtrpdJuLi4jglR0RG4pScr7O3pAX5n8GDpSB8wQKZcrOuXwJkJKpaNUmQg4OBH3+Uf60pJYXikyfL9NYbb0ixuZm527gRoqKiEBgYiPR0GU168EFplRAQcBIxMTGaYwvnNSEhlufWrl077Nu3D6mpkkRWr65tuRAcHIzw8HDNhYjI1zBh8rK//vJ2BOROtWvL9Fx4uKW3kVlIiCRUSgHZ2cDEiXK55hrLMYMHAz//rL3dokXA9Oky+hQYKFO2RiRO5cuXR8uWLbGiUE+EFStWoH379pptzZtLsXhgoNRyffmljLIBwJYtWxAbG4t+/YCNG6WG6803te0MiIh8XZC3AyjL7HWErljR2DjIMx59VNZt+/Zbma4y92Yy27FDEixz/c8bbwCtWgFbtgDpViVCDRrIVJ+5/mnXLuA//5G2BUOGyBTe++/rLyRcUiNGjEC/fv3QqlUrJCQkYNq0aTh8+DAG/dvUacyYMTh69ChmzZqFCROAsLD3UbduHTRr1gQ7d+Zg9uzZmDt3LubOnYthw7T3nZamvZ6XB0yaBGze7P7nQURUUkyYvKjwdIxZ4akc8l9dusgFALZuBT780LLvzBnLzz//DLz0kowcLVsG3HILsH8/EB8PzJ+vPQMPkISqfn0gK0uur1gh2ypVkvvQW4i3OB544AGcOXMGr732Go4fP46mTZtiyZIlqF27NgDg+PHjmp5M+fk5eOmlkTh69ChCQ0PRpEkTLF68GF27dsVff8nIEiCJZOFu4uPGAYlWnQrmzwf693fP8yAiKikWfXvJPffAqq+NFgu+S6/EREluEhKAt97S7itXTkaLIiOB7dul1sesa1fpJm4+bs4caM5cA6Qx5pdfyv2bTMDo0baP4W0//ggcOSLv/8Inw912mxS6A5kAIvD00xn45BNtPdM338hrc9ddcjwRkZs4LfpmwuQFStkfAWCyVHZUrQqcPq2/r2VLoF8/oFMnSyfuzz+X5VeGDJG6p+ho7ftl0iTgxRe193P5soxevfWW9IH65hvbAm1fMXo08PbbgDlhmjEjAwMGWIJ9911psAnI38/y5fL6EBG5gdOEiVNyXmAvWWrVytg4yLsOHJAi76wsOS3/3DnLvs2bLbU8X34p9VBPPqm9/dSplmm8UaOkKLywX3/VdiSPj5cRnoICKcCuUcN2XTxvqV49CVWrJuHy5XxkZQG9e2v3L1xo+bmgQEZorROmggJ5rQ4fBu691/mCyURERcERJoPp9Vwy4+hS2fXqq9JlW09oKFC3riQQjnoyZWXJNJ65Bi4+XqatPv5Ye9zly1JsfvKkXH/xRTlbz1fYWxrlqaeAzz6zHPfpp7LN7LnngA8+kJ+vuUYSzoYNDQqaiPwd+zD5EiZLZM/48TJiMnIkcO212n2XL8uZcW+8IUlV4YV+zSpWlNGjl1+W6audO6V43FqFCnKGnTlZAuRYQO6/alWgTh1ZqsXXvPuu9Khq1Uqe48CB2v3z5ll+vnhRpuysnT4NjBkjr/E//3g8XCIqZTjCZIAdO4Bmzezvz8mRQl4iQBbBbddOejXZ06ULsGSJa/f37LMyVVWxoqx598032rP1TCZJ1rp2tWwLDZVlXy5dkmQuIEBGwUJCivWUiqS4i+926iRLzZgtWybF4YC0ZGjRQgrGASk437HDd+u5iMhwLPr2tocf1i59Udjff8t0C5G1ggKpcerVS0aK9DRoIEunTJlStPs+f16m5Mxr3D7xhJyZN2mS9riLFyWxOHtWrkdHS9NJcw1eTo6lOaU7FTdhOnpUasIOH5Z2BNYd0o8eta3VWr8eaNvWcv2HH4Bp04CYGOCdd6QrOxGVGUyYvCUnR9bRstecEpDGfb5ScEu+6exZ4M47pSeTOcEprFMnWQi4KPU6V65IB/G6dYHWraUAvE0by/7KlYEZM7QF4wCwcqVMDz7zjCR1jRpJQueuvk9A8RMmR3Jy5LkeOybXw8PlNa1aVa5v3izP39xB/dZbgVWrbO9HKcdT60Tkt1jD5A1Dh0pTSkfJUl4ekyVyLjJS/jPPyJARJT2//iqJy4ABrt9vSAhw//2SLAHy7w8/yNTxLbfIdFXhPkmAbBs+3JJY7NkDvPKK/Pzgg0BsLNChg/06K28pX15qmv7v/2RUbulSS7IEANu2aZebSUnR3v70aeDmm2XqvEMH7RmNRFQ2cITJjX77Dejc2bKEhZ6QECniJSqOqVNlNGnPHv39FStKrVHhdeyKa+BAWccOAIYNk7PQAgK0JymYi6+tz2BLSADWrpUvDaNHy3v+9delF5QznhhhcmbfPuCGG6RmCwDuu08SSLMhQ7RnGw4fDrz3nvY+/vxTFl3u2FG7PiAR+QVOyRlh+3b5Vp6R4fi4wYOBjz4yJiYqvXJyZJRj82b7o5itWslp9zfe6P7H79bNUnBevjywdy/QvbuloBqQKb2zZ6XNwfHjsu2aa2T5looVZTQnJUXOCKxUSXv/3kiYAGDTJuCrr6SGacQIbYH7Aw8A339vud6vHzBrluX6G29YRtqaNQP++EOm5InIb3BKzpM2bZJC2Ouvd5wsBQdLbxwmS+QO5cvL6M2FC/YLkzdtkqRp9GiZ/nWnxYuBmTMlQThwQArIO3bUHtO8udQ2mZMlQIrIv/tOEqnISOlmHhlp+btISkpCfHw8WpvnCQ3WqpUU0I8da3s24ODBlm0VKgD/rj181X//a/l5+3bLMjZmmzfLUjaPPCInehCR/+EIUzHMmCFD8vaKcK1Nnqw9W4fInQoKgAULgIcest+GoEoVGfHwdBPHwYOlTqhpU2DuXPn7qFJFe8zKlZIgWU93VaggydTp03LWWlpaJnJzI3D0aAaqV/ed8/737QO2bJFEr1497b6aNeVMPLOff5bpeUCeV/36cnYiIMXne/dqO7MXFAAHD8rrVXjEjYgMwREmdzl3TpanKFcOePxx58nSPfdInQeTJfKkgADpAJ6aan9dtTNnpCi88NIq7vbxxzLitHChJAORkdJoMzBQzizr31+mEgtPI5qLrbt0kdGX3Fy53revZf/DDwONG0s3b2+pX18K5QsnS4BMz0VFyXMdNsySLAGSHJmTJUASI+s1BK9ckd/dddfJFOZPP3nsKRBRCXCEyQGl5INwwgTg0CHXbnPbbbJaPBtRkjd88w3w9NOW5VEKCw2VFgJNmhgbl7UDB+TxzSNiY8bI4sCWeidZfLdOnQwcPBiuqZkCtAXXycnA7NlSQ9i/v6FPQ1denu2afufOScJq7q7euLFM2wUGyvWZM4HHHrMcf911Mppl7fvvJRFt1EimWfn5QuR2HGEqjuXLgZtukg++AQOcJ0smk9Qn5ObKBzg/zMhbHnxQRj8feEB//+XLMmVWuL+SkerVk35IX3whydtbb8n2e+/VHtenj/xbeJmWZcvk3//9T2qnpk+Xv9OhQ7XHrV4NrFnj9vAd0lsAuXJl6ek0cKD0r/r1V0uyBGjbGQC2Z9kuXSq/zzlzgHHjgBdesH2My5elbs26ZoyI3Ewp5eql1Lp4Uanvv1fq9tuVKldOKRlbcn655hqlEhO9HT2RvowMpSpXtv/+DQlRasMGb0ep9fbbSnXsmKEAqIyMDKWUUgkJ2rgfekiObdVKuz0iwnI/LVtatrdpY/zzKIqLFy3PMThYqblztfv/8x/t82zRQrv/zBmlGje2/E5/+sm42IlKEad5UJkdYTp7VgqymzeXIsv775c+Sub6CUfatJGzXrKyZHicyBeFh8s00LPP6u+/ckXey/ZGo7zhpZekiN3aL79IX6cqVaTO6auvZHtMjPY482n8q1fL36fZhg3AunXyc1aWFGibTDI9mZzsiWdRNBUqSMzbt8uyLr17a/e3a6e9npCgvf7FF8Du3fLzlSuyMHFh//ufnAXYoYO2/QMRua7M1DCdOCFD1t9/L0Pcp04V7fa1akmtxVNPuXcZCCIjbNwoHa7N68IVFhIifx/erG0yc7UPU2amTC+mpUmy9NtvkhSsWSM1TdbM68b17g3Mn2/ZHh0tnw2AFKPfdJO0Q6hVS+7HlUabRpg9W2qYGjaUdg7BwZZ9U6ZIXZdZmzbaacx9+4D4eEt7ibg4ScwK++47KUjv0UPqrIjKGOeLHrkyDKX8cEruyBGlli5VavBgpZo1UyogwPWpNvOlTh2l3nhDqfPnvf1siEouP1+pvn0dv+f79/d2lEplZGin5IqjTRvLc0pIsGy/+WbbaXWz//s/x1NfH3yg1LXXynTfP/8UOzS3y8pS6qabJObKlZVas0a7/+efbX/PFy9qjxk71rKvYkWldu/Wf6z9+5U6eNAjT4PI25zmQaUmYfrjD6XmzVOqXz+l4uOLVotkvgQGKnXDDfLBWPgDhai0WLVKEgV7fwdRUUqlpRkf19SpU1Xjxo1VgwYNSpwwKaXU2rVKrV+v3bZsmVImk+W5Dh5s2deoke3rYLZ4sXZflSq2jzdypFK33qrURx+VKOxiyc+XL4mXLtnuO3tWqZo1LbHfdZftMfXqaZ/fpEm2xwwdatn/n/+4/zkQeVnpTJj+/FOp335T6p135JtkgwZFT44AGXWKjpZv1T//rFR2trefGZExcnOVuu02x38fU6d6JzZ3jDA5smGDUs8+q9ScOdrt48drn/+DD1r2DRpk+/pY69VLu++992wf98svlXrsMaWWL3f3M3IuLU2p119X6v33lbp82Xb/nXdq4583T7t/927b53/8uO39fPONUrVrK3XddTLCT+RHnOZBPl3DdPSo1BpVry71Q4cOyem4K1YU/z6jo2Vl9l69gLvukgJQorJqzhxZF63wqe1m7dtLl3AjeWstOQCYNElqhdq3ByZOtGxftEhqe8zMa+WZVaqkXR6pRQvgr78s159/Hnj/fcv1r76SZVKsXbki6/9VrixNco105Ii0PTh4ULrGjxun3Z+aKj2grKWnawvvjxyRtQHNJ85Yrx1o7dIl6W134IAscvzgg25/OkTF4fuL7xYUyLIIYWHArl1y5lpIiBQdDh8uhYqVK0vzt+KoWVN6tXTqJAWd117Lom0ia+fPy3If9tY4K19e/jb1Olx7gjcTJkcmTgQ+/FA+j+bP174eDRtKR2+zvn2liahZZKT2M6x5c1l82CwnR74Ynjkj11u1kkJ9PStWyJl1N91U4qdUJNZJ3/jxkvRY27xZ4raWlmb7pfSxx6RZp9ny5cCdd9o+Xk6O9Ng6d06S+ri4Ej4BIsd8r+h7wQLptbJ1qwyNR0fL8O4dd0hdgHm4NyioaNNrQUEyxXbPPUq9+KJSq1fL3D0Ruebllx3/jb3zjjFxeHpKzhPS0pSqW1f6ILVtK1Oe1mrX1r6WnTtr9ycl2b7eJ07YPo51+UHHjvbjyc2Vz9gLF0r81DQOHZJaKT3Z2Uq1bq2tlSoosD2ucK3YG2/o31/v3pZjYmOVOnXKflwHDkgvvT17iv6ciP5lbA3T/v2SBJk/LJYtU6p5cymkXrFCqYkTLX8AISFKNWlS9LqjwEClypeXws0BA5QaM0apadPkMXNy3PbCEZVJ27c7PmHC+owzT/HHhMmZrVuVCg+X17BGDaXOndPunzPH9rUunOx8953tMX//bftY//xjKeoPCJDbOZKd7b6TXLKylJo+XalZs+zXhD71lCV+k0mp5GTbY3JztcX5gFILF+rf359/Wp5vcLD8X+Msxg0blDp5smjPjUo9zyRM27crNXy4Uq+9Jm8+peTMMvMbvEMHpY4dU6pCBcubvWJFpZo21f4BVK2qvR4Sov1GERpq+Tb29ttKDRsmI0e5ufqFi0RUcrm5cvq8vaTpmms8exZpaUyYXGHd0dz67D2zH36w/V0cOmR73O23a4/RO6PPbOBAy3F33ulanPn5rh1nz5UrSk2YIGc0L1hg/zjr92BAgPy/o+fxx7XPt1cv+/d55IiMBAJKhYXJGaPO5ObKCCK/kJd67kuY1q2Te0xLkyUIzG/Ou++W7YVPU/7wQ9s/7jvu0F4fNMgy9Vavnpw9cv/9Sj36qPT6yMxU6vBh/WFdIvKswYMdj/b+/rtnHresJkxKyTInjpJR6y+d5s/ewsw9mcwX6yVjrB0/bvs7XbLE/mNv3mz5Elyliv5Zcu60Y4ecydm8uYxY2fPii9rn8MQT9o8dM0Z77O23O44hLU2p+vXl2Nq1ZRbFVebBBPIb7kuYqlWTeyz8LScwUBKawmtW/fijdj2n1q0lu+/WTU45HTlSbrdvn1K//irrXhGRb1m71nHS9Npr7n/MspwwuWLNGqW2bLG/f+NG+Vw2/44mTtQ/bs8e29/nzJn279c8MmO+3HKL81j79ZOZh4AApUaPdn58cZw/L4lPYKBS7doplZ5u/9hx47TPwV7SaVb4S0O/fs7jSU+XMhRAmiYfPVq051PSETwqNvclTAEBMiS5c6e2IPuGG+SRfvjBMqX20EOSDGVkSN+PKVNktIiI/M+FCzKlbi9pat/ePY/j7saVZdnRo0p9+qnjxEopy6K95hKJwsXq1gqXUDRr5vi+f//d9r2iV8iu5777JNEKDFTqrbdcu40rMxFnzliSmZgY56/Pk09q43/gAeePMWSI9jYDB7oUvlJKeoGVL69UZCQXUfYC9yVM1ksmLFyoVKdOSvXpo10iICvL9T8IIvIv5v9o9C7WXbFLiiNMxvroIxmFcta4d8IE7e+8cOPPwr780vZ9snWr83gWLrS9XVFrVlu0kNsFByv1v/9p9+XlyVTblSvO72fPHkmszNOQKSnOb/Poo9rY+/Z1LeYNG7S3Cw8v+mjT99/LbE/FilIWU1JlrBzGfQkThwmJ6Pnn7SdNQUHuORmDCZPv+vFHeQ8UXq9Oz8WLljMDAaWqV3ftMd57z/a9VZRprZde0t42NNT12xb2/vuWdUhdPUN00ybL8w4Ls12ex5r1iF7hNf8CAoq2+sSFC5Igmm9vMkm7heIoKFDq6afljNnatWWa150KCiRx9THOZ9pc7ejEZo9ENHkysGCB/r68PCA0FDh82NCQyED33CPvAVeaZlaoAPzzjzQgHj1auoi74tFH5X1kVru2NPV0VVqa9np2tuu3LWzECEsX/HXrgI8/dn6bli2B3buBpUul4WvbtrbH/O9/QLlycqlVS7q833or0KaN5Zhhw6RprKsuXtQ+V6W03ejz8lz/25w/X7rO5+bK7/Cxx1yPw5mFC6UzfkiIbUd5T8nPlwa9JcU0iIiKpEcPx//51a4t/7kQVaoEvPcekJjo+n/+kZGyDNbgwZJo7d9ftMd86SXtF/zbby/a7c0KCmyXDHI14aheHbj7bvtLbz3+uCQwgCR4gwZJApGcLAnFypXAlClFizcmRjrMm3XoANxwg/z8zTdAcLD8bUZGahMpPeaO8/auF1d+viwJlJkpz//116VDfGoq8Oef7nmMwjZsAGJjpUP/3XeXLIH2+tIoROSfzCNK5g/+wmbPBh5+uOj366tLo5D/2LED+OADWbLmhReKfz8tW1rWBAwOluWDijLaZU9wsCz9YtatG/DTTyW/34ICYNkyue+uXS1JauG1Dnv2lFEke06dkpEx8xejiROBF18seXxXrsjIo3Xa0bGjJIgAUL++dokhd2jbVpIms48+kmRch9OlUTjCRETFEhQkQ/aVKunvf+QRGV0gMlrTpsC0aSVLlgAZ/UhMlGnF/fvdkywBssCxWUCA7bp8xRUQIIlSz57aET3zgshmWVmO76dqVWDTJpk6XLvWPckSIKNozz1nud62rSVZAoB9+2Qq0J0uXdJev3ix+PfFESYiKrG6dWUaRc8rrwCvveb6fXGEicqC77+XBZgfe0xGVjxp1CgZJQKAwEBJhszTdd6wZo1My7VtC0RFafd98AHwzDPue6zvv5eR7rw8WTB73TpJCHU4HWFyKWEymUymjIyMAqcHElGZdeutwNat+vueftrygV1YdnY2sq0KCy5cuID4+HikpaUxYSJyk99/lwStZ08gLs7b0Vj06QMsXy4/x8ZKoby5Bi0vT+q9du2SkbOifPGydvCg1Iq1aAGEhekfExEREQHggnKQFLmaMIUDyHB6IBEREZF/ilBKZdrb6fYRpszMTMTFxbnl22Hr1q2xcePGEt2HO+/HnffljvvxxdfanfflS/fD19r1+6lffyMWL9bf/+KLwMsva7cVHmE6fvw42rRpg127dqFGjRpuicnXXiNfuR93vq/dFVNpvZ8TJ4Bt2zJx330le70tIzKt8fbbGzFoUInCsnludetqz6JLTASGDCn6/Vg7e1bqqGJiihZTzZrAhQuW7VFRcqZdjRrA558DjRvbv4+ivrddGWEKciV4R3dgT3h4eIn/AAMDA93yR+yu+3HnfbkzJl96rd15X752PwBfa1fu56efwnHvvcC8ebb7J00C6tRx7QM4LCzM555babwfwD3va8D3npu772fFCmDRIqBLF5kictXUqdq6nIsXw1GzZtHjmj3bMn0FBGL06HC89FKR70aj8GtUuF/R5s0yjfXLL3JK/l13Sf8oZ/djNnkyMHKknBk3YgTw7ruux1Svnkwjmp0+Lf+eOyd9qlzJYV19bzsaWTLz6bPkhg4d6lP34877cmdM7lCaXyO+1sbfz9y5cqq0/jHWH/qe56uvka/cjzv52nNzdD9ZWfIebd5cRiuc3c/HHwOdO0vy062ba//xm40erb3+/POu39bagQOaqKCU7Vlgrvj1V+lX9dVXtq9R4QL0hx4C+veX5969u/Qy0mslovdanz8vo8rmIZfJk4GdO53HZ76vVauAVq2AKlXk92Tt+HHn9+N2rrQD//fiEi5rYBy+1sbha108N99sfymVtDT926SlpSkAKs3eAeQ2/vi+PnTI8Zqlzz4ry6FERNiuI2etXj3t+3HePMePe+212uNr1nQ95goVzLeT1/vOO4v+euflKXX8uHb5k4YNLfvXrVPqjjuUuvNOpf780/79rFhhWe4FUOq//9Xuv3xZqR49lLrxRqW++EKpkydt/3Y3bHAt5tOnbW/rbMFje44cUapaNcv9vPaa4+OL8d5231py1kHfd598EH71le0jXrlyRY0fP15dcWVlQyoRvtbG4WtdfA0a2E+a9NaoPHnypAKgTp48aXywZYyvva9zc5UaM0aSnjNnbPe3bGl57zz6qO3+JUu076/AQPuPZZ00ALKYvCPWjw0o1aSJ688rMdF8uysqMHC82r1b+3pnZCj1yy9K7d1re9uCAqUGDZLnEh0tC+w+/bRSY8da/n7On1eqUiVLbJGRSmVm6scyfLj2edx8s+PYs7Js16jTi9OesWMd/86K4uhRpaZPV2r5cufHFuO97f6EKSFB+8I5ymSVkjfCE08o1b69bSZLRGWD9Yd54Uth/jjqQa7bsUOSjTp1lPr8c+2+uDjL+6JiRVlQ1mzOHNv3TuGk6vXXbY85d04/jqgo7XEffug47kOHlKpcWY6NiFBqz56iPe/331fq//5PEh5r6elK1a0r9xsUpNQ332j3//ijNs66dW3ve8cO2+edmqofx2efaY97+mnnsX/7rSwoHBIiiyMXVWqqUtu3F/12BnN/whQWpn2xZ8xwHMGAAdrjv/7a/rHnzsnxt9yi1AcfFOf5EpGvCgzUT5gKjwIwYfI/S5YoNX68Urt2Wbb17Wv5Yv3SS5btlukpuaxfL9v37LF9b3z6qeV2n35qu//oUW0c+/dr32dVq9qPec8emZarXNm1pKGwzEyZrrK2cqUMKrRvr9Tq1Zbt06dr47ZOmiZN0u5r2lR7n7Nna/dXrmwby5UrMj1nPqZxY6Wys/XjLihQ6tVXlWrXTqknn9Qmpc4UFLh+rB9yf8L04IOWX0pYmFIHDjiOoFUr7S97zBj7x/bpoz124UL9486eVeqxx5Tq0EGyZSLyD/ZGmerXtxzDhMm3WE+bZmQo1bGj1POMHSvbXnzR8nsMCFBq1Sqlfv3V9nd84oTcvvD2V16R+zl3znbfL79YHjs3V1vDkpCgH++GDVLLc++98ngldeyY/D+zeLFl24wZSpUrJ3EMGiTbzp7VDihERFgev0cP7fPq189yXx9/rN3Xrp328c+flwTIvD8xUT/O9HSlXn5ZXk/OZheL+xOm7Gyl3n1X/ki2bnUewfjxtn9M9jRqpH3jvPmm/nE9e2qP+/ln22NycpQaPVr+cN56q9RnxkR+IT/fftLUp89U1bhxY9WgQQMmTF6yZIlSgwfLl9XcXElkzVNF335r+xn9/vsydWa97eabpVi48O93xw55DHsjTErJf/gBATIq1bevbXy5uXLfjoq5i8q6xGXFCqW6d5dam6NHJVmKjbXEOnasxGBd0wMotXat42kx6zqewknPlStKde0q26tVU2rTJtsYMzNlas7VYmsqFvcnTEVVUCDz1C+8IMOVjowYYXlDBQXZr48yz/eaL5Mm2R4zZoz2mKQk22N+/12Sr4cflvlpIvK89evtJ025uRxh8rRduyxJSp8+MsXTpIl8QS080mF9vUIFpcqX12674w6lqlTRbuvWTamLF6XmxbwtLs7y+Nu3y+PVrm1bw+RJBQUySvTtt1LIvGWLnOkGyAjQzp3aROjGG23rfaKjZdCgcMH4ypWS+DRtatl2ww3yxV0p2ffMM0q1aaPUyJHyPi8sM1P/JAgyjHcTpu7du6u4uDgVHBysqlWrph555BF1tPCks5X8fElsRoyQZMaewYMtb8py5ZTavNn2mLvv1r6hn3xSuz8tTalrrrHstz4902z+fKV691bquefsn3HgCw4ePKgef/xxVadOHRUSEqKuvfZaNW7cOJVtbxKbSuSNN95QCQkJKjQ0VEVERHg7HL/00EP2kyYmTO41b55SbdtKUtC9u+V1DgtbpYB7FBCrAKiAgPma30NIiPb3EhRke2p9YqIkIUFBcr1KFUshdkaGJAevvKKfIHhKVpb8e+qUUk89Jc/7p5/k/wBz3C1ayGti/Vys/18BZJRr0SLtthtukPt+4w1tgpiXJ9vPnJEZjcRES7H5W2+9pVq1aqUqVqyoqlatqnr06KH2FLVinFzy0UcfqWbNmqmwsDAVFham2rVrp5YsWeLqzb2bME2ePFmtW7dOHTp0SP3xxx8qISFBJdibeC6C3Fyp1H/2WfuJ1TvvaN/o336r3b9ihe0HtXVStHat9ltE797a2x88KH+A/ftrCx29YenSpWrAgAHq559/VgcOHFALFy5U0dHR6oUXXvBuYKXUuHHj1OTJk9WIESOYMJVAaKh+wlS5MhOm4jh0SGo7R4yw1LyYExn9yxIF/EcBcxUABWgTpnr1tAXUvXtLEtKqlUxTDRumfXxv/LoKCiQhnD5d6mnNo2L162vP6NY74aBOHe31sWO1o2K33SaPMXasjCzdcIP2TK89e5TauNH5qNBdd92lZsyYoXbs2KFSUlJUt27dVK1atVSWObMjt1m0aJFavHixSk1NVampqWrs2LGqXLlyaod5Ptgx7yZMhS1cuFCZTCaVYx6n9LDPPpOCvMKncSolBYiRkZY/jpYttfunTNH+MdWoYdmXna2dFoyJkcI8s82blXr8caWef14ad3nDxIkTVV2980/JbWbMmMGEqYT0/yNnwuTMoUMycvPyy9KyoVw522miol2gateef3V6rVIl6bWzZ4+MEn3xhbefsfQATE+XL8l33CGzCL17W55D4dYVhZPFwgn6m29aEqlataRe6a+/5Iy50aO1n+nuZO4ztspRQS+5TeXKldXnrs39Os2DXFpLzh3Onj2Lr7/+Gu3bt0c5vYVoPODJJ+WiJzoaWL0a+PBDoGJF29b17dsDgYGy0B8A3HKLZd+xY8DBg5brJ04A+/cDLVsCR44AHTsCmf+uSrN2LbB+vfycnQ288w5w+LC0m7/tNvc8Tz0ZGRmIjIz03AMQucGFC7JOFTn200/A9u3AzJnyGXP58r9pThEFBAAF/y6jXrs28PTTwDffyH2/+y5w7722t5k0qUShF9ulS8DHHwMXL8pn5vTpsj04WD5LCyu8BlrlysCpU/JzSIh81o8fL5/Nr7wiS3b07AkcOgQkJMjxsbHAJ5948ElBPpsB8PPZw/Lz8/HDDz/g4sWLSEhIcM+dupJVqRKMML300kuqQoUKCoBq166dOu2tIZdiWLZMTv/8z3+kiNEsO1s7nBsdbfk2UrjJGGA5C+PRRy3brGuvcnKkcH3IEMdnEbpq//79Kjw8XH3GngsexREm9+jcWX+EyWQquyNMS5bIiTI1atgWWhf1UquWnO4eGyvF3r//btuyBYCaP3++V56rUvIZuXSp1EPddZecLWbdQsDVS0SEZTTpxx/l5J/HH3dcE2ukgoIC1b17d3Wzs/baVGzbtm1T11xzjQoMDFQRERFqsXU/CMfcPyU3fvz4f+e77V82btx4NYJTp06p1NRUtXz5cnXTTTeprl27qoJScI7/339LvUC/fpbTZZWS2ibr02abN7fsq1FD+8c9ZYpsf+opbSJlPq30+eeL9lorpdTRo0fVddddp5544gmPPv/Spqjva6WYMLmTXsIEZDhdSaC02LtX/nOvXdtZ3ZHji8kknyFjxkiNkd7ZwXqMTpgyMuTEm/HjZWqtdu3iPd/rrpN6rWrV5EvniRPyRddXz3oeMmSIql27NtdJ9KDs7Gy1b98+tXHjRjV69GgVFRWldu7c6cpNneZBJuX6uK4CgNOnT+P06dMOD6xTpw5CQkJsth85cgRxcXFYu3at+4bIfNDvvwPvvw9ERABvvAFUry7be/QAFi2Sn00mWYn5lluAunVlWNjsnXeAQYOAxo1PIy1NXuv69YF582RI/dtvgbQ04I47gHvvtbzWx44dQ8eOHdG2bVvMnDkTAQEBxj1pP1ec9/XMmTMxfPhwnC88F0DFYjIlAUgCkA9gL4AMAOHFmnryB3PnyufEhg1ATk7x7+eWW4AKFYCxY4Fbby3efZhMJsyfPx89e/YsfiBOnDgh0107dwILFgC5ua7f1noq8bnngNBQoFw5+blKFY+E63bPPPMMFixYgNWrV6Nu3breDqfMuOOOO1CvXj18+umnzg41OTugyDVMUVFRiIqKKurNAMhoFgBk601AlyK33KKteTL78ktgzBiZj3/kEcsxLVpoE6YbbgBSU4G0tCgA8lrv2ye1HlOmAO+9J8fNni2JVJs2wOTJR/H66x3RokVLzJgxg8lSEZXkfU3u0aHDUCQnDwWQCSDi6vb4eGDXLq+F5Vaffgr897/AP/8UrwapYkWgQQOgb19g8GC57qsyM+UL4q+/SoKUlQXk5bl++3LlgOuvl1rSsWOBhg1lW716HgvZI5RSeOaZZzB//nwkJyczWTKYUsptOYfHir43bNiADRs24Oabb0blypXx999/Y9y4cahXr16pHl1ypFIlKWIs7Isv5FvSoUNSDN6pE3DmDBAebikej4mRy9Klltvl5QG//ALMnXsMEyd2AFALycnvYN68U2jdWookc3KqYcAAoHFjjz+9MuPw4cM4e/YsDh8+jPz8fKSkpAAArrvuOlT05f/BfNzKlTLyWtju3cbH4i6XLgGTJ0vB8cmTRb99SIgkCi+/LF+MatVyX2xZWVnYv3//1esHDx5ESkoKIiMjUauYD7Rtm3wmzZkjv7dLl1y/bUCAJIBjxshn3623Ak2bFisMnzJ06FDMmTMHCxcuRFhYGNLT0wEAERERCA0N9XJ0pcvYsWPRpUsXxMXF4cKFC/j222+RnJyMZcuWuecBXJm3U8Uo+t62bZvq2LGjioyMVMHBwapOnTpq0KBB6siRI0W9qzJr7VopgOzWTamUFNlWeL29RYuUuvbaGQ5qb+R027Q0Kfh8/nnpYWVkI7nSpn///rqv9UpnrezJqbQ0bQ2T9Xvdn3z5pSz8WpwapLp1lXrtNVlI1pNWrlyp+z7u37+/y/eRm6vUtm3S3Ne8jEpRLo0bSy+7WbPkvkpBeasNe5/NM5ytXE9F9vjjj6vatWur8uXLq6pVq6pOnTqp5cuXu3pz99cwkXedOweMGAEcOADcfz8wbBgwcCDw+eeWY2bOlG3WNQKTJsm3VPPI5MCBwJtvyrffvDzgmWfc++2VqLhMJvOUnNQwmfl6LVNqKtC7d9GnD8uVA266SabqWreWkRZfVlAA7NkDTJ0qdVhFGTmLiZHR7ueek7YG5ik3Ih/gtIaJCVMpkJUFvPACsGMH0K2bzPe3aAH8O1OEwEDZP3Gi5TY1a8oU4Y4dcr1WLekZ9dFHkpQNGiQfZkRGy8zMRESEbcIE+GbS9MorwNtvF62IOTRUkqtPPvHtOiRrx47JCSmzZknJgCuCgoAaNaRms107+XzSm3Yl8gFMmMqqQ4dkJOrMGRmFqlJFaqPMOnQAkpO1t2nZEti8WX6OiABWrJDRqpwcSbgaNTImdirb/CFhunRJamzMfy+uCAkB+vSRUV1/Ob8gJ0fqrz78UArVXZWQAPTqJY0xw8OdH0/kA5gwkcXnnwNffSWjSYmJUkR6/Ljssy4wN6teXb5VAjKU/uGHwHffAdWqAa+/Lp1xidzNUcIEeDdp2rdPRkrOnnXteJMJuOsuGbn1p5Ojfv9dRs5Wr3bt9Q4KkhUOHn8cuOce/xk1I7LChIns27FDpu9yc4Fx4+Q05a1bZV/58ra9YayXiunUCbjuOqml6tMHeOopY2On0sucMAUEZKCgwDcSpi1bpM7o8mXXjq9VC/jhB/lS4i8uX5YvQklJtl+e9JQrJ0tIDR4M3HknwJU+yM8xYSLXHTkip/SeOwcMHSoX85p511wjazqZlSunrdkYOBD46y/5Fj11qoxIERVFUlISkpKSkJ+fj7179yIjIwMREfrzOUYlTampQPPm+muXFWYyAY89JgmHTt9en7VnjxRhr1jh/HU1mYAmTSRJuu8+WZOTqJRgwkTF9/ffwIQJMtJ0333Aww9bRp0qVpRicz0tWkhyVb68TON16GBUxFQamEeYvJkwZWVJsbIrIy0VKsgCtv/3f56Nyd1+/FESJeuFxO2JipK//2ee8b/GkUQuYsJE7pOcDMyYITVMmZmWVb1NJvv/gYWGSnKVlyetDR5+WLb70zdwMpZ1whQeHq57VtWIEcC773rm8a3PMHUkOlqm6sxLH/mDggIZAR43DsjIcHxsQIBMKY4bB3TpYkx8RF7EhIk8Iz8f+OADqWG64Qb55nnliuPbmEyWy6RJQNeucjZetWqGhEx+wpWECXD/KNOMGVK07EyNGsDevTKy5C9yc2Vdy0mTnNdhhYVJTeJLL3HKjcoUJkxkjL/+ksWBa9SQ5V+2b3ftdkFB0sOmUSNZBoHNM6lwwjRqlLaHmJk7E6bgYOcL4FauLOtA+tMZYHl5Usj99tvO67Bq1ZJjH33UmNiIfAwTJjJeRgbw9ddSGP7TT7IApytCQ+VbbfXqUjPFs27KpsIJE2C/2WFJk6akJOlT5khQkBRG+1vtzgcfAKNGOR/5vfFGaXvQtq0xcRH5KCZM5F35+ZI05eXJMgrffOPa7WrUAHr0kOm+J59kd+CyxKiEKTJSzgh15OOPpeu9P1m8WDprnz/v+LiOHYHPPvO/RJDIQ5gwke9QCli/XgpP335bztJxVDBuFh8vU3avvCIJFJVuegnTkSNAXJztscVJmHJyZArOkQYNpKWAPzlwALjjDuny78hddwHTp8uXEiK6igkT+SalpG1Bbi7Qvz+wYYNMfeTl2b9NuXJAbCzwwAOSPIWFGRcvGUcvYQL0R5lq1JBkylXz5gH33uv4mF27ZIFYf5GTA/TsCSxd6vi4jh1lqjw21pCwiPwNEybyD5cvy6ncAwfKci2uLD1Rv75liRe90QfyL3qNK50lTIDro0x9+gD/+5/9/bGxlqWA/MUXX8h6bY6+aDRrJtN0/BshcogJE/mnSZOAhQtl+RZn/WLKlQMGDJA1rPyteSDZsjfC9MQTkiAU5spHmHkZH3umTZNk3V+cOCFF2o4WxI2Olmlvf1qehciLmDCRfzt6VJKnDRuAdeucH3/99bJK+pgxzutUyDfZS5gA/VEmZx9hzoq7c3NlOthfvPYaMH68/f1BQcCUKcCQIcbFRFQKMGGi0uPHH2UV9WnTnI86hYdLAezEicC11/IsO39S1IQJsJ80Va7s+GwxbyzkW1xnz0qvsuPH7R/TrZucjcovC0RF5vR/iQAjoiByh+7dJQFKTZXpunbt7B+bmSkFvtddJ2fYzZkDXLhgXKzkfTVr2k+WwsL8K1n67DOgShX7yVLlylKs/tNPTJaIPIUJE/mdmBipVVq3Djh5UprzRUXZP37vXlnDrm5dmapz1p+GxLlz59CvXz9EREQgIiIC/fr1w3knL96AAQNgMpk0l3aOMtsiql/fteN69pTpXD0xMa4tqusLCgqA1q1lqRJ7RoyQ0Sd/OrOPyB9xSo5Kje++k47Fa9bIfzT2hIZKJ/HERPaicaRLly44cuQIpk2bBgB46qmnUKdOHfz44492bzNgwACcOHECM2bMuLqtfPnyiCxC23ZHU3KA7bTcW29JImy2Zg1wyy369x0dLQXT/uCff4CGDe0vaVK5MrBtm4ykEVGJsYaJyp6NG6Ux36xZzhcave02+Q+3fXtjYvMXu3fvRnx8PNavX4+2/66ZsX79eiQkJGDPnj1o2LCh7u0GDBiA8+fPY8GCBcV+bGcJEwAEBMiU2sSJwIsvavfZq3MKCXH+fvAVH3/suGj7ySdlmo6I3IY1TFT2tG4NfPKJnEY+ZYrjNelWrQJuuklO0V692rgYfd26desQERFxNVkCgHbt2iEiIgJr1651eNvk5GRER0ejQYMGGDhwIE6ePOnw+OzsbGRmZmouzhQUSMLkarIE+E+ydM899pOlcuWATZuYLBF5AxMmKrViY4FnnwXOnAGWLweaN7d/7IYNMtrUuLE0+fOngmBPSE9PR3R0tM326OhopKen271dly5d8PXXX+O3337Du+++i40bN+L2229Htr15JQCJiYlX66QiIiIQV8wOi1Wr2t/nD79PpWR6bfFi/f2NGslCui1bGhsXEQkmTFQm3HmndBJPTpZ2A/bs2SPf8OPjJckqbSZMmGBTlF34smnTJgCASWe4Rimlu93sgQceQLdu3dC0aVN0794dS5cuxd69e7HYXhYAYMyYMcjIyLh6SUtLK9ZzO31af7s/JEuZmTJlaK9QfcQIYPdumYokIu/wo3ZtRCV3221y2bVLCoV/+km/QHzPHlmkND5eCslvu834WD1h2LBh6Nu3r8Nj6tSpg23btuGETnX0qVOnEBMT4/LjxcbGonbt2ti3b5/dY4KDgxFcwnPh7eVwhafsfNGePfI+s5fY/forcPvtxsZERLaYMFGZFB8vvZzOngWGDgV++AHIz7c9btcuoEMHWV5i+nRpHOjPoqKiEOWoB8O/EhISkJGRgQ0bNqDNv2tr/Pnnn8jIyED7IlTInzlzBmlpaYj14Iqvjz1mf9/EiR57WLf45RcZ/dQTHAykpwOVKhkaEhHZwQFeKtMiI4FvvgH27QMefND+cRs2yCKmnTvbnzYpTRo3boy7774bAwcOxPr167F+/XoMHDgQ99xzj+YMuUaNGmH+/PkAgKysLIwcORLr1q3DoUOHkJycjO7duyMqKgq9evXyWKwzZ+pv9/WpuG++sZ8s1agh9UpMloh8BxMmIkhTyzlzZLX6Xr3s14qsWAHUqiVnMeXkGBuj0b7++ms0a9YMnTt3RufOnXH99dfjq6++0hyTmpqKjH/XqQkMDMT27dvRo0cPNGjQAP3790eDBg2wbt06hIWFeSTGunX1t/t6YfRHHwEPPaS/r1074MgRY+MhIufYh4lIx759wODBUj9iT4UKsjAwFzl1L1f6MJkVdW05XzBlCjB8uP6+vn1l5ImIDMc+TETFUb++1Jds3Gi/bunSJal/ql9fap3IWC+8oL/9lVeMjaMoPv7YfrL04otMloh8GUeYiFywYAEwcKD9U9cBmcqbPVtGnqj4XB1h8rfRpTlzZE1DPW+8AfznP8bGQ0QaHGEicoeePWWl+DfftL8a/Pz5sqL8F18YGlqpkZSUhPj4eLRu3brY9xER4caA3OjXX+0nS//9L5MlIn/AESaiIjp/HhgwQNoS2HPDDTKlV6WKQUGVIq6MMPnT6NK+fUCDBvr73nwTGDvW2HiISBdHmIjcrVIlmaLbuROoU0f/mJQUIDra9/sAkWddumQ/WXrxRSZLRP6ECRNRMcXHAwcPyihBYKDt/oICYNQooF49oJirfVARTJjg7Qhs2eum8OCDTKaJ/A2n5Ijc4Nw5WUpl40b9/QEB0oJgxAhj4/JHxZ2S87XpuMhIeV8U1r498McfxsdDRA5xSo7ICJUrSzfw2bP1i8ILCuQ0+BtvBC5eND4+MtZtt+knS9WrM1ki8ldMmIjc6OGHgZMngbZt9fdv2SIjD7/9ZmxcpZ0vJaFTpgCrV9tuL1eubCyrQ1RaMWEicrPwcGD9emDWLP3appwcoFMn4KmnjI+ttFBKe/GV3ldHj9pvTOlLSR0RFR1rmIg86MIFoFUrYO9e/f21awM7dgAVKxobly8rytIoviYgQL+Wavt2+x3jicgnsIaJyJvCwoDUVDlbTs8//0j905o1xsbli9zRuNKb6tTRT5bGjWOyRFQacISJyCAbNgAdOgCXL+vvHz/eN0+NN5o/jjC9+y4wcqTt9oYNgT17jI+HiIrM6QgTEyYiA+XkAC1byjScno4dWRDubwnTxYv6U6omk5wdSUR+gVNyRL6kfHmpZ3nhBf39K1cC1aoB2dnGxkXFZ2/5m+PHjY2DiDyLCRORF7zzDrBsmX4DxhMnpPaJ3cF938CB+snt6NFATIzx8RCR53BKjsiLTp4E6tcHMjP19ycnSxPEssRfpuQuXJAWEoWFhwMZGcbHQ0Qlwik5Il8WHS0doRs10t/foQPw+eeGhkQuiorS367X4ZuI/B8TJiIvCwgAdu8G7r1Xf//AgcArrxgbEzn26qtSwF/Y++/L75OISh9OyRH5kFdftd9a4Mkngc8+MzQcQyUlJSEpKQn5+fnYu3evz07JKaWfFIWF2Z9aJSKfx7YCRP7m66+BRx7R33f//cB33xkbj9F8vYapVi39gvzsbDkLkoj8EmuYiPzNww9Lsbee778H7rnH0HDIyp49+snSwIFMlohKO44wEfmoHTuAZs3093XpAixZYmw8RvHlEabAQP1mlK5/jBKRj+IIE5G/atoUOHhQv15m6VKgRw/jYyrLpkzRT5a2bDE+FiIyHkeYiHzc0aNSN6P3n/UDDwDffmt8TJ7kqyNMek1GIyOBM2eMj4WI3I4jTET+rkYN4MgR/f+wv/sOeO4542Mqa/r21d+enm5sHETkPUyYiPxAbKz9pOmDD4DERONjKkv0zkxs3x4oV874WIjIO5gwEfmJ6tWBffv0940dC/zvf8bGU1a0a6e//Y8/jI2DiLyLCRORH6lXD9i4UX9fnz7A1q3GxuNOSUlJiI+PR+vWrb0dylVKAX/+abv9ySeNj4WIvItF30R+aMkSoFs3/X0XLgAVKxobjzv5UtF38+bAtm2229lGgKjUYdE3UWnUtSswebL+vipVjI2ltFJKP1kaO9b4WIjI+5gwEfmp558H+ve33Z6TA9SsaXw8pU3btvrb33zT2DiIyDcwYSLyYzNnAtdfb7v96FGgd2/DwylV9GrFRo82Pg4i8g1MmIj83NatQFiY7fb580tfU0uj9Oqlv53tG4jKLiZMRKXA2bP6PZoefLD4najffPNNtG/fHhUqVEClSpVcuo1SChMmTED16tURGhqKDh06YOfOncULwIsWLLDd9uijhodBRD6ECRNRKRAUJIv16omNLd595uTkoE+fPhg8eLDLt5k4cSImT56MqVOnYuPGjahWrRruvPNOXLhwoXhBeMFbb+lv//JLY+MgIt/CtgJEpUhiov5ZXDfeCGzeXLz7nDlzJoYPH47z5887PE4pherVq2P48OEYNWoUACA7OxsxMTF4++238fTTT+veLjs7G9nZ2VevZ2ZmIi4uzmttBfRG6lq2BDZtMjwUIjIO2woQlSVjxgDNmtlu/+svYNYszz72wYMHkZ6ejs6dO1/dFhwcjNtuuw1r1661e7vExERERERcvcTFxXk2UAe2b9ffbq9ZKBGVHUyYiEqZbduAwEDb7f37A1YDOW6X/u9KtDExMZrtMTExV/fpGTNmDDIyMq5e0tLSPBekE61a2W4LC9MfdSKisoUJE1EpdPSo/vbw8AkwmUwOL5tKOPdkKpRdKKVstlkLDg5GeHi45uINBQXSw6owe+v3EVHZEuTtAIjI/WJigPHjgVdf1W7PyRmGhx/ui5dftn/bOnXqFOsxq1WrBkBGmmKtKs1PnjxpM+rki+w1qvSD0InIAEyYiEqpCROAjz4CTp2y3hqFr7+OwocfApUru/fx6tati2rVqmHFihVo0aIFADnTbtWqVXj77bfd+2AeoDew9sEHxsdBRL6JU3JEpZi90qHq1Z3f9vDhw0hJScHhw4eRn5+PlJQUpKSkICsr6+oxjRo1wvz58wHIVNzw4cPx1ltvYf78+dixYwcGDBiAChUq4KGHHnLH0/GYn37S3/7MM8bGQUS+iyNMRKVYQAAwezbwyCPa7VeuAC++CEyaZP+248aNw5dWzYfMo0YrV65Ehw4dAACpqanIyMi4esxLL72Ey5cvY8iQITh37hzatm2L5cuXI0yvFbkP6dnTdtu11xoeBhH5MPZhIioDatQAjh2z3Z6dDZQvb3w8jmRmZiIiIsLQPkx6NelXrgDBwYY8PBF5H/swERFw5Ij+9uJ2AS9N7rlHfzuTJSKyxoSJqAwwmYD33rPdfvYs8OOPxsejJykpCfHx8WjdurWhj7t4se22wmcXEhFxSo6oDLnmGuDSJdvtrn8MeJ6RU3LnzgGRkbbbfen1ICJDcEqOiCzsTc3de6+xcfiKpk1tt3Eqjoj0MGEiKkMqVwb+PcFNY948IC/P8HC8Tq8QfvVq4+MgIt/HhImojFm5Un977drGxuFtu3bpb2/Txtg4iMg/MGEiKoP0OlgfO2Z/yq40atfOdlvVqsbHQUT+gUXfRGVUYKAsOGstOFj6D3mTUUXfer2Xjh1jqwWiMopF30Skb+NG223Z2cAffxgfi9G2btXfzmSJiOzhCBNRGRYWBlgtDQdARl4KjzwZyYgRpogIIDNTu61aNeD4cY88HBH5Po4wEZF9//xju00p4N/1dA1lZOPKwskSYH/UiYgI4AgTUZkXEwOcPGm73VvNGz09wnT0KFCzpu12NqskKtM4wkREjumNMgHeGWUygt4AlkFr/BKRH2PCRFTGhYToFzv37m18LEbQq1Nas8b4OIjIvzBhIiL8/bf+9uRkQ8PwuPx8/e3NmhkbBxH5HyZMRISQEKBKFdvtnToZH4sndeliuy2An4JE5AJ+VBARAODgQdttBQXAoUOGh+IxK1bYbnvvPePjICL/w7PkiOiq0FDbTt+hocClS8bF4Mmz5PS6e/PsOCICz5IjoqLYudN22+XLQG6u8bG4W2k964+IjMGEiYiuuvZa/e1Nmnj+sT3duPKhh2y3sdibiFzFKTki0pgzB3j4YdvtRk1deWpKTm86LisLuOYatz0EEfkvTskRUdHojcQAwPPPGxuHO9lL9pgsEZGrmDARkY1evWy3vf++4WG4jV4SqDfiRERkDxMmIrIxb57+9gMHjI3DXb791nbb6NHGx0FE/os1TESkKzwcuHBBu61CBeDiRc8+ridqmNhOgIicYA0TERXPjh2224zsx+Quly97OwIiKg2YMBGRrlq19Lc//bSxcZTU7bfbbgsONj4OIvJvTJiIyK7+/W23TZvmmcfyVB+m9ettt338sVsfgojKANYwEZFDevU/Fy4AFSt65vHcXcPE+iUicgFrmIioZIKCbLcZ0fnbHbKzvR0BEZUWTJiIyKFly2y3HT5sfBzFcffdttvKlTM+DiLyf0yYiMihTp30t2/fbmwcxZGcbLtt0iTDwyCiUoA1TETkVNWqwOnT2m2hoZ5pM+DOGibWLxGRi1jDRETF8+abb6J9+/aoUKECcnIq2ezX6280YMAAmEwmzaVdu3aeD5aIyMOYMBGRrpycHPTp0weDBw+2u+7aihW22+6++24cP3786mXJkiWeDdQOLn1CRO6kc/4LERHw6quvAgBmzpwJAKhbFzh4UHtMt25ATo52W3BwMKpVq+by42RnZyPb6nS2zMzMYsVb2Lvv2m675x633DURlUEcYSIil2zbZrstN9d2W3JyMqKjo9GgQQMMHDgQJ0+edHi/iYmJiIiIuHqJi4tzS7x5ebbbFixwy10TURnEhImIXGKvUeWaNZafu3Tpgq+//hq//fYb3n33XWzcuBG33367ZgSpsDFjxiAjI+PqJS0tzc2RWwQGeuyuiaiUY8JEVIZMmDDBpii78GXTpk12b1+jhu22zp0tPz/wwAPo1q0bmjZtiu7du2Pp0qXYu3cvFi9ebPc+g4ODER4errmUlL/0iSIi/8EaJqIyZNiwYejbt6/DY+rUqWN3X0qKtBiwpne2nFlsbCxq166Nffv2uR6kG1gncWZhYYaGQESlDBMmojIkKioKUVFRJbi9/vaTJ4HoaNvtZ86cQVpaGmJjY4v9mMWRmmq7bfp0Q0MgolKGU3JEpOvw4cNISUnB4cOHkZ+fj5SUFKSkpCAkJMvqqEYA5qNFCyArKwsjR47EunXrcOjQISQnJ6N79+6IiopCr169vPU0rurTx9sREJE/4wgTEekaN24cvvzyy6vXW7RoAQD44IOVePbZDv9uTQWQgWPHgMDAQGzfvh2zZs3C+fPnERsbi44dO+K7775DGOfDiMjPcWkUIioyTy454o6lUR59FPjqK8v1p54CPv3UPfERUanEpVGIyBgDBng7AotZsyRBuuceYM4cJktEVHIcYSKiInviCeCLL2y3l2SUKSkpCUlJScjPz8fevXvdsvguEZGLnI4wMWEiomLx1LScO6bkiIiKiFNyRGQcNowkotKKCRMRFYveUilt2hgfBxGREZgwEVGxLFtmu+3ECePjICIyAhMmIiqWm27ydgRERMZhwkREbrV0qbcjICJyPyZMRFRscXG223r3Nj4OIiJPY8JERMW2bp3ttitXjI+DiMjTmDARUbHVqOG++0pKSkJ8fDxat27tvjslInITNq4kohLRa2D5889A587Fuz82riQiL2DjSiLyrNhY2229ehkfBxGRJzFhIqISWbXKdtulS8bHQUTkSUyYiKhE6tf3dgRERJ7HhImIPOLUKW9HQETkPkyYiKjEype33da+vfFxEBF5ChMmIiqxpCTbbfv3Gx8HEZGnsK0AEbmFXnsB1z9eLNhWgIi8gG0FiMg/sHElEfkyjjARkVvojTBNnQoMHVq0++EIExF5AUeYiMgY9erZbnv5ZePjICLyBCZMROQWycm224KCDA+DiMgjmDARkVvUrGk7LTd5sndiISJyNyZMROQ2Z88CDRoAVaoAkyYB/fp5OyIiIvfggDkRuU2lSkBqqrejICJyP44wERERETnBhImIiIjICSZMROQT2LiSiHwZG1cSkU9h40oi8gI2riQiIiIqKSZMRERERE4wYSIiIiJyggkTEdk4dOgQnnjiCdStWxehoaGoV68exo8fj5ycHIe3U0phwoQJqF69OkJDQ9GhQwfs3LnToKiJiDyHCRMR2dizZw8KCgrw6aefYufOnXjvvffwySefYOzYsQ5vN3HiREyePBlTp07Fxo0bUa1aNdx55524cOGCQZETEXkGz5IjIpdMmjQJH3/8Mf7++2/d/UopVK9eHcOHD8eoUaMAANnZ2YiJicHbb7+Np59+2qXH4VlyROQFPEuOiNwjIyMDkZGRdvcfPHgQ6enp6Ny589VtwcHBuO2227B27Vq7t8vOzkZmZqbmQkTka5gwEZFTBw4cwIcffohBgwbZPSY9PR0AEBMTo9keExNzdZ+exMREREREXL3ExcW5J2giIjdiwkRUhkyYMAEmk8nhZdOmTZrbHDt2DHfffTf69OmDJ5980uljmEzakW2llM02a2PGjEFGRsbVS1paWvGeHBGRBwV5OwAiMs6wYcPQt29fh8fUqVPn6s/Hjh1Dx44dkZCQgGnTpjm8XbVq1QDISFNsbOzV7SdPnrQZdbIWHByM4OBgF6InIvIeJkxEZUhUVBSioqJcOvbo0aPo2LEjWrZsiRkzZiAgwPGAdN26dVGtWjWsWLECLVq0AADk5ORg1apVePvtt0scOxGRN3FKjohsHDt2DB06dEBcXBzeeecdnDp1Cunp6Ta1SI0aNcL8+fMByFTc8OHD8dZbb2H+/PnYsWMHBgwYgAoVKuChhx7yxtMgInIbjjARkY3ly5dj//792L9/P2rWrKnZZ92KJDU1FRkZGVevv/TSS7h8+TKGDBmCc+fOoW3btli+fDnCwsIMi52IyBPYh4mIfAr7MBGRFzjtw1SUhImIyONMJlM4gAwAEUopNmUiIp/AhImIfIpJehCEAbig+AFFRD6CCRMRERGREzxLjoiIiMgJJkxERERETjBhIiIiInKCCRMRERGRE0yYiIiIiJxgwkRERETkBBMmIiIiIieYMBERERE58f8TMnoAxvACzwAAAABJRU5ErkJggg==\n", "text/plain": [ "Graphics object consisting of 3995 graphics primitives" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plot the computed points on the contour\n", "show(pt, xmin=-3, xmax=3, ymin=-2, ymax=2)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "x^10 + 10*x^9*y + 45*x^8*y^2 + 127*x^7*y^3 + 252*x^6*y^4 + 357*x^5*y^5 + 356*x^4*y^6 + 237*x^3*y^7 + 93*x^2*y^8 + 17*x*y^9 + y^10 + x^9 + 9*x^8*y + 36*x^7*y^2 + 90*x^6*y^3 + 156*x^5*y^4 + 186*x^4*y^5 + 147*x^3*y^6 + 69*x^2*y^7 + 15*x*y^8 + y^9 + x^8 + 8*x^7*y + 28*x^6*y^2 + 61*x^5*y^3 + 90*x^4*y^4 + 86*x^3*y^5 + 49*x^2*y^6 + 13*x*y^7 + y^8 + x^7 + 7*x^6*y + 21*x^5*y^2 + 39*x^4*y^3 + 47*x^3*y^4 + 33*x^2*y^5 + 11*x*y^6 + y^7 + x^6 + 6*x^5*y + 15*x^4*y^2 + 23*x^3*y^3 + 21*x^2*y^4 + 9*x*y^5 + y^6 + x^5 + 5*x^4*y + 10*x^3*y^2 + 12*x^2*y^3 + 7*x*y^4 + y^5 + x^4 + 4*x^3*y + 6*x^2*y^2 + 5*x*y^3 + y^4 + x^3 + 3*x^2*y + 3*x*y^2 + y^3 + x^2 + 2*x*y + y^2 + x + y + 1" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The power series expansion of 1/R begins\n", "(1/R).taylor((x,0),(y,0),10)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1, 2, 6, 23, 90, 357, 1443, 5910, 24426, 101669]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The power series diagonal sequence starts\n", "ser = QQ[x,y]((1/R).taylor((x,0),(y,0),21))\n", "[ser[k,k] for k in range(10)]" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-x^4*y^7 - 4*x^4*y^4 + x^3*y^5 + 4*x^3*y^4 - 6*x^4*y + 3*x^3*y^2 - x^2*y^3 + 12*x^3*y - 3*x^2*y^2 - 6*x^2*y - 2*x^2 - 4*x^4/y^2 + 3*x^3/y + x*y + 2*x + 12*x^3/y^2 - 6*x^2/y - 12*x^2/y^2 + 3*x/y - x^4/y^5 + x^3/y^4 - x^2/y^3 + 5*x/y^2 - 1/y + 4*x^3/y^5 - 3*x^2/y^4 + 2*x/y^3 - 1/y^2 - 6*x^2/y^5 + 3*x/y^4 - 1/y^3 + 4*x/y^5 - 1/y^4 - 1/y^5" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Similarly, the Laurent expansion of 1/R when 1 + |x| + |xy^3| < |y| has terms\n", "(-1/y * add([((1-x-x*y^3)/y)^k for k in range(5)])).expand()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0, 1, -3, 10, -39, 156, -630, 2577, -10647, 44308]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# This Laurent expansion has diagonal sequence beginning\n", "laurent_ser = (-1/y * add([((1-x-x*y^3)/y)^k for k in range(21)])).expand()\n", "[laurent_ser.coefficient(x,k).coefficient(y,k) for k in range(10)]" ] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.2", "language": "sage", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }