{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Code to Compute Asymptotic Contributions of Smooth Contributing Points\n", "The code below computes asymptotic expansions of strictly or finitely minimal smooth contributing points, using Theorem 5.2 and Corollary 5.2 from the textbook. *Note:* The [multivariate_generating_functions (mgf) package](https://doc.sagemath.org/html/en/reference/asymptotic/sage/rings/asymptotic/asymptotics_multivariate_generating_functions.html) has the capability to perform similar computations, however for completeness (and because the mgf package currently needs work to fix various errors) we implement Theorem 5.2 directly using the DifferentialWeylAlgebra package. \n", "\n", "*Requirements: None*" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Set a parameter to help simplify some algebraic numbers\n", "maxima_calculus('algebraic: true;')\n", "\n", "# Procedure to get Hessian appearing in asymptotics\n", "# Input: H, member of the symbolic ring\n", "# r, direction vector (which can contain symbolic entries)\n", "# vars, vector of variables\n", "# CP, a dictionary mapping elements of vars\n", "# Output: The Hessian H defined in Lemma 5.5 of the textbook at the point w defined by CP\n", "def getHes(H,r,vars,CP):\n", " dd = len(vars)\n", " V = zero_vector(SR,dd)\n", " U = matrix(SR,dd)\n", " M = matrix(SR,dd-1)\n", "\n", " for j in range(dd):\n", " V[j] = r[j]/r[-1]\n", " for i in range(dd):\n", " U[i,j] = vars[i]*vars[j]*diff(H,vars[i],vars[j])/vars[-1]/diff(H,vars[-1])\n", " for i in range(dd-1):\n", " for j in range(dd-1):\n", " M[i,j] = V[i]*V[j] + U[i,j] - V[j]*U[i,-1] - V[i]*U[j,-1] + V[i]*V[j]*U[-1,-1]\n", " if i == j: M[i,j] = M[i,j] + V[i]\n", " return M(CP)\n", "\n", "# Procedure to apply differential operator to f and set all variables to zero\n", "# Input: dop, element of a DifferentialWeylAlgebra over a polynomial ring\n", "# f, an element of the base polynomial ring of dop\n", "# Output: dop(f) evaluated when all variables are zero\n", "def eval_op(dop, f):\n", " if len(f.parent().gens()) == 1:\n", " return add([prod([factorial(k) for k in E[0][1]])*E[1]*f[E[0][1][0]] for E in dop])\n", " else:\n", " return add([prod([factorial(k) for k in E[0][1]])*E[1]*f[(v for v in E[0][1])] for E in dop])\n", "\n", "# Procedure to get critical points of rational function with denominator H, in direction r\n", "# Input: H, member of the symbolic ring\n", "# r, direction vector (which can contain symbolic entries)\n", "# vars, vector of variables\n", "# Output: Solutions (if found by solve) of the smooth critical point equations of H in the direction r\n", "def critpt(H,r,vars):\n", " d = len(vars)\n", " criteqs = [H] + [r[j]*vars[0]*diff(H,vars[0]) - r[0]*vars[j]*diff(H,vars[j]) for j in range(1,d)]\n", " return solve(criteqs,vars,solution_dict=true)\n", "\n", "# Procedure to compute asymptotic contribution of a strictly minimal contributing point\n", "# Input: G, member of the symbolic ring\n", "# H, member of the symbolic ring\n", "# r, direction vector (which can contain symbolic entries)\n", "# vars, vector of variables\n", "# CP, a dictionary mapping elements of vars to coordinates of a strictly minimal contributing point\n", "# M, positive integer describing the number of terms in the asymptotic expansion to compute\n", "# g, parametrization of variable vars[-1] near CP, in terms of the remaining variables\n", "# Output: ex, pw, se such that ex*pw*(se+O(n^(M-1)) gives an asymptotic expansion of the r-diagonal of \n", "# G/H in the variables vars, to order M.\n", "# NOTE: Unlike the textbook, M here refers to the number of terms in the expansion\n", "# (not the order of the expansion, so M should be at least 1)\n", "def smoothContrib(G,H,r,vars,CP,M,g):\n", " # Preliminary definitions\n", " dd = len(vars)\n", " field = SR\n", " tvars = list(var('t%d'%i) for i in range(dd-1))\n", " dvars = list(var('dt%d'%i) for i in range(dd-1))\n", "\n", " # Define differential Weyl algebra and set variable names\n", " W = DifferentialWeylAlgebra(PolynomialRing(field,tvars))\n", " WR = W.base_ring()\n", " T = PolynomialRing(field,tvars).gens()\n", " D = list(W.differentials())\n", "\n", " # Compute Hessian matrix and differential operator Epsilon\n", " HES = getHes(H,r,vars,CP)\n", " HESinv = HES.inverse()\n", " v = matrix(W,[D[k] for k in range(dd-1)])\n", " Epsilon = -(v * HESinv.change_ring(W) * v.transpose())[0,0]\n", "\n", " # Define quantities for calculating asymptotics\n", " tsubs = [v == v.subs(CP)*exp(I*t) for [v,t] in zip(vars,tvars)]\n", " tsubs += [vars[-1]==g.subs(tsubs)]\n", " P = (-G/g/diff(H,vars[-1])).subs(tsubs)\n", " psi = log(g.subs(tsubs)/g.subs(CP)) + I * add([r[k]*tvars[k] for k in range(dd-1)])/r[-1]\n", " v = matrix(SR,[tvars[k] for k in range(dd-1)])\n", " psiTilde = psi - (v * HES * v.transpose())[0,0]/2\n", "\n", " # Recursive function to convert symbolic expression to polynomial in t variables\n", " def to_poly(p,k):\n", " if k == 0:\n", " return add([a*T[k]^int(b) for [a,b] in p.coefficients(tvars[k])])\n", " return add([to_poly(a,k-1)*T[k]^int(b) for [a,b] in p.coefficients(tvars[k])])\n", "\n", " # Compute Taylor expansions to sufficient orders\n", " N = 2*M\n", " PsiSeries = to_poly(taylor(psiTilde,*((v,0) for v in tvars), N),dd-2)\n", " PSeries = to_poly(taylor(P,*((v,0) for v in tvars), N),dd-2)\n", "\n", " # Precompute products used for asymptotics\n", " EE = [Epsilon^k for k in range(3*M-2)]\n", " PP = [PSeries] + [0 for k in range(2*M-2)]\n", " for k in range(1,2*M-1):\n", " PP[k] = PP[k-1]*PsiSeries\n", " \n", " # Function to compute constants appearing in asymptotic expansion\n", " def Clj(l,j):\n", " return (-1)^j*SR(eval_op(EE[l+j],PP[l]))/(2^(l+j)*factorial(l)*factorial(l+j))\n", " \n", " # Compute different parts of asymptotic expansion\n", " var('n')\n", " ex = (prod([1/v^k for (v,k) in zip(vars,r)]).subs(CP).canonicalize_radical())^n\n", " pw = (r[-1]*n)^((1-dd)/2)\n", " se = sqrt((2*pi)^(1-dd)/HES.det()) * add([add([Clj(l,j) for l in range(2*j+1)])/(r[-1]*n)^j for j in range(M)])\n", " \n", " return ex, pw, se.canonicalize_radical()\n", "\n", "# Procedure to aid in printing an asymptotic expansion\n", "# Procedure to get critical points of rational function with denominator H, in direction r\n", "# Input: ex,pw,se as returned by smoothContrib(G,H,r,vars,CP,M,g)\n", "# Output: None (function pretty prints the asymptotic expression defined by ex,pw,se, and M)\n", "def disp_asm(ex,pw,se,M):\n", " show(ex*pw,LatexExpr(\"\\\\Bigg(\"), se, LatexExpr(\"+ O\\\\Bigg(\"), n^(-M), LatexExpr(\"\\\\Bigg)\\\\Bigg)\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 5.1 (Asymptotics of Central Binomial Coefficients)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The set of critical points [{x: r/(r + s), y: s/(r + s)}] contains only one point\n", "The asymptotic expansion begins\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\left(\\frac{{\\left(r + s\\right)}^{r + s}}{r^{r} s^{s}}\\right)^{n}}{\\sqrt{n s}} \\Bigg( \\frac{\\sqrt{2} \\sqrt{r + s}}{2 \\, \\sqrt{\\pi} \\sqrt{r}} + O\\Bigg( \\frac{1}{n} \\Bigg)\\Bigg)\n", "\\end{math}" ], "text/plain": [ "((r + s)^(r + s)/(r^r*s^s))^n/sqrt(n*s) \\Bigg( 1/2*sqrt(2)*sqrt(r + s)/(sqrt(pi)*sqrt(r)) + O\\Bigg( 1/n \\Bigg)\\Bigg)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Define parameters\n", "var('x,y,r,s')\n", "G = 1; H = 1-x-y\n", "vars = [x,y]\n", "R = [r,s]\n", "M = 1\n", "\n", "# Get critical point as function of r and s\n", "CPs = critpt(H,R,vars)\n", "print(\"The set of critical points {} contains only one point\".format(CPs))\n", "CP = CPs[0]\n", "\n", "# Parametrize y = g(x) = 1-x on V(H)\n", "g = 1-x\n", "\n", "# Get and print asymptotics\n", "ex,pw,se = smoothContrib(G,H,R,vars,CP,M,g)\n", "print(\"The asymptotic expansion begins\")\n", "disp_asm(ex,pw,se,M)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 5.2 (‘Asymptotics’ in an Irrational Direction)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAGFCAYAAADtt7dbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA7T0lEQVR4nO3deXiU1f3+8feHLezBFRFBQWWJqBQEV1BBUGtdEHcRqbsiVmt/tdhasK3E1l2calWKW92+LiBU3HdZyiKgCOICArIpYoIKAcL5/XEmziSGIROSObPcr+uaK8+ZPDNzk0fMh3POc4455xARERGRytUJHUBEREQknalYEhEREUlAxZKIiIhIAiqWRERERBJQsSQiIiKSgIolERERkQRULImIiIgkoGJJREREJAEVSyIiIiIJZE2xZF5zM7PQWURERCR71Ev2BWbWG/h/QHegFTDAOTcuwfmtgNui5+8L3O2cu7rCOUOAsZW8vJFzbn1VchUVFZGfn09RUVFVThcREZHclHSnSnV6lpoAc4Arq3h+HvA1cFP0dVtTjC++fno45zZUI5+IiIhIjUm6Z8k5NwmYBFCVES/n3GLgN9HzL0h8qluZbB4RERGR2pR0sVSLmprZl0BdYDZwg3Pug62dXFJSQklJyU/t4uLiWg8oIiIiuSddJngvAIYAJwFnAxuA981s3629oLCwkPz8/J8ebdq0SU1SERERySnmnKv+i80c25jgXeH8t4DZFSd4V3JeHWAW8I5zblhl51TsWVq1qpgOHdrw9ttFNG3avGp/ABHJeVOnwsSJsMsucPHF0Lhx6EQiUl2dOlXp73DSE7zTaRjuJ865LWY2HX/3XKXy8vLIy8v7qT17tv965JG1HE5EstYjj4ROICLbY+ZM6Nat5t83LYul6FpJXYEPq/qaDh3817ffhqZNayWWiGSZZ5+FUaNi7RYt4PXXg8URke3UqVPtvG911llqCuwT91Q7M+sKfOucW2JmhUBr59zguNd0jR42BXaJtjc65z6Ofn8EMBX4FGgOXIUvloYCv65KrrJut65doblG4UTS1qpV8MwzsOOOcOaZUCfgzEnn4O9/h9JS3+7Vq3b+VSoima06PUsHAW/GtW+Pfn0YP0m7FdC2wmvi72rrDpwDfAnsFX2uBXA/sBtQFD2/t3Puf9XIJyJpas0a6NEDli717ddegzFjwuXp3h3Gj4eHH4bdd4e//CVcFhFJX9s1wTtFqhSwuLj4pxW8m6trSSQtPfMMnH56rF2/PpSUgDYpEpEUSskK3iKSQbZsgRkzYP780Emg4gofbdqoUBKR9KdiSSSLbdkCAwb4oa+CAhg5Mmyegw+GO++EPff0c4OeeSZsHhGRqtAwnEgWe+89P2m5jBn88AM0ahQuk4hIYBqGE0kHmzfH7rAKKW4pMgDq1YO6dcNkERHJVBlfLEUiEQoKCujRo0foKCKAH2Zq3BiaNAl7pxf44bdh0TXw69WDe++FBg3CZhIRyTQahhOpQV9+Ce3a+fV7wBcoq1b5NYVCWrPG9zJpwVYREQ3DiQS1bl2sUAI/HPfDD+HylNlpJxVKIiLVpWJJMt577/m7rLp1g0mTwmbZbz846aRYe9Cgn98uLyIimUXDcJLRfvwRWreG777z7UaNYNEiaNkyXKYtW+DNN/0QXO/eWkdIRCTNJP1/5bTcSFekqr75JlYoAaxfD199FbZYqlMH+vYN9/kiIlKzNAwnSVu5Eq64AgYPhlmzwmbZYw849NBYu6DAP0RERGqKhuEkaQceCHPn+uP8fL+NRqtW4fJ8/z08+CBs2gQXXhj+zjMREUlrSQ/DqViSpKxbBxV/vC+/DP37h8kjIiKSJC0dkI1KS/3ihn//OyxZEjZLs2bQpUus3bx5+baIiEi2Uc9SBjjvPHjsMX/csiXMmRN2AvPy5X5D1u+/h2uu8atEi4iIZAgNw2VjsdSggZ+PU+bJJ+HMM8PlERERyWC5NwxXW3vDzZoFEyZAUVGNvm217L137NgM2rcPl0VERCTXqGepEnffDb/5jT/eZx+YOtVvFxHK/Plw2WV+TaHf/AYuuSRcFhERkQynYbiaKJb22MMvbFjm/vvh4ou36y1FREQkPeTeMFxtyM8v327RIkgMERERSQMqlioxZoxfZLFOHRgyBAYODJ1IREREQtHecJU45BB/e/zmzX4zVBEREcld6llKQIWSiIiIqFgSERERSUDFkoiIiEgCKpZEREREElCxJCIiIinx/fdwzjl+wecLL4SSktCJqibjpzBHIhEikQilpaWho4iIiEgCf/oTPPGEP/78c9hzT/jzn8NmqoqM71kaOnQoH3/8MdOnTw8dRURERBJYtKh8e/HiIDGSlvHFkoiIiFTOOSgshGOPhRtu8OsHhnTWWX5DePALP59+etg8VZXxw3AiIiJSuXvugeuv98evvAJ168LIkeHynH02tGwJM2bAEUfAYYeFy5IMFUsiIiJZaubM8u1Zs8LkiNenj39kEg3DiYiI1KAXXoDLLoNIxA+DhXTMMeXbffuGyZHp1LMkIiJSQ158EU4+Odb++uuww16DBvm5QW+/DT16wEUXhcuSydSzJCIiUkNeey1xO4RzzoF//UuF0vZQz5KIiGS0jz6CiRP9QoennRY2S7du5du/+EWYHFKzVCyJiEjG+vBDOOQQ+PFH377xxrCLHA4aBN9844fjunSBUaPCZZGaYy707LNtq1LA4uJi8vPzKSoqonnz5rWdSURE0sBNN/lVocsUFMC8eeHySEawZF+gOUsiIpKUr7/2W1a8807oJLD33uXb7duHySHZLeOH4bQ3nIhI6qxa5e+qWrrUt0eNguHDw+U56yzfk/R//+cLp/vvD5dFspeG4UREpMr+9S+/hlCZVq1g+fJweUSqQcNwIiLZZuNGeO89+Pjj0Elgl13Kt3fdNUwOkVRSsSQiksY2bICjjoJevWC//eD228PmGTAAhg2DJk2gQwd46KGweURSQcNwIiJp7Pnn4dRTY+0mTeD778PlEckCGoYTEakJixbBypWhU0CjRuXbjRuHySGSy1QsiYjEcQ7OO8/fgt66NYweHTbPscfCkCH+uHFjePDBoHFEclLSxZKZ9TazCWa23MycmZ2yjfNbmdnjZvaJmW0xszu3ct5AM/vYzEqiXwckm01EZHtNngyPPeaPt2yB3/7WT7AOxQzGjoXvvoO1a+Gkk8JlEclV1elZagLMAa6s4vl5wNfATdHX/YyZHQo8BTwKHBj9+rSZHVyNfCKSgdavh5KS0Cl8z1LFdjpM7czPhwYNQqcQyU1JF0vOuUnOuT85556r4vmLnXO/cc49AhRt5bSrgVedc4XOuQXOuULg9ejzIpLlbrrJT1xu2hTuvTdslsMP9wsdgu/VueUWyMsLm0lEwtquu+HMzAEDnHPjqnj+W8Bs59zVFZ5fAtzhnLsj7rlrgKudc20re6+SkhJK4v4ZWlxcTJs2bXQ3nEiG+ewz2HffWLtuXVi9GnbcMVwmgIUL/RyhPfYIm0NEalzG3g23G7CqwnOros9XqrCwkPz8/J8ebdq0qdWAIlI7ynaLL1Na6tcWCq1DBxVKIuKlS7EEP19PySp57ifDhw+nqKjop8fSso2KRGSbnnzS99zk58MDD4TNsv/+MHBgrH3JJbD77uHyiIhUlC4b6a7k571Iu+J7lyodhsvLyyNPEwlEkvbdd3D++bE7vC6/HI47DkJ1zprB00/DlClQrx4crNs6RCTNpEvP0hSgX4Xn+gOTA2QRyWrFxeVvhS8t9QVUSHXq+InVKpREJB1VZ52lpmbW1cy6Rp9qF223jX6/0MweqfCasvObArtE2wVxp9wF9Dez68ysk5ldBxwD3Jn8H0kk/bz7LhQUwJ57wpgxYbO0bVt++4xjjvHZRESkcknfDWdmRwFvVvKth51zQ8zsIWAv59xRca+p7EO+dM7tFXfOacDfgPbA58Afo8sTaG84yWilpX5n9m+/9e06dfzu8R07hs00aRJs3gwnnAD164fLIiKSYknfDaeNdEVq2bp1UPE/yTfegKOPDpNHRCTHZezSASI1asEC6NMHDjzQbxURUrNm5Ye9OneGnj3D5RERkeSoZ0myUufOvmACf7fVjBnQrVu4PJs3w1NPwfffw5lnQosW4bKIiOS4pHuW0mXpAJEa45xfFTq+/cUXYYulevXg3HPDfb6IiFSfhuGkRqxc6ffT6tULHnoobBYzOO20WLtlS59LRESkOjQMJzWiTx94M3qPpBm88w4ccUS4PJs3+1v016zxPTp77hkui4iIpJXcG4aLRCJEIhFKS0tDR8lpH34YO3YO5s0LWyzVqweXXhru80VEJHuoZylDrVsH118Pn3/uJwyff37YPEOGwMMP++PGjWHWrLDrCImIiGxF7vUs5arLLoPHH/fHkyZBq1bQv3+4PA884G/TX74czj5bhZKIiGQPFUsZasaM8u2ZM8MWS/XrwzXXhPt8ERGR2qK74apo0ya48UY4/fTYcFNIffrEjuvUgSOPDJdFREQkm6lnqYr+8Ae4/XZ//MwzP1+VOdXuvhv22svPWTrtNDjssHBZREREspmKpSqaMqV8e/LksMVS/fpw3XXhPl9ERCRXaBiuig4/vHw75G3xIiIikjrqWaqiwkK/n9dHH8EJJ8App4ROJCIiIqmgdZZEREQklyS9zpKG4UREREQSyPhiKRKJUFBQQI8ePUJHERERkSykYTgRERHJJRqGExERkczywgswdCjcd5/fjD3d6G44ERERCWbiRDj55Fj766/hhhvC5amMepZEREQkmNdeK99+9dUwORJRsSQiIpJj3n4brr4aIhHYsiVslq5dE7fTgYbhREREcsiUKdC3L5SW+vaiRXDrreHyDBkC33wDL70E++/vF4FON7obTkREJIf85S8wYkSs3aULfPhhuDwBJH03nHqWREREatn06fDkk9CmDVx5JdQL+Nv3wAPLt/ffP0yOTKJiSUREpBbNmwe9e8OGDb794YcwZky4PCefDKNHwzPPwL77hh2CyxQahhMREalF99wDw4bF2rvvDl99FS6PaBhORESEefPgoYdgl13gqqugYcNwWSoOc3XpEiaHVF/GF0uRSIRIJEJp2bR+ERHJaV9+CYcfDkVFvj11Kjz3XLg8Rx7ph90eesjPWbrjjnBZpHo0DCciIlnl8cfh3HNj7by82HwhETQMJyIiIXz2Gdx/PzRtCtdcA82ahctSUAB16sQWW9Swl2wvFUsiIrJdvv7aD3utXu3br7/uV4gOpWtXeOIJvzr1LrvA7beHyyLZQcNwIiKyXV56CY4/vvxzP/wAjRuHySOyDRqGExHJBUuW+LVy6teH3/4Wdt45XJYOHfy8oJIS3957bxVKkl1ULImIZJh166BXL18wAUycCB98AHXrhsnTvj2MGwf/+Ac0bw633BImh0ht0TCciEiGmT4devYs/9zSpbDHHmHyiGSYpIfh6tRGChGRbLNiBVx7rb/Tq6xHJ5S99vI9OGVatfITmUWkdmgYTkRkGzZuhKOOgoULffv55+Hjj8PNy9llF5g0ye8e36ABjBrl5wyJSO3QMJyIyDZ8/jnss0/55+bMgQMOCJNHRLZL7g3DRSIRCgoK6NGjR+goIlKDVq/2Q16XXup7cUJq1Qpatoy1d9wR9twzXB4RSS31LIlIWura1ffeAOy0E8yfH3Zezocfwp//DKWlMGIEdO8eLouIbJeke5ZULIlI2lm71vfexHvjDTj66DB5RCSr5N4wnIjUjDVr4Kqr4Lzz4P33w2Zp0aL8HKFmzaBTp2BxRCTHqWdJRAA44ohYkdS4sR92at8+XJ7Fi+GGG/y2Gb//PRxySLgsIpJVNAynYkkkec5BvXqxXdoBnnkGBg4Ml0lEpJZoGE4kU6xdC1dcAaecAi+8EDaLWfmem0aN/ARrERFRz5JIMMcf73drB7+n17RpYe+w+uYbf5fX2rVw+eV+7zERkSxU+z1LZtbbzCaY2XIzc2Z2ShVec6SZzTSzDWb2hZldVuH7Q6LvVfHRMNl8Ipnif/+LHZeWwqxZ4bKA37U+EoHHH1ehJCISrzrDcE2AOcCVVTnZzNoBLwLvAr8ARgF3m1nF2RDFQKv4h3NuQzXyiVTqu+/goougTx944IHQafz2GWUaNIBDDw0WRUREEkh6bzjn3CRgEoBZlXqyLgOWOOeujrbnm9lBwO+AZ8u/tVuZbB6Rqrr4Yj9pGeDNN6FNGzjuuHB5HnkEOneG5cvh/POhS5dwWUREZOtSsZHuocArFZ57GbjQzOo75zZFn2tqZl8CdYHZwA3OuQ+29qYlJSWUlJT81C4uLq7R0JJ9Zs8u354zJ2yx1KQJ/O1v4T5fRESqJhV3w+0GrKrw3Cp8obZztL0AGAKcBJwNbADeN7N9t/amhYWF5Ofn//Ro06ZNjQeX7VNUBIMH+0nLI0eGTgPHHhs7rlfPD8eJiIhsy3bdDWdmDhjgnBuX4JyFwFjnXGHcc4cD7+HnJf1s6M3M6gCzgHecc8Mqe9/KepbatGmju+HSyJAh8PDDsfbYsf65UDZvhnvugS++gNNOg969w2UREZFgkr4bLhXDcCvxvUvxdgU2A2sqe4FzbouZTQe22rOUl5dHXl5ejYWUmrdgQfn2/PlhcpSpVw+uvjpsBhERyTypGIabAvSr8Fx/YEbcfKVyzM8c7wqsqN1o2aW4GM44A/be209m3rgxbJ6TT44d160Lv/pVuCwiIiLVlXTPkpk1BeK2uKSdmXUFvnXOLTGzQqC1c25w9Pv3AVea2e3AA/gJ3xfi5yaVvecIYCrwKdAcuApfLA0Ffp1sxlw1fDj83//54y++8Pt6DR8eNs+ee8K8eX4BxiOOCJdFRESkuqozDHcQ8GZc+/bo14fxk7RbAW3LvumcW2RmvwTuwBc/y4GrnHPxywa0AO7HD9cVAR8AvZ1zccv2ybYsWpS4HcI554ROICIisn203cl2WLcOzj0X3n3X76v1xBPQokVKPrpSjz0G553nj+vW9VtpHHNMuDwiIiJpKOkJ3iqWtsPvfw+33BJrDxsGd9+dko/eqtdf99tm9O4NBx8cNouIiEgaSsu74bLWygqLHqxIg+noffv6h4iIiNSMVNwNl7WGDPF7eoG/Lf2CC4LGERERkVqgnqXt0KcPzJjhd4/v3h26dg2dSERERGpaxhdLkUiESCRCaWlpkM/ff3//EBERkeykCd4iIiKSS5Ke4K05SyIiIiIJqFgSERGRtLNxI3zwASxfHjqJiiURERFJM99/D4cdBt26Qbt28MwzYfOoWBIREZG08vjjMHOmP964Ef7wh7B5VCyJiIgIP/wAkyb5JXFCq1cvcTvVVCyJiIjkuLJhr1/+Enr0gMLCsHnOPdevZQjQtCncdVfYPFo6QEREJMc99RScdVas3awZFBeHywOwZQssWwY77ugLphqkveFEREQywXffwXPP+cJk4ECoE3Csp2IfQ35+mBzx6tSBtm1Dp/A0DCciIpJi69bBoYfChRfCGWfA+eeHzXP88XDZZb5A2XFHeOihsHnSjYbhREREUuyll3yBUsYM1q+HvLxwmQA2bYL69cNmSIHcG4YLvTeciIhkhlWr4P77oUEDuOIKP/wVyu67+wKprL9ip53CF0qQE4VStahnSUREst4PP0DXrvDZZ7598MEweXLYeUL33AOjRvn5Qg88AL16hcuSY5LuWVKxJCIiWW/aNDjkkPLPffWV7+GRnKONdEVEJD0sXuwnDV96KXz+edgsbdtCo0ax9i67+KEvkarI+DlLIiKSfjZsgKOP9gUTwIsvwoIF0KRJmDytWsG4cTBihJ+zdNtt6TFHSDKDhuFERKTGLVwIHTuWf+7DD6FLlzB5ROJoGE5EJFd9/DH86lfQty+8/nrYLHvsUX4+UMuWsOee4fKIbA/1LImIZIHSUthrL789BEDjxr53p3XrcJk++QT++le/bcUf/wj77Rcui0ic3FtnSURE/NYZZYUSwI8/+knVIYuljh3hscfCfb5ITdEwnIhINU2fDgce6O+0Gj06bJaddoKePWPt1q3hgAPC5RHJJhqGExGppjZtyvfmzJgB3buHy1NUBHfe6XuVrrhCc4REtiL3huG03YlIbnEOSkqgYcOwOUpLYcWK8s8tWxa2WMrP97fGi0jNyvhhuKFDh/Lxxx8zffr00FFEpJa99RbsvLNfXHDIkNi+WiHUrQuDBsXa7drBkUeGyyMitUfDcCKSMfbeG774ItZ+7jkYMCBcni1b4Omn/eTqgQP9qtAikvZybxhORGrX6tV+uKlzZ7/ycUjr1iVup1qdOnDWWWEziEjty/hhOBGpPRMn+knCXbv6TUiLi8Pm+eMfY8ddusAppwSLIiI5RMNwIrJVXbrAvHmx9ujRcOWV4fIAzJ0LK1fC4YeH22dMRDKahuFEMt3s2X7D0V69wi4oCH4Sc6J2CAccoPWDRCS1NAwnkkYef9zfen722b4gWLgwbJ7bb4emTf3x4YfD+eeHzSMiEoKKJZE0Mnq0v8MK4Ntv4dFHw+bp2xeWL/d3oL39tt9vTEQk12gYTnLeCy/AzJnQp0/4dXJ22ilxO4RmzfxDRCRXaYK35LT774dLL/XHder4u7+OPz5cni++gFNPhfnz4cQT/bBc6Nv1RUSyTO5N8NZ2J7I9nnsudrxlC4wbF7ZYat/eT/AWEZH0kfFzlrTdSWZxDu64A047DW69Nex2FQAdOiRui4iIZHzPkmSW0aPht7/1x88+C2Zw7bXh8hQW+lWgZ83yk5mvvjpcFhERSU8qliSlpkxJ3E61Jk1g7NiwGUREJL1l/DCcJLZ5MwwbBh07wplnQlFR2Dy9eiVui4iIpBv1LGW50aPhnnv88cKFkJ/v7wAL5Yor/NDbu+/CYYfB0KHhsoiIiFSFiqVasmEDNGwYOgV8/nnidgiXX+4fIiIimUDDcDVswwZ/63mjRrDXXuU3IQ3htNOgXlxJfNZZ4bKIiIhkoqSLJTPrbWYTzGy5mTkzO6UKrznSzGaa2QYz+8LMLqvknIFm9rGZlUS/Dkg2Wzq47z546SV//OWXcNVVYfMcdRS8956/6+vFF+Hii8PmERERyTTVGYZrAswBxgLPbutkM2sHvAg8AAwCDgf+aWZfO+eejZ5zKPAUcAPwPDAAeNrMjqjKCuMbNsCMGdX4k9SCdevKt4uLw+SId/DB/iEiIiLJ267tTszMAQOcc+MSnPN34CTnXOe45+4DDnTOHRptPwU0d84dH3fOS8Ba51zCgaOiIujdG+bOLQbyGTu2iCFDwm13smQJHHqo33y0Xj144gk/FCYiIiJpIentTlIxZ+lQ4JUKz70MHGRm9bdxzmFbe9OSkhKKi4sZM6Y4Wij5Lpy//rVGMldb27Ywd64f8vroIxVKIiIimS4Vd8PtBqyq8Nyq6GfvDKxIcM5uW3vTwsJCbrzxxp89n5e3PVFrxk47hd1fTERERGpOqu6GqzjWZ5U8X9k5Wx0jHD58OEVFRaxeXcTRRxcBSwH4xz+2M6mIiIhInFT0LK3k5z1EuwKbgTXbOGcV0LayN83LyyMv2o30+uvw2Wd+E9TevWsqtoiIiEhqepamAP0qPNcfmOGc27SNcyZX5QPMoGXL7cooIiIiaW7WLLjwQrjmGlizZtvn15Ske5bMrCmwT9xT7cysK/Ctc26JmRUCrZ1zg6Pfvw+40sxuxy8fcChwIXB23HvcBbxjZtcB44GTgWOAIwAtoygiIpLjli2Do4+OLckzZQpMnZqaz65Oz9JBwAfRB8Dt0eO/RNutiBs6c84tAn4JHAXMxq+ldFXZGkvRcybji6JfA3OBIcCZzrlp1cgnIiIiWWb27PJrF06bBhs3puazt2udpRSpUsDi4mLy8/MpKiqiefNw6yyJiIhki/HjYcQIf6f5HXf4DdBDWbwYCgpg/XrfPvBAX0BVQ9LrLKlYEhERkZ9Ztgz23jvWe7Pjjn7B5ZBL9Lz7Ltx1F+Tn+3UVd9+9Wm+TdLGUirvhREREpIo2bPBfGzYMm2Pp0vLDXN9+C999F/aGql69/CPVUrXOkoiIiGzDXXdB06bQpAncfnvYLAceCJ06xdpHHgm77houT0gahhMREUkDK1ZA69ZQ9mvZzPfutG4dLtOaNTB2rO/luvBCaNQoXJYalHvDcJFIhEgkQmlpaegoIiKSgSZP9kNfvXv7DdBDWb8+ViiBP/7xx3B5wG/f9bvfhc2QDtSzJCIiOWvoUPjnP/3xscfCf/8LdeuGy3P++fDII/743HPhscfCZcliuhtOxZKIiFTF2rX+Dq94U6bAIYeEyVNm6lTfq3TIIX4oTmpc7g3DiYhI5igpgXvvhaIi34uy117hsuTlQYMG5e/4Sod/a4cu1uTndDeciIikzGmn+X29Ro70RcHXX4fL0rgxjBnjJy/XqQN/+Ytf9FCkIg3DiYhISmze7Hty4n/tjB8PJ50ULhNAaal/NGgQNoekTNLDcOpZEhHJYitWwDnnQN++8PTTYbPUq+dXhC5Tty7ss8/Wz0+VunVVKElimrMkIpLFzjgD3nvPH7/1FrRvDwcdFC7PhAlw1VV+Jehrr9Wwl2QGFUsiIrXAufS4k2nOnNjxli3w0Udhi6VOneCVV8J9vkh1aBhORKQGzZnje2/q1/d3e23ZEjbPccfFjps0gSOOCJdFJFNpgreISA3q2ROmT4+1H3kEzjsvXJ4NG+DOO2HVKp+jW7dwWUTShNZZEpHcs2QJLF4MXbuGXydn7dry7W+/DZOjTMOG8Ic/hM0gkukyfhguEolQUFBAjx49QkcRkQD++1/o0MHviH7AAbB8edg811wTm6vUpo2fYC0imU3DcCKS0Q491G8PUWbkSBgxIlgcAGbOhC+/9AXcTjuFzSIiP6NhOBGpfePGwQcf+LV7evcOm6VRo8TtELp39w8RyQ7qWRKRpNx3H1x+uT+uUwdefNHv1h7K7Nlw/PGwcqW/02vSJGjaNFweEUl7WsFbRGrXc8/Fjrds8dtVhNS1KyxbBt98A+++q0JJRGqeiiWRNLd5M/z+93D44f7rpk1h83TokLgdQt26mhskIrVHc5ZE0tzNN8Mtt/jjyZP9rfF/+lPYPN9/D7NmQb9+MGxYuCwiIqmgYklkK378ERo3Dp0CPvywfPujj8LkKNO0KTz0UNgMIiKppGE4kQq++84PeTVp4vexWrQobJ4TTijfPv74MDlERHKVepZEKrj1Vj/cBfDJJ37146eeCpdn8GBfuL3/vi/iBg4Ml0VEJBepWJK0sGoVTJkC++4L++0XNktxceJ2CAMHqkgSEQkl44fhtN1J5lu0yG9TMWAAHHhg2F4c8GsIld1Z1bAh/O53YfOIiEhYWpRSghsxAv7yl1i7Z0+YNi1cHvA9XbNmQefOsNdeYbOIiEiN0nYnUjVvvglPPukLgWuvhQYNwmXZYYfE7RBattREahER8VQs5aAZM6B/f7/YIcDnn8ODD4bLc/nlvnibMMHPWbr77nBZREREKlKxlIPefjtWKAG8/nq4LAB5eX7LjNJSvxKziIhIOsn4Cd6Z4uabYc89/XycefPCZunWLXE7FBVKIiKSjjTBOwXeeAP69o21998f5s4NlwfgscfgP//xc5Zuvhny88PmERERSRFN8C6zZQt8/LHfR6tt27BZli4t316yJEyOeIMG+YeIiIgklpXDcKWlcPLJvgenXTsYPTpsnv79YbfdYu3zzw+XRURERJKTlcNwU6c259hjY9/Ly/ObotYJWBp+9RWMG+dvSR84ECzpTkARERGpARqGg59PFA5ZJJVp3RqGDg2dQkRERJKVBmVEzevTB84+2x/XqweRSHoUTCIiIpJ5Mn4YLhKJEIlEKC0tZeHCheXuhlu61O/WvuOOKckpIiIi6S/pYbiML5bKpPPSASIiIpI2ki6WNDglIiIikoCKJREREUnK11/D8uWhU6SOiiURERGpsttu88vgtG4Nw4aFTpMamrMkIiIiVVJcDDvs4HfJKDNnDhxwQLhM1aB1lkRERLLJ5s3w6KNQVOSXxWnZMlyWLVvKF0rg82W7ag3DmdkVZrbIzDaY2Uwz67WN84ea2XwzW29mn5jZ4ArfH2JmrpJHw+rkExERyRbnnAMXXADXXAMHHwxr14bL0qIFjBgRaw8eDN26BYuTMkkPw5nZmcCjwBXA+8ClwEVAgXPuZ1vEmtnlwN+Bi4HpQE/gAeAc59yE6DlDgLuAjvGvdc6tRMNwIiKSozZuhIYNIf5X9cSJcMIJ4TIBfPEFbNgABQVhc1RTSpYO+C0wxjn3oHNuvnPuamApcPlWzj8P+Jdz7inn3BfOuSeBMcB1Fc5zzrmV8Y9qZBMREdkuS5bAgAFw2GF++CukBg1g991j7Tp1oG3bcHnKtG+fsYVStSQ1Z8nMGgDdgZsrfOsV4LCtvCwP2FDhufVATzOr75zbFH2uqZl9CdQFZgM3OOc+2FqWkpISSkpKfmoXFxdX9Y8hIiKyVaefDv/7nz+eOhU6dYIePcLlGT8eLrvMz1m67jrYf/9wWXJVsj1LO+OLmVUVnl8F7LaV17wMXGRm3c07CLgAqB99P4AFwBDgJOBsfHH1vpntu7UghYWF5Ofn//Ro06ZNkn8UERFJF99/77eoqjh5OIR582LHzsH8+eGyAHTvDtOnw8KFcOGFYbPkququs1RxHpFV8lyZvwKTgKnAJmA88FD0e6UAzrmpzrnHnHNznHPvAmcAC4GtruAwfPhwioqKfnosXbq0mn8UEREJ6bXXoFUrP7x0zDF+LkxIJ54YO27eHHr3DpdF0kOySwd8gy9wKvYi7crPe5sAcM6tBy4ws0uBlsAK4BJgXfT9KnvNFjObDmy1ZykvL4+8vLwk44uISLr5zW98zxLAm2/CY4/BRReFy/Pww/6us1WrYNAg2GuvcFkkPSRVLDnnNprZTKAf8Hzct/rhe4wSvXYTsAzAzM4CJjrnKu1wNTMDugIfJpNPRESqZvx4+OQTOO648AsKbtqUuJ1qDRrA1VeHzSDppTrDcLfj5yBdYGadzewOoC1wH4CZFZrZI2Unm1kHMxtkZvuaWU8zexLoAlwfd84IMzvWzNqbWVf83XJdy95TRERqzi23wCmn+MnCBx8MH2z1VprUKCz0BQr4NXsGDQqbR6SipFfwds49ZWY7AX8GWgEfAb90zn0ZPaUVvngqUxe4Fr+G0ibgTeAw59ziuHNaAPfjh/eKgA+A3s65/yWbT0REEnv88djxhg0wbhz84hfB4jBwICxa5Ie9CgpAMywk3WhvOBGRWrZund9wdPZs6N8fbr7Zr5cTysCB8NxzsfbYsTBkSLA4IqmmveFERNLNddf5ScPgNx1t0ybsbu3//KdfGfqTT/xw3Pnnh8sikglULIlIVvr+e1i8GNq1gyZNwmZZuLB8+5NPwuQo07IlTJgQNoNIJgnYESwiUjsWLIAOHfxKxx07wqefhs0zYEDsuE4d35sjIplDPUsiknVGjYIVK/zxV1/5OUJjxoTLM3So399r9mzo21eLHIpkGhVLIlIjXnkFXn/d3/p95plhs5glbocwYED5HiYRyRwZXyxFIhEikQilpaWho4jkrAkT4OST/T5aAF9/DVdeGS7PH/8Ib7wBy5b5LTSuv37brxER2RotHSAi2+3SS+H++2PtY4+Fl14Klwdg/XpYssQXS40ahc0iImkl6b5mTfAWyUClpXDZZX7z0T59YvNzQikoKN/u3DlMjniNGvnJ3SqURGR7ZfwwnEguuv9++Ne//PHKlX7I69lnw+UZNsyvvlw2Z2nUqHBZRERqmoolkSoqKoJp0/yCgqF7TpYtS9xOtTp1VCCJSPbSMJxIFaxe7XtMjj0WunSBRx7Z9mtq0xlnQOPGsba2qhARqT2a4C1SBXfeCddcE2t37OgXPgxpwQJ/x1enTn7ekoiIVIn2hpPs8eqr8O9/w267wciRkJ8fLkuzZonbIXTq5B8iIlK7VCxJWpo7F044ATZt8u2FC+G//w2XZ/BgmDgRxo3zxds//xkui4iIpJaKJSln/XqoXx/qBf4vY/r0WKEE8P774bKA/5k8/zz8+KO/FT0dVoQWEZHU0ARv+cmwYX7ScIsWMH582Cw9evgCpczhh4fLEq9xYxVKIiK5JuMneMdvd7Jw4UJN8K6mt96Co4+OtVu0gLVrQ6XxXn0Vxo71w14jRoSdsyQiIlkj6X/yZnyxVCYT74Zbuxb+8x8/rHPeedCgQbgsEyfCiSfG2vXrw4YNfv0cERGRLKK74TLFDz/4oaX58337uefCTmDu1w8OOwwmT/btG25QoSQiIgLqWQrm3Xehd+/yz61ZAzvuGCYPwMaNfiJ1ixbwi1+EyyEiIlKLtJFuIi+/7IuAbt38Yn4htW5d/o6zHXYIv3ZPgwZ+3pIKJRERkZicKZbWrIFTT4XZs+GDD+CUU/xeX6G0bw8PPwz77gsHHAAvvFD+7i8RERFJDzkzDDdvnt/TK97Chb5YERERkZyhYbit2Xff8sNLPXtCu3bh8oiIiEhmyJm74Ro08GsJ/fvfflHBCy8Mv0q1iIiIpL+cGYYTERERQcNwIiIiIjVLxZKIiIhstx9/hNWrQ6eoHRlfLEUiEQoKCujRo0foKCIiIjlp/HjYeWdo2RLOPBO2bAmdqGZpzpKIiIhsl912g1WrYu3x4+Gkk8Ll2QbtDSciIpLttmyBBx+EJUvgtNOga9eweUpKErczXcYPw4mIiOSaq66CSy+Fm24qvyl7KIWFsc3XDz8cTjwxbJ6apmJJREQkw4wfHzv+8Ud47bVwWQAuu8zvijF1Krz5JjRsGDZPTVOxJCIisg2ff+43YW/aFAYPhs2bw+bp1ClxO4S994aDD87OfU5VLImIiGzDFVf4Tdh/+AEefRTGjAmb55FH/ObwBx0Eo0dDv35h82Q7TfAWEZG0NG0azJ0LvXtDx45hs1RcPyj+zq8QWrWCZ58NmyGXqGdJRETSzhNPwKGHwiWX+E3QZ8wIm2foUL+vKMCOO8I554TNI6mlniUREUk7DzwAZcsArl8Pjz3mh5xCuegi6NIFPv0UjjoK2rQJl0VSTz1LIiLC+vXw619D586+Nyf0Ojm7716+3apVmBzxDjkEzjtPhVIuyviepUgkQiQSobS0NHQUEZGMNXIkPPSQP16wwBcrI0eGy3PbbbB8OcyZA8ceC1dfHS6LiLY7EREJZMUKmDUL9tsP9torbJaBA+G552LtQYP8XV8iWSjp7U40DCciEsCHH0JBAfzqV/7rm2+GzXP66bFjM7+Fhoh4GT8MJyKSie69F777zh+vXw933AFHHx0uz1ln+V3jp03z21UcdVS4LCLpRsWSiOSM226DCRP8sNc//gFNmoTL0qxZ+XZ+fpgc8Y45xj9EpDzNWRKRnPCf//h5OGUuuQT+9a9wedau9UNwkyf74m3SJN1lJZIiSc9ZUs+SiNSa9ev93JzWrf0jpLlzy7fnzAmTo8wOO8D778OGDdm36ahIttEEbxGpFWvXQo8efmPNvfeGcePC5unXL7YCM/jb0dOBCiWR9FetYsnMrjCzRWa2wcxmmlmvbZw/1Mzmm9l6M/vEzAZXcs5AM/vYzEqiXwdUJ5uIpIeHH4Z58/xxSQn86U9h8xxzjB/qGjYM7rsv7BpCIpJZkh6GM7MzgTuBK4D3gUuBSWZW4JxbUsn5lwOFwMXAdKAn8ICZrXXOTYiecyjwFHAD8DwwAHjazI7IgDlVImnjxRf9bui77w5//Su0aBEuS4MGidshHHts+vQoiUjmSHqCt5lNA2Y55y6Pe24+MM45N7yS8ycD7zvn/l/cc3cCBznnjoi2nwKaO+eOjzvnJWCtc+6squTSBG/JdR98AD17wubNvv3LX8J//xsuz/r1cPzx8PbbvmgbP97vHi8iEljtTvA2swZAd+DmCt96BThsKy/LAzZUeG490NPM6jvnNgGHAndUOOdl4OqtZSkpKaEkbvOi4uLibcUXqXHOwWefQfPm0LJl2CwzZsQKJYCpU8NlAWjUyC+0uGKFn8zcqFHYPCIi1ZXsnKWdgbrAqgrPrwJ228prXgYuMrPu5h0EXADUj74f0dcm854UFhaSn5//06ON7rmVFCsthVNPhQ4d/J1eY8aEzXPIIVC/fqzdK+FMwtQw80OCKpREJJNV9264imN3VslzZf4KTAKmApuA8cBD0e/F736bzHsyfPhwioqKfnosXbq0itFFasarr8bu8Cothd/8Jmgc9t8fXn4ZBg+G666Dxx4Lm0dEJFskO8H7G3yBU7HHZ1d+3jMEgHNuPXCBmV0KtARWAJcA66LvB7AywXu2rex98/LyyMvLSzK+ZLq5c+HOO6FxY3931W5b7Xusfel478HRR4fdMkNEJBsl1bPknNsIzAT6VfhWP2DyNl67yTm3zDlXCpwFTHTObYl+e0ol79l/W+8puWX1al8IjB0LkYhfNydkwdK/P5x8sj+uW9cXcSIikn2qs4L37cCjZjYDX+Rcgu/9uQ/AzAqB1s65wdF2B/xyAdOAHYDfAl2A8+Pe8y7gHTO7Dj9MdzJwDHAEvrCSgFasgDp1wk9gnjcPvv021v7oI78R6Q47hMlTty48/3z6TPAWEZHakfScJefcU/i71P4MzAZ6A790zn0ZPaUV5YfO6gLXAnOAV4GGwGHOucVx7zkZXxT9GpgLDAHOdM5NSzaf1Kw//MFP0N1tN79uT0idO5ffbLRjx7DrCIGfwLzvviqURESymTbSla36/HPYZ5/yz61YEXae0IwZcOutfrf4kSO18aiIiCRNG+lmuq++gsJC2LQJrr3W35YeSmlp1Z5LpYMOgiefDJtBRERyi4qlNLJ5M/TpAwsX+vYLL8CCBeWHnlKpQwe4/HK4917f/t3vwu8cLyIikmoahsNPGl6/PnwhsGzZz4eV/vc/v3N7SAsW+MnM++4bNoeIiEgNSHoYrrqLUmaN+++HXXeFPfaACy4Im2XXXaFt3NT4HXeEvfcOl6dMp04qlEREJHdlfLEUiUQoKCigRzW6XzZvhmHDYvNwxo6F99+v4YBJaNAAXn8dzjkHTjsNXnvNF0wiIiISTk4Pw23a5Pesip+0/NZbcOSR1YkpIiIiGUDDcMmoX9/feWbRH9uAAemx+aiIiIikj5zuWSqzaBF8/z106RIrnERERCQraZ2l6mjXLnQCERERSVc5PQwnIiIisi0qlkRERKTGbdgAw4f7u7ufeCJ0mu2jYTgRERGpccOGwYMP+uNnn4WddoL+/cNmqi71LImIiGSJb76BmTPhhx9CJ4EpU8q3p04Nk6MmqFgSERHJAu+/73d9OOggf3f3smVh88QvxWMGhx8eLsv20jCciIhIFrjxRigu9seLF8M998DNN4fLc9ddsNtu8OmncOqp0LdvuCzbS8WSiIhINWzeDCNG+A3Pe/WCP/0J6gQcr6lfv3y7XuDf8A0a+J9PNsj4YikSiRCJRCiN37NERESklt10E4wa5Y9few2aNYNrrgmXZ9QoP19p1SrYf3+4+upwWbKNVvAWEZGMsXo1zJkDnTvDHnuEzXLKKTB+fKw9aBA8+miwOIC/XX/1ath99/A9S2lMe8OJiEh2mjfPF0n9+0OnTn5Cc0j9+iVuh9CwIbRtq0KppunHKSIiGWH0aPj2W3/8ww9w661h77AaOhSaNIFp06B3bzj77HBZpHapWBIRkUo5BzfcABMn+lvRIxHIzw+Xp0mT8u2mTcPkiDdkiH9IdtOcJRERqdSDD8LFF8fav/41/Pvf4fJ88w0cd5yfxNyxI7z8Muy5Z7g8krGSnrOkniURkTSyfDm88Qa0bw+HHRY2yyeflG8vWBAmR5mdd4YZM/xaQvo3saSSJniLiKSJxYuha1c47zw/F+ef/wyb58QToW7dWHvAgHBZ4qlQklRTz5KISJp4+mn4+utY+5574IorwuXp3RveftsPd3XpAmecES6LSEhZUyz9+KP/Ont2ekz6ExFJVtn/x8o0agSzZoXJEp/hlFP8cegsItvSqRM0blzz75s1E7zfeaeYI4/MB4oA9dGKiIjkmpkzoVu3bZ6WexO8y7Y72bTJb4rz9tvqWRIREclFnTrVzvtmTc+Slg4QERGRKtB2JyIiIiI1ScWSiIiISAIqlkREREQSyIQ5S1ViZs3xt8LlO+eKQ+cRERGR7JBNxZIBzYB1Llv+UCIiIhJc1hRLIiIiIrVBc5ZEREREElCxJCIiIpKAiiURERGRBFQsiYiIiCSgYklEREQkARVLIiIiIgmoWBIRERFJQMWSiIiISAIqlkREREQSULEkIiIikkC90AESidvvTURERKSmJLWPbFoXS/hCqSh0CBEREckq+UBxVU9O6410k+xZagYsA/YA1tVaKPgf0LMW3z/bPidV1wWy52eWqs/ItmuTTddf/z9Lz8/Jtr8zqfqcdPw7kz09S9E/SJUqP19XAf4HUOVqMVlmtqU23z/bPidV1yX6WVnxM0vhZ5QdZsW1ybLrX3ao/5+l0edk29+ZVH1ONvyd0QTv5EX0OWkrm35m2XRdILt+Ztl0bbLtZ6Zrk56fk/HXJa2H4ZJhZs3x85vyU1GNS9XouqQvXZv0pWuTnnRd0ldtX5ts6lkqAW6MfpX0oeuSvnRt0peuTXrSdUlftXptsqZnSURERKQ2ZFPPkoiIiEiNU7EkIiIikoCKJREREZEEVCyJiIiIJKBiSURERCSBrCiWzOwKM1tkZhvMbKaZ9QqdKdeY2XAzm25m68xstZmNM7OOFc4xMxtpZsvNbL2ZvWVm+4XKnIui18mZ2Z1xz+m6BGJmrc3sMTNbY2Y/mtlsM+se931dmxQzs3pm9rfo75T1ZvaFmf3ZzOrEnaPrkgJm1tvMJkR/zs7MTqnw/W1eBzPLM7PRZvaNmf1gZi+Y2R7JZsn4YsnMzgTuBG4CfgG8C0wys7Yhc+WgI/GrtB4C9MNvpfOKmTWJO+f3wG+BK4EewErgVTOr6v5/sh3MrAdwCTC3wrd0XQIwsx2A94FNwPFAAXAt8F3cabo2qXcdcBn+Z94Zfw3+HzAs7hxdl9RoAszB/5wrU5XrcCcwADgLOAJoCkw0s7pJJXHOZfQDmAbcW+G5+UBh6Gy5/AB2ARzQO9o2YAVwXdw5efhfDJeGzpvtj+j/IBYCxwBvAXfqugS/JjcD7yb4vq5NmOsyERhT4blngUd1XYJeFwecEtfe5nUA8oGNwJlx5+wOlALHJvP5Gd2zZGYNgO7AKxW+9QpwWOoTSZz86Ndvo1/bAbsRd62ccyXA2+hapUIE+K9z7rUKz+u6hHMSMMPM/i86dP2BmV0c931dmzDeA/qaWQcAMzsQ3yPxYvT7ui7poSrXoTtQv8I5y4GPSPJa1dvOsKHtDNQFVlV4fhX+hygBmN/++XbgPefcR9Gny65HZddqz1Rly0VmdhbQDd9NXZGuSzjtgcvxf1dGAT2Bu82sxDn3CLo2ofwd/4+9BWZWiv8d80fn3BPR7+u6pIeqXIfdgI3OubWVnJNUjZDpxVKZinu2WCXPSercAxyA/9dYRbpWKWRmbYC7gP7OuQ0JTtV1Sb06wAzn3PXR9gfRyamXA4/Enadrk1pnAoOAc4B5QFfgTjNb7px7OO48XZf0UJ3rkPS1yuhhOOAb/NhjxQpxV35ebUoKmNlo/PDC0c65ZXHfWhn9qmuVWt3xP+OZZrbZzDbjJ+NfFT0u+9nruqTeCuDjCs/NB8puTtHfmTBuAW52zj3pnPvQOfcocAcwPPp9XZf0UJXrsBJoEL2ZYmvnVElGF0vOuY3ATPzdV/H6AZNTnyh3RW/hvAc4FejjnFtU4ZRF+P9w+8W9pgH+F7euVe15Hdgf/6/jsscM4D/R4y/QdQnlfaBjhec6AF9Gj/V3JozGwJYKz5US+32p65IeqnIdZuLvNo0/pxXQhSSvVTYMw90OPGpmM4Ap+Fuj2wL3BU2VeyL4buuTgXVmVlbtFznn1jvnytb2ud7MPgU+Ba4HfgQeDxE4Fzjn1uEnM/7EzH4A1pTNJ9N1CeYOYLKZXQ88jZ+zdEn0gf7OBDMB+KOZLcEPw/0Cf3v6v0HXJZXMrCmwT9xT7cysK/Ctc27Jtq6Dc67IzMYAt5nZGvwNR7cCHwIVb3ZJLPTtgDV0S+EVwGKgBF9J9g6dKdce+PHfyh5D4s4xYCR++GED/q6FLqGz59qDuKUDdF2CX4tfRf/HvQE/BHdxhe/r2qT+mjTDr83zJbAe+Bz4G9BA1yXl1+Korfxeeaiq1wFoCIwG1uALqQlAm2SzWPTNRERERKQSGT1nSURERKS2qVgSERERSUDFkoiIiEgCKpZEREREElCxJCIiIpKAiiURERGRBFQsiYiIiCSgYklEREQkARVLIiIiIgmoWBIRERFJQMWSiIiISAL/H7qyzZRw/FjUAAAAAElFTkSuQmCC\n", "text/plain": [ "Graphics object consisting of 102 graphics primitives" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Get asymptotic formula from previous calculation in \"(pi, 1)-direction\"\n", "ASM = (ex*pw*se).subs(r=pi,s=1)\n", "\n", "# Define oscillating factor\n", "def M(n): return (pi/(pi+1))^(-round(n*pi)+n*pi)\n", "\n", "# Plot oscillating factor, together with minimum and maximum\n", "pt = plot(sqrt(pi/(pi+1)),n,0,100)\n", "pt += plot(sqrt((pi+1)/pi),n,0,100)\n", "for k in range(100):\n", " pt += point([k,M(k)])\n", "pt" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAGGCAYAAABmPbWyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJrElEQVR4nO3dd7hU1bnH8e9L74gFERARsaJYUDCIolHEa8NyjSWAifGqgFGDsTdiCZZEMVcsUaNXE6+a2AuKRo0FsaBgiIANUIoErCD98N4/1p579hlOmXOYPfucmd/neeY5e+1Zs9c7W+S8rLX2WubuiIiIiEj1GqUdgIiIiEhDoKRJREREJAdKmkRERERyoKRJREREJAdKmkRERERyoKRJREREJAdKmkRERERyoKRJREREJAdKmkRERERyoKRJREREJAd1SprMbKSZzTazlWY2xcz2raH+KDObYWYrzGyWmQ3Per+pmV1uZp9G15xmZodk1WlrZuPMbG50nUlmtlclbe1oZk+a2XdmttTMJptZt7p8TxEREZGMWidNZnY8MA64BtgdeA2YUFViYmYjgLHAGKAXcAUw3syOiFW7Gjgd+CWwE3A78JiZ7R6rcxcwCBgG7AJMBF40sy6xtrYBXgdmAvsDuwJXAStr+z1FRERE4qy2G/aa2VvAe+4+InZuBvC4u19USf1JwBvufl7s3DhgT3cfEJUXANe4+/hYnceBZe4+1MxaAkuBIe7+TKzOVOBpd780Kj8IrHH3YbX6UiIiIiI1qFVPk5k1A/oQenniJgL9q/hYc9bv6VkB9DWzpjXUGRAdNwEaV1fHzBoBhwEfmdnzZvZvM3vLzI6q5vs0N7N2Wa9Nzcyq+oyIiIiUpia1rL8pIXlZlHV+EdCpis88D5wa9Ry9R0i6TgGaRtdbGNUZbWavAp8CBwJDorZw96Vm9iZwWdSrtQg4EegHfBy10xFoA1wIXApcABwCPGpmB7j7PyqJ7SLCcGG29sD3Vd6FytWuy05EREQKJS+dIXV9ei47QbBKzmVcBUwAJgNrgCeAe6P3yqKfZxOSn5nAauAW4J7Y+xDmMhkwH1gFnAU8EKuT+S5PuPtN7j7V3a8FngbOqCK2sYQEKfPqWkU9ERERKXG1TZqWEJKU7F6ljqzf+wSAu69w91OAVkB3oBswhzBHaUlUZ7G7HwW0BrYCdgCWAbNj1/nU3QcSepO2dPe+hN6qTJ0lwFrgw6wQZkRtVhbbKnf/PvOKYhIRERFZT62SJndfDUwhPMUWNwiYVMNn17j7PHcvA04gTOBel1VnpbvPJwwbHkvolcq+zg/uvtDMOgCDM3Wi2N4Bts/6yHbA3By/ooiIiEilajunCeBG4H4zexd4EziN0JNzO4CZjQW6uPvwqLwd0Bd4C+gAjAZ2Bk7OXNDM+gFdgKnRzzGEhO76WJ3BhOG5WUBP4Ibo+J5YbDcAD0Vzo14mzGk6grD8gIiIiBTI55/DuHHQqBH8+tfQqaqZzw1IrZMmd3/IzDYBLge2AKYDh7p7pjdnCyoOhzUGziX0AK0hJDP93X1OrE4LwlpNPQjDcs8Cw9z921id9oQ5SF2Br4FHgEvcfU0stsfM7AzCBO8/EJKqY9399dp+TxEREambH36A/faDuVFm8PTT8MEH0KxZunFtqFqv01TMzKwd8B3QPprjVBu6kSIiIsC0abDbbhXPffIJbLNNKuFAyk/PiYiIiFRqq61g443Ly506QefO6cWTL0qaREREZINMmwa9e8Nmm8Gll8JGG8HEiXDEEXDUUfDCC9CyZX7bLCuDzz6D72s7LrQBNDwXo+E5ERGR2ttlF5g+vbw8YQIcckhy7S1fDoMHw+uvQ5s28OijMCj7uf6KNDyXL2Y2ysw+BN5OOxYREZGGZsGCiuWFC5Nt7777QsIEsGwZjB6dbHsZSpoAdx/v7jsRlkYQERFp0D77DG66CR5+uDDtnX56+XGXLnDoocm2V1ZWfTkpGp6L0fCciIg0dHPmQJ8+8PXXoXzuufC73yXf7oQJoYfpsMNg882TbWvpUth/f3jvPWjRAh56CI48stqP5GV4TklTjJImERFp6G67DUaOLC936pT8cFkaVq+GmTPD9+vYscbqeUma6rIiuIiIiNRTW29dsdy9eyphJK5Zs/DEXiFpTpOIiEgD9sILYSHJ3XYLj/kfcghcc01YSHK//eDPf85/m889B2PGwIsv5v/a9ZmG52I0PCciIg3JN9/AlluGbUsAWreGL76ADh2Sa/OBB+CnPw3HZvC3v8ExxyTXXp5oyQEREZFStmhRecIE4XjRomTbfOSR8mP3sEZSqVDSJCIi0kD17Al77lle3nPP5Pd322676svFTMNzhMUtgVGEJHJ7NDwnIiJ1UFYGV14Jr70Ge+8djpsk/MjV0qVwzz3h+Oc/h7Ztk21v5Uo46yyYPBn23TesB9WsWbJtAnz7LTRqBO3a1enjWnIg3zSnSURENsTvfgfnnVdevvJKuOyy9OIpFpdfDlddFZKm3/8ezjmn1pfQnCYREZH6ZOrU6stSe598EhImgHXrwmKdX32VTixKmkREROpozRqYMiU8sQZw8MEV388u58Prr4cVsNNKHApt1aqK5XXrwsKWadDiliIiInWwYgUceCC8+WaYt3T33TB8ODRvXj6naejQ/LZ5/fVwwQXhuFs3eOednFbDbtB69QpLHPzlL6F81lmwxRbpxKI5TTGa0yQiIrn6859h2LDyciG2K+nSBRYsKC/fcQecdlqybdYX770XktM6rgKubVRERETS0rRp9eUkbLppxaRp002Tb7O+2GOPtCPQnCYREZE6OfZYOPzwcNyqFYwfn3yb99wT1mFq0SJsynv00cm3uXw5vPpqmJBd6jQ8F6PhORGRhu2ll2D2bBg8GLp2Tb4999Dz0749tGmTfHuF9t13MGAATJ8OjRuHeVsnn5x2VHWidZryRYtbiog0fPFJ0ptuCu++C1ttlW5MDd0dd8AZZ5SXu3cPSWkDpHWa8sXdx7v7TkDftGMREZG6ueOO8uMlS0prT7SktGpVfbnUKGkSEZGi0LlzxXK+H0ufNy887n7mmQ22t6XWTjwRjjwyHG+0Edx6a6rhpE7DczGa0yQi0jAsXx7WRHrllbBJ7QMPhN6lE0+EOXPCuj433wyWl0GZsJhir17lk6G33BI+/LA45zFV5ttvw3dNeh+9BGlOU74paRIRaRgye5FlnHEG3HZbcu19+in07Fnx3NSpsOuuybUpeaU5TSIiUpriaxVVVs63zp0rDv9tummYFC2lRUmTiIg0OMOHQ7Nm4bhxYzjllGTba9kSXnwRfvKTsD7Tiy+GZQaStHx5GGbs3h1OOgl++CHZ9uqTd96Bp5+GpUvTjqQiDc/FaHhORKThmD4dJk0KK0XvuWfa0eTfRRfBtdeWl88/H667Lr14CuW66+DCC8PxjjuGvf3ykKBqGxURESldO+8cXsVqzpzqy8UqnhjOmAFPPZX/jY/rSsNzIiIi9dBJJ0Gj6Ld0o0bhycBSsNFG1ZfTpOE5tCK4iEh9t3o1TJkSJmBvu23a0RTOpElheGrvvWGffdKOpjBefRWOOw6++gpOPTU8FZmHpSO05EC+aU6TiEj9s2IF/PjHMHly6HG55RYYMSLtqCRpa9fmdV0oLTkgIiKF99lnoedjs83C6thJ/9v7ySdDwgSwbh1cckmy7Un9UB8X0lTSJCIitfJf/wVvvRVW4B4/Hv7852Tba9Gi+rLk17//DYsXpx1F/VSnpMnMRprZbDNbaWZTzGzfGuqPMrMZZrbCzGaZ2fCs95ua2eVm9ml0zWlmdkhWnbZmNs7M5kbXmWRme2XVudfMPOs1uS7fUUREKlfohSWPOAKOPz4ct24Nf/xjsu1BWAG8f3/o2jU8+l8qLrsMNt8cOnasuOK6BLWe02RmxwP3AyOBN4DTgVOBndz980rqjwCuA/4LeAfoC9wJnOTuT0V1rgOGRnVmAoOBG4H+7v5+VOchYGdgBLAgqv+rqN35UZ17gc2Bn8dCWO3uX+f43TSnSUSkBjfeCOeeG447dIC3315/i5EkLFkS9j8rRE/TgAHwxhvl5UcegWOOSb7dNM2ZA1tvXfHc/Pnrb4TcQKUzEdzM3gLec/cRsXMzgMfdfb183MwmAW+4+3mxc+OAPd19QFReAFzj7uNjdR4Hlrn7UDNrCSwFhrj7M7E6U4Gn3f3SqHwvsJG7H1WrL1V+PSVNIiI5+PvfQ2/MwQcX53YiW20Fn8e6AW6+Gc46K714CuGTT9Z/MnHuXOjWLZ148qzwE8HNrBnQB5iY9dZEoH8VH2sOrMw6twLoa2ZNa6gzIDpuAjSuoU7G/mb2bzP7yMzuNLOOVX0fERGpmwMPhNNOK86ECeDnsfGKjTcOQ4TFrmfPsPFxxllnFU3ClDe16mkys87AfGAfd58UO38xcLK7b1/JZ35LGC47HHiPkHQ9A3QEOrv7QjN7ANgVOAr4FDgQeAJo7O7No+tMAlYDJwGLgBOB+4CPM+1GQ4fLgLnA1sBVhISrj7uvqiS25oSELaMtMA/1NImIlLynnw49LYcdVrzJYWU+/DAs7bDDDmlHklepbqOSnSBYJecyrgI6AZOjeouAe4HzgbKoztmEeU4zo+t8CtxDxblJw4A/EZK2MkIC9gCwx/8H5f5QrP50M3uXkEAdBjxaSWwXAVdU+S1FRIS1a8Pclk6doHnzmusXi8MPTzuCdOy0U9oR1F+1fXpuCSFh6ZR1viMhGVqPu69w91OAVkB3oBswhzBHaUlUZ3E0D6k1sBWwA6HHaHbsOp+6+0CgDbClu/cFmsbrVNL2QkLSVNX6sWOB9rFX16quJSJSihYtgt69Q09Ljx6hF0KkVNUqaXL31cAUYFDWW4OASet/osJn17j7PHcvA04gTOBel1VnZfQkXBPgWMIQXfZ1foiG9DoQnrJbr06GmW0CbAksrCKmVe7+feZFSOREROq1774LT5IVwu9+FzZNhbC0wGWXFabdUrR4Mdx/P7z0UtqRSFXqsk7TjcCpZnaKme1oZjcReo9uBzCzsWZ2X6aymW1nZkPNbFsz62tmDxKWDrg4VqefmR1jZj2iNZ+ei2K7PlZnsJkdYmZbm9kg4GVgFmEYDzNrY2a/M7MfmVl3M9sfeIrQm/VYHb6niEi9c9ttsMkmYTXuc85Jvr01a6ovJ2XKlLAKeKns9LVoEfTpA8OHh0n2V16ZdkRSmVonTdG8oXOAy4GpwH7Aoe4+N6qyBSGJymgMnAtMA14AWhDWX5oTq9MCuBr4kJDgzAcGuPu3sTrtgfGEeU/3Aa8DB7t75n/hMmAXQs/TR8D/RD9/5O7qQRKRBm/5cvjlL6Esmg16883w/vvJtvmrX8GWW4bjjTYqTE/TL38Je+4JP/oR/OQnpZE4PfUUfPFFeXn8+KrrSnq0YW+M1mkSkfps6VJo167iucmToV+/5Nv96KOw8OHGGyfb1qJFYcJ53Pvvw267Jdtu2p56Co48sry8yy7wwQfpxVOEtGGviEgpadu2Yk/Pf/4n9O1bmHb79Ek+YYLwdF7jxhXPtWqVfLtpO+IIGD06JMU77AD33VfzZ4rB2rVw6aUwaBBcfXXYkLk+U09TjHqaRKQhmDkzDNXtvjtYXv79XL/cdlsYolu3Dq64IrykOP3mNzBmTHn5ppsSm6uX6jpNIiKSkiJbdHA9I0bAz34W5m61aZN2NJKk996rWE56jt6G0vAcYGajzOxD4O20YxEREWjZUglTKRiUtYDRQQelE0euNDwXo+E5EalvvvoKFi4MG6mW0mrcaXAvzuHO+u6ee+Dtt2G//eDEExNrJi//ZZU0xShpEpH65O9/hyFD4IcfwtNU//gHdOiQdlTF59tv4Zhj4NVXw8T6xx+Hjtrqvdjo6TkRkWJ20UUhYQL45z/hzjvTjadYXXMNvPxymEP15pvhaS6RyihpEhHJ0dKl8NBD8NxzhWmvUaPqy0m4887Qq7X//uEpvVLw1VfVl4tZWRksW5Z2FA2HkiYRkRz88APssw+ccAL8x3/AqFHJt3nDDdC+fTju0wdOOy3Z9t59F04/HaZPD0OBxxyTbHv1xemnl68F1bw5jByZbjyF8sorYTuetm3DXKL6vkZSfaA5TTGa0yQiVXnmGTj88PJy48awciU0SXjhlqVLwyrZ3bsn39bDD8Pxx5eXmzWDVauSbbO++OyzsN/drrvCdtulHU1hbLcdfPxxefnhh+G449KLJ2Fap0lEpFA226xiuUOH5JMYCL0Abdsm3w6Ep5c23zwkaVDUv0DX06NHeJWS7GE5DdPVTMNzIiI56Ns3TBhu2zZsYPvww2lHlH+dOsFbb8HYsWFu0733ph2RJOnyy8uXWNh5Zzj22HTjaQg0PBej4TkRESklH34IX34Je+9d9Hv8aZ2mfDGzUcAoQs/b9ihpEhERKSZapylf3H28u+8EFGC/cBERyZgzJ8yl2nJLOP/8tKMRqZ56mmI0PCcicW++GZYaGDgQmjZNO5ridOCB8NJL5eUHH6z4BJ9InqinSUQkKaNHQ//+YUPRgw+GNWvSjqg4ffFF9eVi9c474c/WQQeF5FwaBvU0xainSUQAVqxYf1Lsyy+HVbKT9MwzMHEi9O4Nv/hFsm3VF9dcU75tSfv2IZnYdtt0Y0rasmWw1Vbw9deh3L59GKbcaKM0oyp6WqdJRCQJTZtCy5Yhecpo1y7ZNrMXz1y8GC68MNk264NLLoHddoPZs8NK69tsk3ZEyVuwoDxhAvjuO5g/X0lTQ6DhORFpEF5+GX79a/jTn5Jvq0kT+J//gdatw35vl1wCe+yRbJvPP1+xXKj97eqDww6DM88sjYQJwuruvXqVl7fbrnS+e0OnniYRqfdeeSXM/cjsjTVvXliYL0nHHRcW+1u7NmwnkrRdd61Y7t07+TYlHc2ahT/Tf/gDuIeEsUWLtKOSXGhOU4zmNInUTxdcANdfX17u2zesXF1sxo4NPUy77grXXReGCEUkLzSnKV+yFrcUkXpm552rLxeLiy4KL5Fi9e238MYbYV2uhtibqp6mGPU0idRfY8fC00/DjjvCTTcVbhPbUuBevgeZSFIWLw7btXz2WfjzdtttcPrpBWte26jkm5ImESklf/sbnHoqrFoVHv0fPTrtiKSYjR8f5m9lbL11SKAKRItbiohI3SxfDsOGhcfdV64MTybOmpV2VIVx553hachDDw3rI0lhtG9ffbkhUNIkIlKCVqwIyVKGO3zzTXrxFMqkSXDaafD++zBhAvzkJ2lHVDpOPDG8zKBzZ/jjH9OOqPaUNIlIvbNqVVjs8dVX046keG2yCfzsZ+XlgQNhzz1TC6dgPvqoYrlUetfqg8aN4YEHwv/f8+fDXnulHVHtKWkSkXpl1So44ICwOvbAgXDOOYVp98MP4Ykn4N//Lkx79cGf/hS2bXnyyfCzSQk8T33AARVX3j766NRCKVkNefNrTQSP0URwkfS9+GLYyDTDLMy/SXLxv4cfhpNOgrIy6NgxbKDao0dy7Um6Zs2CBx+EzTcPE+FLIVkUPT2Xd0qaRNI3ZUrFYaI2bcJk5UYJ9ovvvXfFxTIvvxx+85vk2hORgtPTc/liZqPM7EPg7bRjESl1ffqEpKVp0zCM8pe/JJswwfobpWrjVCk2S5bAvvuG/68GDqy4YbDkTj1NMeppEqncww/DP/8JgwfDgAGFaXPduuSTpYxZs+DII+Hjj8PmsX/7GzRvXpi2RQrhzDPDOkkZZ50FN9+cXjwp0DYqIpK8G2+Ec88Nx2PHwgsvhMm0SStUwgSw/fYhcSorC0/4iBSbr76qviy50fCciFTrscfKj8vK4Kmn0oslaYVMmNatg3/9C774onBt1gfr1sGiRbB2bdqRlJYRI8ofpmjZEs44I914GiolTSJSre23r74stVdWFoYDd94ZuneHW29NO6LCWLQIdt0VOnWCnj3DcKgUxn77hSH2hx6CDz4o3DB7salT0mRmI81stpmtNLMpZrZvDfVHmdkMM1thZrPMbHjW+03N7HIz+zS65jQzOySrTlszG2dmc6PrTDKzKpfGMrM7zMzN7Jy6fEcRCW68MWy3sdtucNFFYTVl2TAvvBAW74TQ83LuuWFF7mJ3/fUwfXo4njsXLr003XhKTc+eYQX0nj3TjqThqvWcJjM7HhgHjATeAE4HJpjZTu7+eSX1RwBjgf8C3gH6Anea2TfununovxoYGtWZCQwGHjOz/u7+flTnLmBnYBiwIKr/YtTu/Kw2jwL6RfVEZAO0awf33Zd2FMUle75Wqcyjim/bUllZpL6r9dNzZvYW8J67j4idmwE87u4XVVJ/EvCGu58XOzcO2NPdB0TlBcA17j4+VudxYJm7DzWzlsBSYIi7PxOrMxV42t0vjZ3rArxFSLyeAca5+7gcv5uenhORxK1bF/bgevjhsLDinXdW3NKkWH30URgmWrQoJOPPPx/WyBIpgMI/PWdmzYA+wLVZb00E+lfxseZA9r8nVgB9zaypu6+ppk5m1LUJ0LiGOphZI+B+4AZ3/5dZ9ffIzJpHbWe0rfYDIiJ50KhRmFtyww1h8c6NN047osLYbjuYMSNsWdOzZ1iRu5SsXavVxxu62s5p2pSQvCzKOr8I6FTFZ54HTjWzPhbsCZwCNI2ul6kz2sy2NbNGZjYIGAJsAeDuS4E3gcvMrLOZNTazoYQhuC1ibV0ArAX+kOP3uYjQs5R5zcvxcyJFyx0WL4Y1a9KOpPh161Y6CVNGhw6wzz6llTC9/TZ06RLW/ho+PPQ0SsNU16fnsoeirJJzGVcBE4DJwBrgCeDe6L2y6OfZwMeE+UyrgVuAe2LvQ5jLZMB8YBVwFvBApo6Z9Ymu8zPPfcxxLNA+9uqa4+dEitKyZWH4pGNH6No1bGmStNdeg223Db9Ex41Lvj2RQjv1VFiwICRL998fehmlYapt0rSEkKRk9yp1ZP3eJwDcfYW7nwK0AroD3YA5hDlKS6I6i939KKA1sBWwA7AMmB27zqfuPhBoA2zp7n0JvVWZOvtGcXxuZmvNbG10rd+b2ZwqYlvl7t9nXlFMIiXr9tvh9dfD8b//DaNHJ9/mscfCJ5+E9n71K5g6Nfk2RQrpu+8qlr/9NpUwJA9qlTS5+2pgCjAo661BwKQaPrvG3ee5exlwAmEC97qsOiujJ+GaAMcSeqWyr/ODuy80sw6Eyd6ZOvcDvYHdYq8FwA1RPRGpQaGfblq9OuyJFbdwYbJtihTa+edDZortNtvAccelG4/UXV2enjuekKCcQZhndBphqYBe7j7XzMYCXdx9eFR/O8IyA28BHYDRhCSrj7vPier0A7oAU6OfY4CtgT3c/duozmDC8NwsoCchGVoFDIgmk1cW6xz09JxIzr78Evr3h9mzw/yLRx4Je7El6ZRT4J57wvGOO8Jbb0FbPZIhReb992HevLBprjaETkU6e8+5+0NmtglwOWES9nTgUHefG1XZgjAEl9EYOBfYnjCn6WWgfyZhirQgrNXUgzAs9ywwLJMwRdoT5iB1Bb4GHgEuqSphEpHa69QJpk0Lr+7dw7ympN19NxxxRBjCOPpoJUxJuv9++PTTcL/79Ek7mtKy++7hJQ1brXuaipl6mqQ+KisLw1gtW6YdSfGYNw8mTQo9W7vsknY0hXHZZXD11eG4efPw/ffYI92YRAooLz1N2ntOpB57+unQld+qVdhwUzbczJnQuzccf3z4l/9f/5p2RIXx6KPlx6tWwbPPpheLSEOlpEmkHvvZz8IyABCebPv731MNpyjcey988004LiuDP+S6qlsDp42XRTac1iYVqafc4YcfKp7LJFBSd5tsUn25WN1xR9jj7pNP4D//s3Se4HKHl18Oq3EfeGDp7PMnydCcJsDMRgGjCD1v26M5TVJP/Pa3cMkl4bhfP3jlFWjRItWQGryVK8PQ3LPPQq9e8PjjYdK7FKehQ+EvfwnHhx4KTz21/obJUhLyMqdJSVOMJoJLfTR1Knz9ddh6onnzGqtLjtzL186R4jRvHmy5ZcVz06aFOW1SctJZckBECmu33dKOoDgpYSp+rVtD06bl+yiaQbt26cYkDZs6KUVEpCh16AB33hmePm3WDG66SUOxsmE0PBej4TkpNevWwTvvhF8qpbJekZQe9/DSXKaSpnWaRIrN5Mnwxz/CjBnJt7VuHQwZAnvvHeZ4ZCacixQbMyVMkh/qaYpRT5Ok6cEH4aSTwr+IW7QIj0nvvXdy7b3xBgwYUPHcsmVhHoiISJFRT5NIMbn77pAwQXgsPvOYdFKyly5o2hSa6NGQREycCFttBZttBrfemnY0IlJXSppE6onOnasv51ufPjB6dDhu2jQsfqglDfJv7dqwkOTnn8OSJXDmmYUZfhWR/NPwHFrcUuqHxYvD8NzUqTBoENxzT2GSmO++C08WFWpD4E8+gY8+gr32Cj0vxe7776F9+4rnXnkFBg5MJRyRUqXFLfNNc5pEkvXUU2ELj9WrYfPNYdIk6NEj7aiSd9JJ8L//G4533z1871JY2f2hh8Jw5Kabwrhx6y80KVJASpryTUmTSLIGDoRXXy0vX3ghjB2bXjyFsm4dPPFE2Evw6KNLY7L9++/DnnuG7w7h+J130o1JSppWBBeRhqVt24rlUlmduVGjkCyVkhkzyhMmgOnT04tFJF80EVykCvfdB127wjbbwHPPpR1NcbjxxnA/AX78YzjrrHTjkeTss0/FpPjQQ9OLRSRfNDwXo+E5yZgzB3r2hLKyUG7TBv7978JNli52K1eWxpyeUjd9evjHx2abhQRZT2dKijSnKd+UNEnG229Dv34Vz335ZZi8LCIiDY4WtxRJyq67homrGYcfroRJRKTUqacpRj1NEvfDD/DXv4YhheOO02rZIiINmIbn8kWLW4qIiBQ1JU35pp4mKTRNiE7e2rXw5pvhSa5dd007GhFJieY0iTRUq1bBEUeEp/G23BKmTUs7ouK0Zg0cfDDstx/sthuMGZN2RCLSkKmnKUY9TVJWBo0bJ9/OrbfCqFHl5QED4LXXkm+31LzwQkiaMho3hhUrwgbFxW7RIpgyBXbYoTS2qhGpgXqaRPLllVfC03EtWsDZZyff3tKlFcvf1zZFl5y0alWx3KJFYZLitM2cCb16wWGHhZ8vvJB2RCLFQUmTCDBsWFi8cu1a+MMfkv8lM2xY+ealjRuHPdgk//bZp3zV8RYt4E9/CluaFLvbb4evvgrHK1fC73+fbjwixaIE/voQqdk331RfzrfOncM8pgkTwqrJJ56YbHsAH3wAP/pRGK65667k26svbr459Ox99x385CdpR1MY2Xv6lcoefyJJ05ymGM1pKl1XXglXXBGOd94ZJk1af3PZhq5HD5g9OxybwXvvhcnRUny+/z4syPraa7D99iE533rrtKMSSZWWHMg3JU2lbfJkWLwYDjgg7DVXTNatC5Of47vOP/44DBmSWkhSACtWaL9EkYiSpnxT0iTF7KST4H//Nxx37QpTp8Imm6QakohIoShpyhetCC6lYO1auP/+MF/rhBPCvCoRkRKhpCnf1NMkIiJSlLROk4iIiEihKGmSeul3vwtPdg0ZAgsXph2NiIhIHZMmMxtpZrPNbKWZTTGzfWuoP8rMZpjZCjObZWbDs95vamaXm9mn0TWnmdkhWXXamtk4M5sbXWeSme2VVWeMmc00sx/M7Bsze9HM+tXlO0p6nn4azjsvrGP05JPw85+nHZGIiEgdkiYzOx4YB1wD7A68Bkwws25V1B8BjAXGAL2AK4DxZnZErNrVwOnAL4GdgNuBx8xs91idu4BBwDBgF2Ai8KKZdYnV+Qg4M3p/ADAHmGhmm9X2e0p6Pvqo+rI0HF99Bc89Bx9/nHYkIiIbrtYTwc3sLeA9dx8ROzcDeNzdL6qk/iTgDXc/L3ZuHLCnuw+IyguAa9x9fKzO48Aydx9qZi2BpcAQd38mVmcq8LS7X1pFrJmJ3Qe5+99z+G6aCF4PfPgh7LUXLF8eyqNHaxuIhmju3LAC+cKFYY2ohx6Co49OO6rkLVsG558Ps2aF73vmmWlHJCLkaSJ4k1q1aNYM6ANcm/XWRKB/FR9rDqzMOrcC6GtmTd19TTV1BsTibFxDncpiPY2QBE2rIjaph3baCd58Ex59FLp3h5NPTjsiqYs//al8PtqaNXDddaWRNJ15JvzP/4Tjl14KSzscc0y6MYlIftQqaQI2JSQvi7LOLwI6VfGZ54FTo56j9whJ1ylA0+h6C6M6o83sVeBT4EBgSNQW7r7UzN4ELot6tRYBJwL9gAod/2Z2OPAg0Cq69iB3X1JZYGbWnJCwZRTZxhkNV+/e4ZWkH36A1q2TbaOUZW9DUyr7n02dun5ZSZNIcajr03PZQ1FWybmMq4AJwGRgDfAEcG/0Xln082xC8jMTWA3cAtwTex/CXCYD5gOrgLOAB7LqALwM7Ebo+XoOeNjMOlYR20WEnqjMa14V9aSIzJ0bNq1t0wb69YOvv047ouI0ciQcdFA47tYNxo1LNZyCGTy4/LhRIzjwwPRiEZH8qtWcpmjIazlwnLs/Fjt/M7Cbuw+s5rNNgc0JvT+nAdcBG7n7ulidFsAmwALCEODh7t4r6zqtgXbuvtDMHgLauPth1bT7MfAndx9byXuV9TTNQ3OaitpPfwoPPFBePv/8MHQkySi1/c/WrYNbbw1zmo48EgYNSjsiESGNOU3uvtrMphCeYnss9tYgQg9SdZ9dQ9STY2YnECZwr8uqsxKYHyVYxwIPV3KdH4AfzKwDMBg4v4awjYqJUfxaqwi9VkRx1XApKQZLl1Ysf1/b9FhqpZQSJgi9S5r8LVKc6jI8dyNhjtIpZrajmd0EdCMsE4CZjTWz+zKVzWw7MxtqZtuaWV8zexDYGbg4VqefmR1jZj2iNZ+ei2K7PlZnsJkdYmZbm9kgwjDcLMIwHmbW2sx+a2Z7m9lWZraHmd0FdAX+WofvKQUyZw6cfTacey58+WXy7Z17LrRqFY433hhGjUq+zWeegb59Yb/94N13k29PpNAmTQp/xnfdNayvJlKM6rT3nJmNJPTwbAFMB37l7q9G790LdHf3/aPyjoS5R9sT5jS9DFzg7rNi1xsI3Ab0AJYBzwIXuvuCWJ2fENZ76gp8DTwCXOLu30Xvt4ja6UeYYP4V8A5wtbu/k+P30pIDBbZsWXha7osvQnmHHeCDD8Ij6kn6/HOYMSOsOr755sm2NW8e9OwJq6I+zc02g/nzk/+OIoWyalV4SjAzP7B5c/j0U+jSpfrPiRSQNuzNNyVNhffee9CnT8Vzs2eHpQaKxaRJsM8+Fc8tWgQdq3o8QaSBWbgwJE1xkyeHBy1E6glt2CsNX/fu0KFDeblzZ+hU1eIVDVTv3rDNNuXlAQNCb5NIsejUKQw9Z+ywA+yyS3rxiCRFPU0x6mlKx7vvwlVXQZMmcOWV0KtXzZ9paBYtgrvughYt4IwztD6UFJ/ly8OCpitXwimnhPmCIvWIhufyTUmTiIhIUdLwXL6Y2Sgz+xB4O+1YREREpH5ST1OMeppERESKknqaRETinnsuTEpu0wbGrrcHgIjIhlFPU4x6mqRYLF0a9n774IOwF9q114aVqovZunXhScz4Cu/vvQe7755eTCJSbxR+GxURaRh+/Wv485/D8QcfwFZbFWbl8zStXr3+ljhffZVOLCJSnIr8355SF9Onw333wcyZaUcidfXRR9WXi1GLFnD66eXlPfZYf1FREZENoZ4mqeCFF+Cww2DNmrAVwsSJFRetk4bhmGPglVfCcaNGcOSRqYZTMLffDkcfHYYnDz20dDYLXr06rMDdoYMWlRRJkuY0xWhOExx3HPztb+XlYcNCr5M0PI88EobmDjoI9t037WgkKStXwoEHhu16IMxfu+CCdGMSqYf09JzkX/Z+aElvZivJOfZY+M1vlDAVu4kTyxMmgDFjQP8WFkmGkia0uGXcVVeFf7W2bg2HHAKXXZZ2RCJSnVat1i9bXv5NLSLZlDQB7j7e3XcC+qYdS9o23hhefBGWLYMJE6Bdu/xef8qUsJlnhw5w8cX5vbZIKTroIDj11HDcqlXY/01EkqE5TTGa05S8HXes+FTe88/DwQenF49IsVi6NDxB2LRp2pGI1Eua0yQNz5dfVl9OwsUXQ/fucMAB8Pnnybcnkoa2bZUwiSRNSZMU1IgR5cdbbRUeC0/SX/8attOYOzc8gv+LXyTbnoiIFC+t0yQF9dvfwv77hx6mQw+FTTdNtr05cyqW585Ntj0RESle6mkqcXfdBQMGwIknwqJFhWnz4INh+PDkEyaAIUMqTmYfOjT5NkVEpDhpInhMqU0Ef/ll+PGPy8sHHxwmZhebjz+GZ58N85qGDEk7GhERSYEmgsuGmT69+nKx2HZbOPtsJUxJe/JJ+NGPwiPwxfpnSURKm+Y0lbADDwx7c61YEcqHHZZuPNJwffZZ2IJn9epQ/o//CE8qapFFESkmSpoIK4IDoyixnreddoJXX4WHH4auXWHkyLQjkobqs8/KEyaAefPCAqlt26YXk4hIvmlOU0ypzWmSZCxdGpY5+OqrsFLzXnulHVHyvvkGevcOyRKENbFeeindmArlzTfD/m+9e8PRR6cdjYhUIS/93kqaYpQ0ST4MHhx+iULYw++DD6BHj3RjKoQvvoC774Y2bUKvZfaeaMXo1VfDMPfataE8blyYPyci9Y6SpnxT0iT50KwZrFlTXn7wQTj++PTikeSMHg033VReHjAAXnstvXhEpEp6ek6kPurTp/y4WbMwbCPFadttK5a32y6dOESkMDQRXCTPHnss7He3ZEkYptpxx7QjkqScfjrMng0TJsAuu8CNN6YdkYgkScNzMRqeExERKUoanhMREREpFCVNIiIiIjlQ0kRY3NLMPgTeTjuWVavCXBgRERGpX5Q0Ae4+3t13AvqmGcff/w4dO8Jmm8GRR5av/SIiIiLpU9JUj4wYAd9H08+fegoeeijdeERERKSclhyoR1aurFjObKSbL2Vl8MILYRPVQYOgkVJmERGRnOnXZj1y1VXQuHE43nXX/K4i7Q7HHBN2nz/kkLAjfdKrTZSVwTXXhKHG3/8++fZERESSVKekycxGmtlsM1tpZlPMbN8a6o8ysxlmtsLMZpnZ8Kz3m5rZ5Wb2aXTNaWZ2SFadtmY2zszmRteZZGZ7ZV3jOjP7p5n9YGYLzOw+M+tcl++YhpNPho8+gtdfh8mT87tD/KxZ8OST5eVHHw2L8iXp2mvh0kvDUOOvfw233ppse6VOSamISLJqnTSZ2fHAOOAaYHfgNWCCmXWrov4IYCwwBugFXAGMN7MjYtWuBk4HfgnsBNwOPGZmu8fq3AUMAoYBuwATgRfNrEv0fitgD+Cq6OcxwHZALFWo/3r0gH32gRYt8nvddu3Ke7EAmjTJb1JWmbffrr4s+bFmTeg5bNoUdtgBZs5MOyIRkeJUl56m0cDd7n6Xu89w93OAL4ARVdQfBtzh7g+5+2fu/iBwN3BBVp3fuvuzUZ3bgOeBcwHMrCVwLHC+u7/q7p+4+xhgdqZdd//O3Qe5+8PuPsvdJxOSsD5VJXSlpHPn0NPTqhW0bg133BGe0kvS/vtXLA8cmGx7perOO+FvfwvDobNmwahRaUdUOG++CXffDZ99lnYkIlIKajUR3MyaAX2Aa7Pemgj0r+JjzYGsKc6sAPqaWVN3X1NNnQGxOBvXUKcy7Qnbm3xb2Ztm1jxqOyPhvpd0nXYa/Nd/hWPLy4Ly1TvnHGjePPQw7bcfnHJK8m2Wom++qVj++ut04ii0e+6BX/wiDEu2bQuTJsHOO6cdlYgUs9r2NG1KSF4WZZ1fBHSq4jPPA6eaWR8L9gROAZpG18vUGW1m25pZIzMbBAwBtgBw96XAm8BlZtbZzBqb2VCgX6ZONjNrQUjuHqhmH7mLCHvNZV7zqv/6DZ9ZYRKmTFsjR8K99yphStJPfwqbbx6OGzWCX/0q3XgK5Y9/LJ/HtXQpPPBAuvGISPGr69Nz2VNOrZJzGVcBE4DJwBrgCeDe6L2y6OfZwMfATGA1cAtwT+x9CEN4BswHVgFnAQ9k1QnBmDUFHiR8v5HVfI+xhN6ozKtrNXUTsXx5oVuUYtO9O0ybFobo3nsPhg+v8SNFoVOn6ssiIvlW26RpCSFJyf7rqSPr9z4B4O4r3P0UwkTt7kA3YA6wNLoe7r7Y3Y8CWgNbATsAywhzljLX+dTdBwJtgC3dvS+ht6rCM2BRwvQwsDUwqJpeJtx9lbt/n3lFMRXEokWwxx5hftFuu8GXXxaq5eK2ahX89a/w+ONhjk+p2HxzOPbYsFRFqfjv/4b+/WGjjWDYsLA4rIhIkmqVNLn7amAK4Sm2uEHApBo+u8bd57l7GXAC8LS7r8uqs9Ld5xPmMB1L6JXKvs4P7r7QzDoAg+N1YgnTtsBB7v5Vbb5fIV15Jbz/fjieNg3GjEk1nKKwdm1YtPMnP4Gjjw4/pXh17QpvvBHmdN13X3h6UEQkSXVZEfxG4H4ze5cwz+g0Qu/R7QBmNhbo4u7Do/J2hD3d3gI6EJ6+2xk4OXNBM+sHdAGmRj/HEBK662N1BhOG52YBPYEbouN7ovebAH8jLDdwONDYzDI9Yl9HCV+9sTSrT+v7KvvDJFcffACvvVZefvRRWLAgPDkoIiKyoWqdNLn7Q2a2CXA5YRL2dOBQd58bVdmCkERlNCYsHbA9YU7Ty0B/d58Tq9OCsFZTD8Kw3LPAMHf/NlanPWEOUlfga+AR4JLo6Tui80dGx1Ozwj4AeKW23zVJZ50VhpCWLg1P/pxzTtoRNXwbbxwmQq+L+i9btEh+LSoRESkd5lpG+P+ZWTvCU3Ttq5sLVYVa38j580PvyC67hKEG2XB33gkXXADNmsFtt4VhOhERKXl5eW5cSVNMoZMmERERKYi8JE3asFdEREQkB0qaRERERHKgpAkws1Fm9iGgLWVFRESkUprTFKM5TSIiIkVJc5pERERECkVJk4iIiEgOlDSJiIiI5EBJk8gGmj0bJkworU2Xv/gCTjoJDjsMXngh7WhERApDE8FjNBFcauull0LisHJl2Mbl1VehV6+0o0pe797wz3+G4+bNYfp06Nkz3ZhERKqhieAiafv970PCBPD113D77enGUwhr15YnTACrVsGHH6YXj4hIoShpKlJz58Khh8Juu8Ett6QdTfFq165iuRQ2CG7SBAYOLC9vtBH06ZNaOCIiBaPhOcLilsAoQhK5PQkMz/3973DvvdC5M1x2GbRpU7dYc7XPPjBpUnn5pZfggAOSbbMUzZkD//EfMHMm9O8PzzwTkohi9/33cP318O23cPrpYdNpEZF6TBv25ltSc5qmTYO99oI1a0L56KPh0UfrHGZOOnaExYvLy3fcAaedlmybpWzFCmjZMu0oRESkCprT1FBMnlyeMAG89lrybR5zTPlxu3Zw0EHJt1nKlDCJiBS/JmkHUAr22gsaN4ayslDee+/k27z11tDuggVw3HHQo0fybYqIiBQzDc/FJLnkwDPPlM9puvJKaN++rlGKiIhILWlOU75pnaYN8/jj8Oc/Q7du8JvflMaTZCIi0iDkJWnS8JzkxaRJcOyxsG5dKM+bBw8/nG5MIiIi+aSJ4JIXb79dnjBBmPwuIiJSTJQ0SV707x8mu2fsu296sYiIiCRBc5oozOKWpeC55+CBB8Kcposvhlat0o5IREQE0ETw/NNEcBERkaKkxS1FRERECkVJk4iIiEgOlDSJiIiI5EBJk0gDU1YW9i989920IxERKS1KmkQakLIyOPJI2G+/sLfgueemHZGISOnQ03MxenpO6rs33wxrYsV99x20a5dOPCIiDYSenhN57TV45JGQOJSC1q0rlps2DS8REUmekiZpsMaMCcNU//mf0K8ffPtt2hElr3fvsHAoQLNmcOed0LJlujGJiJQKDc+hFcEbqvbt4fvYf6W//AVOOim9eApp+XJo0iQkTiIiUiMNz+WLu493952AvmnHIrnbZJPqy8WsVSslTCIihaakSRqs+++Hrl2heXP41a9g8OC0IxIRkWKm4bkYPT0nIiJSlNIbnjOzkWY228xWmtkUM9u3hvqjzGyGma0ws1lmNjzr/aZmdrmZfRpdc5qZHZJVp62ZjTOzudF1JpnZXll1jjGz581siZm5me1Wl+8nIiIikq3WSZOZHQ+MA64BdgdeAyaYWbcq6o8AxgJjgF7AFcB4MzsiVu1q4HTgl8BOwO3AY2a2e6zOXcAgYBiwCzAReNHMusTqtAbeAC6s7fcSERERqU6th+fM7C3gPXcfETs3A3jc3S+qpP4k4A13Py92bhywp7sPiMoLgGvcfXyszuPAMncfamYtgaXAEHd/JlZnKvC0u1+a1WZ3YDawu7tPrcV30/CciIhI8Sn88JyZNQP6EHp54iYC/df/BADNgZVZ51YAfc2saQ11BkTHTYDGNdQRERERSUxth+c2JSQvi7LOLwI6VfGZ54FTzayPBXsCpwBNo+tl6ow2s23NrJGZDQKGAFsAuPtS4E3gMjPrbGaNzWwo0C9Tpy7MrLmZtcu8gLZ1vZaIiIgUt7ouOZA9FGWVnMu4CpgATAbWAE8A90bvlUU/zwY+BmYCq4FbgHti70OYy2TAfGAVcBbwQFad2rqIMByXec3bgGuJiIhIEatt0rSEkKRk9yp1ZP3eJwDcfYW7nwK0AroD3YA5hDlKS6I6i939KMJE7q2AHYBlhHlJmet86u4DgTbAlu7el9Bb9f916mAs0D726roB10rVunWweHH4KSIiIvlXq6TJ3VcDUwhPscUNAibV8Nk17j7P3cuAEwgTuNdl1Vnp7vMJc5iOJfRKZV/nB3dfaGYdgMGV1anF91nl7t9nXoRErsH54gvYcUfo2BF69YL589OOSEREpPg0qcNnbgTuN7N3CfOMTiP0Ht0OYGZjgS7uPjwqb0fYnuQtoAMwGtgZODlzQTPrB3QBpkY/xxASuutjdQYThudmAT2BG6Lje2J1No5i6Ryd2t7MAL509y/r8F0bhCuvhI8+CsczZ8JVV8Htt6cbk4iISLGpddLk7g+Z2SbA5YRJ2NOBQ919blRlC0LiktEYOJewEe4a4GWgv7vPidVpQVirqQdhWO5ZYJi7fxur054wnNYV+Bp4BLjE3dfE6hxJLIkCHox+/oaQiBWlFSuqL4uIiMiG0zYqMQ11nab334cDD4RvvoGNN4aXX4bevdOKRkREpN7JyzpNSppiGmrSBGES+IwZsNNOsOmmNdcXEREpIUqa8q0hJ01SGA88AA8/DD16hLljrVunHZGIiOQgL0lTXSaCi5SkF1+En/60vLx4Mdx/f3rxiIhIYdV1ccuiYmajzOxD4O20Y5H66513Kpbf1p8WEZGSoqQJcPfx7r4TYWkEkUrtuy80iv0fM3BgerGIiEjhaU5TjOY0SU2eeQb++lfYZhu44AJo1iztiEREJAeaCJ5vSppERESKUl6SJg3PSZ28/z4cfTQceyz8619pRyMiIpI89TTFqKcpN999F4anvvoqlDt1gs8+g5Yt041LRESkCuppknTMmVOeMAF8+SUsWJBaOCIiIgWhpElqrWdP6BbbXXDbbWHLLdOLR0REpBC0uKXUWuvW8OqrcOON4RH8887TU2QiIlL8NKeJsLglMIrQ87Y9mtMkIiJSTLTkQL5pIriIiEhR0kRwERERkUJR0iQiIiKSAyVNIiIiIjlQ0iQiIiKSAyVNIiIiIjlQ0iQiIiKSAyVNIiIiIjlQ0iQiIiKSAyVNhBXBzexD4O20YxEREZH6SSuCx2hFcBERkaKkFcGldLjDzTfD0KFw111pRyMiIqWoSdoBiOTi+uvhwgvD8V/+Ao0bw89/nm5MIiJSWtTTJA3Ca69VXxYREUmakiZpEPbeu2K5X7904hARkdKl4TlpEC66KAzJvfMO7L8/nH562hGJiEip0dNzMXp6TkREpCjp6TkRERGRQlHShBa3FBERkZppeC5Gw3MiIiJFScNzIiIiIoWipKkBuuIK6No1PIY/a1ba0YiIiJQGDc/FNIThuWeegcMPLy/vtRe8rZlYIiIi1UlveM7MRprZbDNbaWZTzGzfGuqPMrMZZrbCzGaZ2fCs95ua2eVm9ml0zWlmdkhWnbZmNs7M5kbXmWRme2XVMTMbY2YLojqvmFmvunzH+uqLL6ovi4iISDJqnTSZ2fHAOOAaYHfgNWCCmXWrov4IYCwwBugFXAGMN7MjYtWuBk4HfgnsBNwOPGZmu8fq3AUMAoYBuwATgRfNrEuszvnAaOBMYC/gS+AFM2tb2+9ZXx12GGy2WXn5Zz9LLRQREZGSUuvhOTN7C3jP3UfEzs0AHnf3iyqpPwl4w93Pi50bB+zp7gOi8gLgGncfH6vzOLDM3YeaWUtgKTDE3Z+J1ZkKPO3ul5qZAQuAce5+XfR+c2ARcIG735HDd6v3w3MAn38OTz8NXbrAkCGFalVERKTBysvwXK22UTGzZkAf4NqstyYC/av4WHNgZda5FUBfM2vq7muqqTMgFmfjGupsDXSKYgHA3VeZ2T+i2NZLmqKkqnnsVIPokerWDUaOTDsKERGR0lLb4blNCcnLoqzziwgJS2WeB041sz7RnKM9gVOAptH1MnVGm9m2ZtbIzAYBQ4AtANx9KfAmcJmZdTazxmY2FOiXqRNrvzaxXUToWcq85lX91UVERKSU1XXJgeyhKKvkXMZVwARgMrAGeAK4N3qvLPp5NvAxMBNYDdwC3BN7H8JcJgPmA6uAs4AHsurUNraxQPvYq2sV9URERKTE1TZpWkJIUrJ7bjqyfg8PAO6+wt1PAVoB3YFuwBzCHKUlUZ3F7n4U0BrYCtgBWAbMjl3nU3cfCLQBtnT3voTeqkydL6OftYltlbt/n3lFMYmIiIisp1ZJk7uvBqYQnmKLGwRMquGza9x9nruXAScQJnCvy6qz0t3nE+YwHUvolcq+zg/uvtDMOgCDY3VmExKn/48tmoM1sKbYRERERGpSq4ngkRuB+83sXcI8o9MIvUe3A5jZWKCLuw+PytsBfYG3gA6EJQF2Bk7OXNDM+gFdgKnRzzGEhO76WJ3BhKG2WUBP4Ibo+B4Ad/foqbyLzexjwnDfxcBywjCeiIiISJ3VOmly94fMbBPgcsIk7OnAoe4+N6qyBSGJymgMnAtsT5jT9DLQ393nxOq0IKzV1IMwLPcsMMzdv43VaU+Yg9QV+Bp4BLgkevou43qgJXArIUF7Czg4mkguIiIiUmfaRiWmoazTJCIiIrWS3jYqIiIiIqVGSZOIiIhIDpQ0iYiIiORASRNgZqPM7EPg7bRjqY+WLYOLLw6bA//jH2lHIyIikg5NBI/RRPDKHXUUPBGthtWsGUyZAjvvnGpIIiIitaGJ4FIYr75afrx6Nbz1VnqxiIiIpEVJk9Sob9/y4yZNYI890otFREQkLRqei9HwXOW++QYuvRQWLoRf/AIOOyztiERERGolL8NzSppilDSJiIgUJc1pEhERESkUJU0iIiIiOVDSJCIiIpIDJU1ocUsRERGpmSaCx2giuIiISFHSRHARERGRQlHSJCIiIpIDJU0iIiIiOVDSJCIiIpIDJU0iIiIiOVDSJCIiIpIDJU0iIiIiOVDShBa3FBERkZppccsYLW4pIiJSlLS4pYiIiEihKGkSERERyYGSJhEREZEcKGkSERERyYGSJhEREZEcKGkSERERyYGSJhEREZEcKGlCi1uKiIhIzbS4ZYwWtxQRESlKWtyyVKxcGV4iIiKSHiVN9dxNN0Hr1uF1441pRyMiIlK6NDwXU9+G5xYsgK5dIfOfyAzmzYPOnfPdkoiISFFLb3jOzEaa2WwzW2lmU8xs3xrqjzKzGWa2wsxmmdnwrPebmtnlZvZpdM1pZnZIVp0mZnZ11O4KM/ss+kyjWJ3NzexeM1tgZsvN7Dkz27Yu37E+WLGiPGGCcLxiRXrxiIiIlLJaJ01mdjwwDrgG2B14DZhgZt2qqD8CGAuMAXoBVwDjzeyIWLWrgdOBXwI7AbcDj5nZ7rE6FwBnAGcCOwLnA+dFn8HMDHgc6AEMiWKbC7xoZq1r+z3rg222geGx9HLYsHBORERECq/Ww3Nm9hbwnruPiJ2bATzu7hdVUn8S8Ia7nxc7Nw7Y090HROUFwDXuPj5W53FgmbsPjcpPA4vc/RexOo8Ay919mJltB8wCdnb3f0XvNwb+DVzg7nfl8N3q1fBcxuTJ4efeeyfVgoiISFEr/PCcmTUD+gATs96aCPSv4mPNgexnv1YAfc2saQ11BsTKrwMHRskRZrZr9P6zsWsQv467lwGrs67T4Oy9txImERGRtNV2eG5ToDGwKOv8IqBTFZ95HjjVzPpYsCdwCtA0ul6mzmgz29bMGpnZIMIQ2xax61wH/C8w08zWAO8D49z9f6P3ZxKG48aaWQcza2ZmF0Zxxa/z/8ysuZm1y7yAtrneCBERESktdV1yIHsoyio5l3EVMAGYDKwBngDujd4ri36eDXxMSHxWA7cA98TeBzgeGAqcBOwBnAz82sxOBnD3NcCxwHbA18ByYP+o7fh14i4iDMdlXvOq+sIiIiJS2mo1pykanlsOHOfuj8XO3wzs5u4Dq/lsU2BzYCFwGqHnaCN3Xxer0wLYBFgAXAsc7u69ove+AK7Nmvd0KTDU3XfIaqs90MzdF0dzsN5191GVxNSc8mE9CD1N86jbnCYREREpYrXqaXL31cAUYFDWW4OASTV8do27z4vmGZ0APB1PmKI6K919PtCE0Gv0ROztVkCF+oQepPW+g7t/FyVM2wJ7Zl0nXm+Vu3+feRGStfbA0uq+i4iIiJSeJnX4zI3A/Wb2LvAmodeoG2GZAMxsLNDF3YdH5e2AvsBbQAdgNLAzYXiNqE4/oAswNfo5hpAMXR9r9yngEjP7HPgXYUmB0cCfYtc5DlgMfA7sAtxMeKove+J6pTx0u6mHSURERNZT66TJ3R8ys02AywkTrKcDh7r73KjKFoQkKqMxcC6wPWFO08tAf3efE6vTgrBWUw9gGeGJuGHu/m2szi8J86NuBToSeoXuAK6M1dmCkNRlhgHviz4jIiIiskG0jYqIiIhIDrRhr4iIiEgOlDSJiIiI5EBJk4iIiEgOlDSJiIiI5EBJk4iIiEgOlDSJiIiI5EBJk4iIiEgOlDSJiIiI5KAu26hIFjMzwma/IiIiUj8t9Q1c0VtJU360Bb5LOwgRERGpUns2cH9ZbaOSB7XoaWoLzAO6AksTCOVtwubISWnI19e9T+/6uvfpXV/3Pr3r696nc/3q7rt6muqD6D9CjdlryK2A8B9ug7LdKq6/LonrFsP1de/Tu77ufXrX171P7/q69+lcP+n7rongxWW8rp+ahn5vdO+L9/pJauj3Rve+eK+fCA3PFZCZtSPMfWqfZAYv69O9T4/ufXp079Oje5+OpO+7epoKaxXwm+inFJbufXp079Oje58e3ft0JHrf1dMkIiIikgP1NImIiIjkQEmTiIiISA6UNImIiIjkQEmTiIiISA6UNG0gM9vPzJ4yswVm5mZ2VNb7ZmZjovdXmNkrZtYrq05zM/tvM1tiZj+Y2ZNm1rWgX6QBqu7em1lTM7vOzP4Z3dMFZnafmXXOuobufR3U9Oc+q+4dUZ1zss7r3tdBLvfezHaM7ud3ZrbUzCabWbfY+7r3dZDD3/dtzOwWM5sX/X0/w8xGZNXRva8lM7vIzN6J/iz/28weN7Pts+oU5HetkqYN1xqYBpxZxfvnA6Oj9/cCvgReMLP4tivjgKOBE4ABQBvgaTNrnFDMxaK6e98K2AO4Kvp5DLAd8GRWvXHo3tdFTX/uAYh+qfQDFlTy9jh07+ui2ntvZtsArwMzgf2BXQn/H6yMVRuH7n1d1PTn/ibgEGAosGNU/m8zGxKrMw7d+9oaSFgMc29gEGE3k4lm1jpWpzC/a91drzy9AAeOipUNWAhcEDvXHPgWOD0qtwdWA8fH6nQGyoDBaX+nhvLKvvdV1NkrqtdN9z75ew90IewB1QuYA5wTe0/3PqF7DzwI3F/NZ3Tvk7v304HLss5NAa7Svc/rvd8suv/7ReWC/a5VT1OytgY6ARMzJ9x9FfAPoH90qg/QNKvOAsL/fP2RfGpP+B/t26ise58QM2sE3A/c4O7/qqSK7n0Covt+GPCRmT0fDWW8lTWMpHufnNeBI82sSzRcdAChh/v56H3d+/xoH/38OvpZsN+1SpqS1Sn6uSjr/KLYe52A1e7+TTV1ZAOZWQvgWuABL19aX/c+ORcAa4E/VPG+7n0yOhKGHC4EngMOBh4DHjWzgVEd3fvknAV8SOhhXU34bzDS3V+P3te930BmZsCNwOvuPj06XbDftU1qF67UUfay61bJuWy51JEcmFlTwpBFI2BkLh9B977OzKwPcDawh0d94LX5OLr3GyLzD+En3P2m6HiqmfUHziD8y7squvcb7izCvJsjgbnAfsCtZrbQ3V+s5nO697m7BehNmJOULfHfteppStaX0c/sLLYj5Rnxl0AzM+tQTR2poyhhepjQfTvIK27gqHufjH0J9/BzM1trZmuBrYDfm9mcqI7ufTKWEHr4Psw6PwPIPD2ne58AM2sJ/BYY7e5PufsH7n4L8BDw66ia7v0GMLP/JiSkB7j7vNhbBftdq6QpWbMJ/6EGZU6YWTPCkwCTolNTgDVZdbYAdo7VkTqIJUzbAge5+1dZVXTvk3E/4V+Cu8VeC4AbgMFRHd37BLj7auAdYPust7Yj9HyA7n1SmkavdVnnyyj/Xat7XwfR/LBbCE9B/9jdZ2dVKdjvWg3PbSAzawP0jJ3a2sx2A75298/NbBxwsZl9DHwMXAwsBx4AcPfvzOxuwr/CvyJMbPsd8E+guu7cklfdvSf8kv4bYbmBw4HGZpb5V8jX7r5a977uavpzD3yVVX8N8KW7zwL9ud8QOdz7G4CHzOxV4GXCI/BHEJYf0L3fADn8ff8P4AYzW0FIUgcCwwmPwuve19144CRgCLA09nf5d+6+wt29YL9r0350sKG/CH8ReSWve738UcgxhMchVxLmFOycdY0WwH8TftEsB54Ctkz7u9X3V3X3HuhexXsO7K97n9y9r6L+HGJLDujeJ3vvgVMIvzhWAFOBIbr3yd97wvDQPcD86N7PJCRMpnu/Qfe9qr/LfxarU5DftRZdSERERESqoTlNIiIiIjlQ0iQiIiKSAyVNIiIiIjlQ0iQiIiKSAyVNIiIiIjlQ0iQiIiKSAyVNIiIiIjlQ0iQiIiKSAyVNIiIiIjlQ0iQiIiKSAyVNIiIiIjlQ0iQiIiKSg/8D7JuLg5JPQOsAAAAASUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 100 graphics primitives" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Plot ratio of series coefficients and asymptotic behaviour\n", "PTS = [binomial(round(pi*k)+k,k)/ASM.subs(n=k).n()/M(k) for k in range(100,201)]\n", "\n", "pt = plot([])\n", "for k in range(100):\n", " pt += point([k+100,PTS[k]])\n", " \n", "pt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 5.4 (Vanishing of Asymptotic Terms)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The main diagonal of F(x,y) = -1/(x + y - 1) has expansion beginning\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{4^{n}}{\\sqrt{n}} \\Bigg( \\frac{1}{\\sqrt{\\pi}} - \\frac{1}{8 \\, \\sqrt{\\pi} n} + \\frac{1}{128 \\, \\sqrt{\\pi} n^{2}} + O\\Bigg( \\frac{1}{n^{3}} \\Bigg)\\Bigg)\n", "\\end{math}" ], "text/plain": [ "4^n/sqrt(n) \\Bigg( 1/sqrt(pi) - 1/8/(sqrt(pi)*n) + 1/128/(sqrt(pi)*n^2) + O\\Bigg( n^(-3) \\Bigg)\\Bigg)" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "The main diagonal of F(x,y) = (2*y^2 - x)/(x + y - 1) has expansion beginning\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{4^{n}}{\\sqrt{n}} \\Bigg( \\frac{1}{4 \\, \\sqrt{\\pi} n} + \\frac{3}{32 \\, \\sqrt{\\pi} n^{2}} + O\\Bigg( \\frac{1}{n^{3}} \\Bigg)\\Bigg)\n", "\\end{math}" ], "text/plain": [ "4^n/sqrt(n) \\Bigg( 1/4/(sqrt(pi)*n) + 3/32/(sqrt(pi)*n^2) + O\\Bigg( n^(-3) \\Bigg)\\Bigg)" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "For the main diagonal of F(x,y) = -(x - y)/(x + y - 1) we conclude only\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{4^{n}}{\\sqrt{n}} \\Bigg( 0 + O\\Bigg( \\frac{1}{n^{3}} \\Bigg)\\Bigg)\n", "\\end{math}" ], "text/plain": [ "4^n/sqrt(n) \\Bigg( 0 + O\\Bigg( n^(-3) \\Bigg)\\Bigg)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Define parameters\n", "var('x,y')\n", "Ga = 1; Gb = x-2*y^2; Gc = x-y\n", "H = 1-x-y\n", "vars = [x,y]\n", "R = [1,1]\n", "M = 3\n", "\n", "# Get critical point as function of r and s\n", "CPs = critpt(H,R,vars)\n", "CP = CPs[0]\n", "\n", "# Parametrize y = g(x) = 1-x on V(H)\n", "g = 1-x\n", "\n", "# Print asymptotics\n", "print(\"The main diagonal of F(x,y) = {} has expansion beginning\".format(Ga/H))\n", "ex,pw,se = smoothContrib(Ga,H,R,vars,CP,M,g)\n", "disp_asm(ex,pw,se.expand(),M)\n", "\n", "print(\"The main diagonal of F(x,y) = {} has expansion beginning\".format(Gb/H))\n", "ex,pw,se = smoothContrib(Gb,H,R,vars,CP,M,g)\n", "disp_asm(ex,pw,se.expand(),M)\n", "\n", "print(\"For the main diagonal of F(x,y) = {} we conclude only\".format(Gc/H))\n", "ex,pw,se = smoothContrib(Gc,H,R,vars,CP,M,g)\n", "disp_asm(ex,pw,se,M)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 5.5 (Asymptotics of a Multinomial Expansion)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The [a, b, c] diagonal of F(x,y) = -1/(x + y + z - 1) has expansion beginning\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\left(\\frac{{\\left(a + b + c\\right)}^{a + b + c}}{a^{a} b^{b} c^{c}}\\right)^{n}}{c n} \\Bigg( \\frac{\\sqrt{a + b + c} \\sqrt{c}}{2 \\, \\pi \\sqrt{a} \\sqrt{b}} - \\frac{{\\left(a + b\\right)} {\\left(a + c\\right)} {\\left(b + c\\right)}}{24 \\, \\pi \\sqrt{a + b + c} a^{\\frac{3}{2}} b^{\\frac{3}{2}} \\sqrt{c} n} + O\\Bigg( \\frac{1}{n^{2}} \\Bigg)\\Bigg)\n", "\\end{math}" ], "text/plain": [ "((a + b + c)^(a + b + c)/(a^a*b^b*c^c))^n/(c*n) \\Bigg( 1/2*sqrt(a + b + c)*sqrt(c)/(pi*sqrt(a)*sqrt(b)) - 1/24*(a + b)*(a + c)*(b + c)/(pi*sqrt(a + b + c)*a^(3/2)*b^(3/2)*sqrt(c)*n) + O\\Bigg( n^(-2) \\Bigg)\\Bigg)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Set basic info\n", "var('x,y,z,a,b,c')\n", "G = 1; H = 1-x-y-z\n", "r = [a,b,c]\n", "vars = [x,y,z]\n", "CP = critpt(H,r,vars)[0]\n", "M = 2\n", "g = solve(H,vars[-1])[0].rhs()\n", "\n", "# Determine asymptotics\n", "ex,pw,se = smoothContrib(G,H,r,vars,CP,M,g)\n", "\n", "# Simplify and print asymptotic expression\n", "var('N')\n", "se_simp = add([k.factor()/n^p for [k,p] in se.subs(n=1/N).series(N,2).coefficients()])\n", "print(\"The {} diagonal of F(x,y) = {} has expansion beginning\".format(r,G/H))\n", "disp_asm(ex,pw,se_simp,M)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 5.6 (Lattice Path Asymptotics)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Set basic info\n", "var('x,y,t')\n", "G = (1+x)*(1+y); H = 1-t*x*y*(x+y+1/x+1/y)\n", "r = [1,1,1]\n", "vars = [x,y,t]\n", "sigma = critpt(H,r,vars)[0]\n", "tau = critpt(H,r,vars)[1]\n", "M = 4\n", "g = solve(H,vars[-1])[0].rhs()\n", "\n", "# Determine asymptotics\n", "ex1,pw1,se1 = smoothContrib(G,H,r,vars,sigma,M,g)\n", "ex2,pw2,se2 = smoothContrib(G,H,r,vars,tau,M,g)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The asymptotic contribution of {x: 1, y: 1, t: 1/4} to diagonal asymptotics is\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{4^{n}}{n} \\Bigg( \\frac{4}{\\pi} - \\frac{6}{\\pi n} + \\frac{19}{2 \\, \\pi n^{2}} - \\frac{63}{4 \\, \\pi n^{3}} + O\\Bigg( \\frac{1}{n^{4}} \\Bigg)\\Bigg)\n", "\\end{math}" ], "text/plain": [ "4^n/n \\Bigg( 4/pi - 6/(pi*n) + 19/2/(pi*n^2) - 63/4/(pi*n^3) + O\\Bigg( n^(-4) \\Bigg)\\Bigg)" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "The asymptotic contribution of {x: -1, y: -1, t: -1/4} to diagonal asymptotics is\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\left(-4\\right)^{n}}{n} \\Bigg( \\frac{1}{\\pi n^{2}} - \\frac{9}{2 \\, \\pi n^{3}} + O\\Bigg( \\frac{1}{n^{4}} \\Bigg)\\Bigg)\n", "\\end{math}" ], "text/plain": [ "(-4)^n/n \\Bigg( 1/(pi*n^2) - 9/2/(pi*n^3) + O\\Bigg( n^(-4) \\Bigg)\\Bigg)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "print(\"The asymptotic contribution of {} to diagonal asymptotics is\".format(sigma))\n", "disp_asm(ex1,pw1,se1.expand(),M)\n", "print(\"The asymptotic contribution of {} to diagonal asymptotics is\".format(tau))\n", "disp_asm(ex2,pw2,se2.expand(),M)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[(4)^{n}n^{-1}\\Bigl(1 - \\frac{3}{2}n^{-1} + \\frac{19}{8}n^{-2} - \\frac{63}{16}n^{-3} + \\frac{871}{128}n^{-4} + O(n^{-5})\\Bigr), (-4)^{n}n^{-3}\\Bigl(1 - \\frac{9}{2}n^{-1} + \\frac{107}{8}n^{-2} - \\frac{525}{16}n^{-3} + \\frac{9263}{128}n^{-4} + O(n^{-5})\\Bigr)\\right]\n", "\\end{math}" ], "text/plain": [ "[4^n*n^(-1)*(1 - 3/2*n^(-1) + 19/8*n^(-2) - 63/16*n^(-3) + 871/128*n^(-4) + O(n^(-5))),\n", " (-4)^n*n^(-3)*(1 - 9/2*n^(-1) + 107/8*n^(-2) - 525/16*n^(-3) + 9263/128*n^(-4) + O(n^(-5)))]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Note that additional terms in the expansion can be computed very efficiently \n", "# from P-recurrence satisfied by the diagonal sequence / D-finite equation of the GF\n", "from ore_algebra import *\n", "Pols. = PolynomialRing(QQ)\n", "Shift. = OreAlgebra(Pols)\n", "\n", "# Import recurrence satisfied by diagonal sequence (obtained from kernel method and creative telescoping)\n", "rec = (-n^2 - 7*n - 12)*Sn^2 + (8*n + 20)*Sn + 16*n^2 + 48*n + 32\n", "\n", "# The asymptotic expansion of the diagonal will be a C-linear combination of these terms\n", "# Looking at the dominant terms allows one to compute the coefficients (here 4/pi and 1/pi)\n", "# This approach can compute about 20 terms in the expansion in one second on a modern laptop\n", "show(rec.generalized_series_solutions(5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 5.7 (Lonesum Matrices)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Set basic info\n", "var('x,y,a,b,r,s')\n", "G = exp(x+y); H = exp(x)+exp(y)-exp(x+y)\n", "R = [r,s]\n", "vars = [x,y]\n", "CP = dict([(x,a),(y,b)])\n", "M = 1\n", "g = x - log(exp(x)-1)\n", "\n", "# Get asymptotics\n", "ex,pw,se = smoothContrib(G,H,R,vars,CP,M,g)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1/2*sqrt(2)*sqrt(pi)*sqrt(-((2*a - 1)*e^a + 1)*b*r*s*e^b + (b^2*(e^a - 1) + b*(e^a - 1))*r^2*e^b + (a^2*e^(a + b) - a^2*e^a)*s^2)*sqrt(b)*s*sqrt(e^a - 1)*e^(a + 1/2*b)/((pi*a - (2*pi*a^2 - pi*a)*e^(2*a) + 2*(pi*a^2 - pi*a)*e^a - (pi + (pi - 2*pi*a)*e^(2*a) - 2*(pi - pi*a)*e^a)*log(e^a - 1))*b*r*s*e^b + ((pi*a*e^(2*a) - 2*pi*a*e^a + pi*a - (pi + pi*e^(2*a) - 2*pi*e^a)*log(e^a - 1))*b^2 + (pi*a*e^(2*a) - 2*pi*a*e^a + pi*a - (pi + pi*e^(2*a) - 2*pi*e^a)*log(e^a - 1))*b)*r^2*e^b - (pi*a^3*e^(2*a) - pi*a^3*e^a - (pi*a^3*e^(2*a) - pi*a^3*e^a - (pi*a^2*e^(2*a) - pi*a^2*e^a)*log(e^a - 1))*e^b - (pi*a^2*e^(2*a) - pi*a^2*e^a)*log(e^a - 1))*s^2)\n" ] } ], "source": [ "# The leading constant is very complicated, as Sage doesn't know the relationship\n", "# between (x,y) and (a,b).\n", "print(se)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\log\\left(2\\right)^{-2 \\, n}}{\\sqrt{n}} \\Bigg( -\\frac{i \\, \\sqrt{\\pi} \\sqrt{\\log\\left(2\\right) - 1}}{\\pi \\log\\left(2\\right)^{2} - \\pi \\log\\left(2\\right)} + O\\Bigg( \\frac{1}{n} \\Bigg)\\Bigg)\n", "\\end{math}" ], "text/plain": [ "log(2)^(-2*n)/sqrt(n) \\Bigg( -I*sqrt(pi)*sqrt(log(2) - 1)/(pi*log(2)^2 - pi*log(2)) + O\\Bigg( 1/n \\Bigg)\\Bigg)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# It can be hard to get Sage to fully simplify in the general case,\n", "# but easy to simplify for explicit points. Here is one example.\n", "sbs = [r==1,s==1,a==log(2),b==log(2)]\n", "disp_asm(ex.subs(sbs), pw.subs(sbs), se.subs(sbs).canonicalize_radical(),M)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|C|\\phantom{\\verb!x!}\\verb|=| \\frac{1}{\\sqrt{-\\pi {\\left(\\log\\left(2\\right) - 1\\right)}} \\log\\left(2\\right)}\n", "\\end{math}" ], "text/plain": [ "'C = ' 1/(sqrt(-pi*(log(2) - 1))*log(2))" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# The leading constant has been published in the literature as\n", "C = 1/sqrt(pi*(1-log(2)))/log(2)\n", "show(\"C = \",C)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# This is the same as our computed constant\n", "(C/se.subs(sbs)).canonicalize_radical()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 5.9 (Asymptotics of Apéry Numbers)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The Apéry numbers for ζ(3) have asymptotic expansion\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{{\\left(12 \\, \\sqrt{2} + 17\\right)}^{n}}{n^{\\frac{3}{2}}} \\Bigg( \\frac{\\sqrt{2} \\sqrt{17 \\, \\sqrt{2} + 24}}{8 \\, \\pi^{\\frac{3}{2}}} + O\\Bigg( \\frac{1}{n} \\Bigg)\\Bigg)\n", "\\end{math}" ], "text/plain": [ "(12*sqrt(2) + 17)^n/n^(3/2) \\Bigg( 1/8*sqrt(2)*sqrt(17*sqrt(2) + 24)/pi^(3/2) + O\\Bigg( 1/n \\Bigg)\\Bigg)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Apéry sequence for ζ(3)\n", "var('x,y,z,t')\n", "G = 1; H = 1 - t*(1+x)*(1+y)*(1+z)*(1+y+z+y*z+x*y*z)\n", "r = [1,1,1,1]\n", "vars = [x,y,z,t]\n", "CPs = critpt(H,r,vars)\n", "CP = CPs[1]\n", "M = 1\n", "g = solve(H,vars[-1])[0].rhs()\n", "field = QQbar\n", "\n", "# Determine asymptotics\n", "ex,pw,se = smoothContrib(G,H,r,vars,CP,M,g)\n", "\n", "# Print results\n", "print(\"The Apéry numbers for ζ(3) have asymptotic expansion\")\n", "disp_asm(ex,pw,se.expand(),M)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The Apéry numbers for ζ(2) have asymptotic expansion\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{{\\left(\\frac{5}{2} \\, \\sqrt{5} + \\frac{11}{2}\\right)}^{n}}{n} \\Bigg( -\\frac{\\sqrt{2} \\sqrt{-11 \\, \\sqrt{5} + 25}}{2 \\, \\pi {\\left(11 \\, \\sqrt{5} - 25\\right)}} + O\\Bigg( \\frac{1}{n} \\Bigg)\\Bigg)\n", "\\end{math}" ], "text/plain": [ "(5/2*sqrt(5) + 11/2)^n/n \\Bigg( -1/2*sqrt(2)*sqrt(-11*sqrt(5) + 25)/(pi*(11*sqrt(5) - 25)) + O\\Bigg( 1/n \\Bigg)\\Bigg)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Apéry sequence for ζ(2)\n", "var('x,y,z')\n", "G = 1; H = 1 -z*(1+x)*(1+y)*(1+y+x*y)\n", "r = [1,1,1]\n", "vars = [x,y,z]\n", "CPs = critpt(H,r,vars)\n", "CP = CPs[0]\n", "M = 1\n", "g = solve(H,vars[-1])[0].rhs()\n", "\n", "# Determine asymptotics\n", "ex,pw,se = smoothContrib(G,H,r,vars,CP,M,g)\n", "\n", "# Print results\n", "print(\"The Apéry numbers for ζ(2) have asymptotic expansion\")\n", "disp_asm(ex,pw,se,M)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1/20*sqrt(110*sqrt(5) + 250)/pi" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# In the textbook example we use the following value of the constant\n", "C = (sqrt(250+110*sqrt(5))/20/pi)\n", "C" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[1]" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# This is equivalent to what we calculate here\n", "(se/C).minpoly().roots(multiplicities=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 5.10 (Asymptotics of Bar Graphs)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}[x^iy^j]B(x,y) = [x^iy^jz^j] \\frac{{\\left(x y z + 2 \\, x z + y z + x - 1\\right)} z}{x y z + x y + x z + y z + x - 1}\n", "\\end{math}" ], "text/plain": [ "[x^iy^j]B(x,y) = [x^iy^jz^j] (x*y*z + 2*x*z + y*z + x - 1)*z/(x*y*z + x*y + x*z + y*z + x - 1)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Recompute the rational function encoding the generating function from Example 3.12\n", "var('x,y,z,a')\n", "P = z - x*y - (x+y+x*y)*z - x*z^2\n", "R = (z^2*(diff(P,z)/P).subs(y=y*z)).factor()\n", "show(LatexExpr(\"[x^iy^j]B(x,y) = [x^iy^jz^j]\"),R)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left\\{x : \\frac{\\sqrt{a^{2} + 4} - 2}{a}, y : -\\frac{1}{2} \\, a + \\frac{1}{2} \\, \\sqrt{a^{2} + 4}, z : -\\frac{1}{2} \\, a + \\frac{1}{2} \\, \\sqrt{a^{2} + 4}\\right\\}, \\left\\{x : -\\frac{\\sqrt{a^{2} + 4} + 2}{a}, y : -\\frac{1}{2} \\, a - \\frac{1}{2} \\, \\sqrt{a^{2} + 4}, z : -\\frac{1}{2} \\, a - \\frac{1}{2} \\, \\sqrt{a^{2} + 4}\\right\\}\\right]\n", "\\end{math}" ], "text/plain": [ "[{x: (sqrt(a^2 + 4) - 2)/a,\n", " y: -1/2*a + 1/2*sqrt(a^2 + 4),\n", " z: -1/2*a + 1/2*sqrt(a^2 + 4)},\n", " {x: -(sqrt(a^2 + 4) + 2)/a,\n", " y: -1/2*a - 1/2*sqrt(a^2 + 4),\n", " z: -1/2*a - 1/2*sqrt(a^2 + 4)}]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# There are two critical points, one of which is has positive coordinates\n", "G = R.numerator(); H = R.denominator()\n", "r = [a,1,1]\n", "vars = [x,y,z]\n", "CPs = critpt(H,r,vars)\n", "show(CPs)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{x: (sqrt(a^2 + 4) - 2)/a,\n", " y: -1/2*a + 1/2*sqrt(a^2 + 4),\n", " z: -1/2*a + 1/2*sqrt(a^2 + 4)}" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Select the critical point with positive coordinates\n", "[CP] = [k for k in CPs if x(k).subs(a=1)>0]\n", "CP" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}-\\frac{{\\left(a^{3} t^{2} - \\sqrt{a^{2} + 4} a^{2} t^{2} + 2 \\, a^{2} t^{2} - 2 \\, \\sqrt{a^{2} + 4} a t^{2} + 4 \\, a t^{2} - 2 \\, \\sqrt{a^{2} + 4} t^{2} - 2 \\, a t + 4 \\, t^{2} + 2 \\, \\sqrt{a^{2} + 4} t - 2 \\, a - 4 \\, t\\right)} {\\left(t - 1\\right)}}{2 \\, a}\n", "\\end{math}" ], "text/plain": [ "-1/2*(a^3*t^2 - sqrt(a^2 + 4)*a^2*t^2 + 2*a^2*t^2 - 2*sqrt(a^2 + 4)*a*t^2 + 4*a*t^2 - 2*sqrt(a^2 + 4)*t^2 - 2*a*t + 4*t^2 + 2*sqrt(a^2 + 4)*t - 2*a - 4*t)*(t - 1)/a" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# This point is minimal as long as the following polynomial in t has no root in (0,1)\n", "var('t')\n", "p = H.subs([v == v*t for v in vars]).subs(CP).factor()\n", "show(p)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-a^3*t^2 + sqrt(a^2 + 4)*a^2*t^2 - 2*a^2*t^2 + 2*sqrt(a^2 + 4)*a*t^2 - 4*a*t^2 + 2*sqrt(a^2 + 4)*t^2 + 2*a*t - 4*t^2 - 2*sqrt(a^2 + 4)*t + 2*a + 4*t" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Thus we need to prove the following quadratic has, for any a>0, no root of t in (0,1)\n", "q = p/(t-1)*2*a\n", "q" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The value of q(t) at the vertex t = -1/2*(a^2 + sqrt(a^2 + 4)*(a + 1) + a + 2)/a is a - sqrt(a^2 + 4)\n" ] } ], "source": [ "# The vertex of the parabola in t defined by q(t)=0 has negative coordinates when a>0\n", "vertex = q.diff(t).solve(t)[0].rhs()\n", "print(\"The value of q(t) at the vertex t = {} is {}\".format(vertex,q.subs(t=vertex).simplify_full()))" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2*a" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Thus, since q(0)>0 the factor q(t) has no root with t>0 when a>0\n", "q.subs(t=0)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "# Determine asymptotics (this takes a *long* time) -- uncomment last line to run\n", "# First order term vanishes so second order term needed\n", "M = 2\n", "g = solve(H,vars[-1])[0].rhs()\n", "# ex,pw,se = smoothContrib(G,H,r,vars,CP,M,g)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 5.11 (Non-Finitely Minimal Critical Point)" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{4^{n}}{\\sqrt{n}} \\Bigg( \\frac{1}{2 \\, \\sqrt{\\pi}} - \\frac{1}{8 \\, \\sqrt{\\pi} n} + \\frac{1}{256 \\, \\sqrt{\\pi} n^{2}} + \\frac{5}{256 \\, \\sqrt{\\pi} n^{3}} + \\frac{819}{65536 \\, \\sqrt{\\pi} n^{4}} + O\\Bigg( \\frac{1}{n^{5}} \\Bigg)\\Bigg)\n", "\\end{math}" ], "text/plain": [ "4^n/sqrt(n) \\Bigg( 1/2/sqrt(pi) - 1/8/(sqrt(pi)*n) + 1/256/(sqrt(pi)*n^2) + 5/256/(sqrt(pi)*n^3) + 819/65536/(sqrt(pi)*n^4) + O\\Bigg( n^(-5) \\Bigg)\\Bigg)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "var('x,y')\n", "G = 1; H = (1+2*x)*(1-x-y)\n", "r = [1,1]\n", "vars = [x,y]\n", "CPs = critpt(H,r,vars)\n", "CP = CPs[0]\n", "M = 5\n", "g = solve(H,vars[-1])[0].rhs()\n", "\n", "# Determine asymptotics\n", "ex,pw,se = smoothContrib(G,H,r,vars,CP,M,g)\n", "\n", "# Print result\n", "disp_asm(ex,pw,se.expand(),M)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 5.12 (Local Central Limit Theorem for Weighted Walks)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{x: 1, y: 1/4}" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Define the function F = G/H and find the minimal critical point\n", "var('x,y')\n", "G = 1; H = 1-y*(x+1+2/x)\n", "m = [-1/4,1]\n", "vars = [x,y]\n", "CPs = critpt(H,m,vars)\n", "[CP] = [k for k in CPs if x(k)>0]\n", "CP" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{2 \\, \\sqrt{11} \\sqrt{2} 4^{n}}{11 \\, \\sqrt{\\pi} \\sqrt{n}}\n", "\\end{math}" ], "text/plain": [ "2/11*sqrt(11)*sqrt(2)*4^n/(sqrt(pi)*sqrt(n))" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# The largest coefficient of [y^n]F(x,y) as a function of x has the following asymptotic behaviour\n", "g = solve(H,vars[-1])[0].rhs()\n", "ex,pw,se = smoothContrib(G,H,m,vars,CP,1,g)\n", "asm = ex*pw*se\n", "show(asm)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "e^(-1/22*(4*ep + 1)^2*n)" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# After scaling by asm, the coefficients of [y^n]F(x,y) approach the following normal distribution\n", "var('n,ep')\n", "hes = getHes(H,m,vars,CP)[0,0]\n", "norm = exp(-n*(ep-m[0])^2/2/hes)\n", "norm" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAGGCAYAAACJ/96MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAA9hAAAPYQGoP6dpAABTiElEQVR4nO3dd3QUVR/G8e8kJDQpIlJFpEvoJYAgoqgoVYoKYuO1AIJKR4N0pEU6hC5SBYTQO6IUadK7UgVFivQaSJn3jw2RIJAQkr1bns85e8jOzmSeHSabX+69c8eybRsRERERuTcf0wFEREREXJ0KJhEREZE4qGASERERiYMKJhEREZE4qGASERERiYMKJhEREZE4qGASERERiYMKJhEREZE4qGASEZdkOaS1LMsynUVEJNlDbKspwkUkyVy8eJF06dJx8eJF01FExLPF648ytTCJiIiIxEEFk4jEafXq1dSsWZNs2bJhWRZz5syJc5tVq1ZRqlQpUqRIQe7cuRk5cmTSBxURSSIqmEQkTlevXqVYsWIMGzYsXusfOXKEatWqUbFiRbZt20aHDh34/PPPCQ0Njf9Or19PYFoRkcRn2XaChyJpDJOIF7Isi9mzZ1O7du17rvPFF18wb9489u3bF7OsadOm7Nixg/Xr199/B+fOQdOmXAoNJV1UFBfr1yftmDGQJk0ivQMRkVjiNYbpYQZ9i4jc1fr166lSpQqcOQOjR8Pvv/NKypR8u3kz4eHh+Pn5xVo/IiKC/fv3c+rUKc63b8+FzZs5Hf3aD9Onk+HECR7/+mty5MhB9uzZ/7O9iEhSU8EkIonu5MmTZE6dGkqVgmPHAMgMRAB/HzvG6XPn2Lp1K9u2bWPr1q3s2rWLsLCwu36vjwFWr4bnngPAx8eH/PnzU6xYMUqXLk3VqlUJCAhAsw+ISFJSwSQiScLasCGmWAoDVkcvL1KkCJevX8fX15eCBQtSsmRJGjZsSP78+Uln26SvUYP0QBTwJHAC8PPx4dSmTfz5zz8cPXqUXbt2sWPHDjp37ky7du3ImTMn1V56ieqRkbxw8SKpSpSATz6BjBlNvHUR8UAqmEQk0WXJkoUTR4+yGJgEzAeuRL/WskwZqvftS9GiRUmZMuV/N65SBZYt41L001RA2rp1eaxkSQLuWDUsLIyVK1eycM4cFo4bx4jwcFIAdWfPpsnw4VTcvRvrsceS5k2KiFfRVXIikqiuXr2Kv78/I48epRqwG/gCqA+UA7q3a0fZsmXvXiwBTJjgKJpuqVoVRo2666opUqTg1VdfZWjp0hwKD2cf0BXYBFQ6eZJCAQGMHTuW8PDwxHuDIuKVVDCJSJyuXLnC9u3b2b59O+CYNmD79u0ci+5yCwoK4q233qJfv37kypWLVatWEQXUT5OGaUA2YBbQtmRJRwF0P1mywNKlcOiQ4/m0aZAhw/232bcPC3gaR3H2O7ACKJgyJR9//DEFChTgu+++IyIiImEHQES8ngomEYnT5s2bKVGiBCVKlACgdevWlChRgs6dOxMVFcXKlSsJDQ0lKCiI2rVrc/DgQVasWMH+3LkpkSwZPdKmZUjjxtT79VfwiefHzoOMPypbNtZTC6gMhH75JTt37qRkyZJ88MEHFCxYkIkTJ6pwEpEHpnmYRCRu338PjRrBra6tJ5+En35ix5UrNGnShI0bN/L+++/TpUsXcuXKlSi7vHTpUsy95NKmTXv/lSMioHp1WLbs32XPPAMrVkB019/27dvp2rUrc+fOpVSpUowZMyamABQRrxavS2xVMInI/V29CtmywaVLMYuuAF3z5WPQ4cMUKFCAkSNHUrFixUTd7QMVTACRkTB/PmzZAoULQ926cJf5mjZs2EDjxo3Zu3cvrVq1olu3bqRKlSpRs4uIW1HBJCKJYPVqqFQp5uk84DPgNNC5Vy/atGmDv79/ou/2gQumBxAeHk6/fv3o1q0b2bJlY+TIkY6JNkXEG8WrYNIYJhG5vyefBMviItAAeA0IAPYEBBAUFJQkxVJS8/PzIygoiF27dpErUyZeeeUVPixYkKuLFpmOJiIuSgWTiNzfU0+x+dVXKQksBr4HFgG5u3UzmysR5Dt0iB+3bGEsMO233wisXp3drVubjiUiLkgFk4jck23bDBkyhPI//kiGHDnYVqYMb732GtbSpfD666bjPbygIKyICD4ENuOYyTdw4EDGDh3KQwxXEBEPpDFMInJXYWFh/O9//2PatGm0bNmSPn36kDx5cqftPynHMMXw9YWoqJin14FWwCjgo48+IiQkxC27HEXkgWgMk4gkzD///EPlypWZM2cOM2bMYODAgU4tlpymZMlYT1MCI9OmZfzIkUyYMIEqVapw5swZM9lExKWoYBKRWH777TfKlSvH4cOHWbVqFa97QtfbvQQHQ4oUsZf17Mn7TZrw008/sWfPHsqWLcvevXvN5BMRl6EuORGJ8dNPP1GvXj2yZ8/OwoULyZkzp7EsTumSAzh8GCZNgmvXHOOyAgNjXvrjjz+oWbMmx44dY9asWbz44otJl0NETNE8TCJyF8eOOW5me/w4vPQSvPUW+Pry3Xff0bhxYypXrswPP/xAunTpjMZ0WsEUh8uXL/Pmm2/y008/MW3aNOrUqWMsi4gkCRVMInKHffugfHm4cOHfZfXr802pUrRv354mTZowdOhQ/O4yQ7azhISEEBISQmRkJPv37zdeMAHcvHmTd999l5kzZ/Ltt9/SqFEjo3lEJFGpYBKRO7z7LkyeHGtRT6Aj0LFjR7p3745lxeuzI8m5SgvTLZGRkXzyySeMGTOGgQMH0rJlS9ORRCRxxOtDL1lSpxARF7J7d8yXNtAV6A50r1uXTj16GArlHnx9fRk1ahSPPvoorVq14vz583Tt2tVlCkwRSVoqmES8SenSsH07NvAV0BvoA3yhYileLMuib9++pE+fng4dOmBZFl27djUdS0ScQAWTiDf56ivsBQtod/Ik/YEBQKtmzSAgwHQytxIUFATXr9OhWzf8Ll3iqwEDTEcSkSSmeZhEvMlTT9GxYUP6A0NfeolWy5dDSIjpVO7n228J6t+fbkDHgQP5pnhxCA83nUpEkpAKJhEv0qdPH3oNGED//v35dPlyx7QC8mCOHoUmTeDaNTrjGDDffscOhjdsaDqZiCQhFUwiXmLEiBEEBQXRuXNnWrdubTqO+1q8GCIjY552B1oAn86cSWhoqLFYIpK0VDCJeIFZs2bRvHlzWrRooUHKDytDhlhPLRxjwernyEHDhg1ZtWqVkVgikrRUMIl4uLVr1/L222/zxhtvMGDAAF0G/7Bq1YI7bhnj4+PD+G+/pWLFitSqVYudO3caCiciSUUFk4gH+/3336lVqxZlypRhwoQJ+PjoR/6hpUgBP//suKVM1qyOmdPnziX5yy8za9Ys8uTJw6uvvsoff/xhOqmIJCLN9C3ioU6dOsUzzzxDihQpWLt2LY8++qjpSA/E1Wb6jq9Tp05Rvnx5kiVLxrp163jsscdMRxKR+4tXs7v+3BTxQFevXqVGjRpcv36dxYsXu12x5M4yZ87M0qVLOXfuHHXr1uXGjRumI4lIIlDBJOJhoqKieOedd/jtt99YtGgROe8YbyNJL2/evMyZM4cNGzbQpEkTHqIlX0RchAomEXd1/ToEB8Mrr8AHH8COHQB06tSJuXPn8v3331OiRAnDIb1XhQoVGDduHBMmTCA4ONh0HBF5SLo1ioi7qlEDfvrp3+dTpzK1Y0d69epFcHAwNWvWNJdNAHj77bfZt28fQUFBFClShGrVqpmOJCIJpEHfIu7ol1+gYsVYi3YC5Xx9eb1hQyZMmOD20we466DvO0VFRfHaa6+xZs0afv31V/Lnz286kojEpkHfIh7r0KFYTy8A9YD8/v6MHDnSrYulkJAQAgICCAwMNB0lUfj4+DB58mSyZMnCa6+9xqVLl0xHEpEEUAuTiDs6dAjy5QPbJgqoA6wCtjRqRJ7vvjMcLnF4SgvTLb///jtlypShUqVKzJkzR3NiibgOtTCJeKw8eaBzZwD6AvOAyU88QZ6+fY3GknsrUKAAU6ZMYf78+fTr1890HBF5QCqYRNxV164sHzWKjpZFp3r1qHHoEGTKZDqV3EeNGjX48rPP6NChA2t0zzkRt6IuORE3dezYMUqWLEnp0qVZuHAhvr6+piMlKk/rkuPmTfj4YyImT6ZyVBSHfH3ZNmcOmWrUMJ1MxNupS07EU4WFhVGvXj0eeeQRpkyZ4nHFkkfq0wcmTiRZVBTTgIjISN6uW5fIq1dNJxOReFDBJOKGWrZsya5duwgNDdW9ytzFDz/EfJkN+B5YER5OzyZNjEUSkfhTwSTiZqZPn86oUaMYOnQopUqVMh1H4itFilhPXwS6Al2//54VK1aYSCQiD0BjmETcyJEjRyhevDjVqlXj+++/d+v5luLicWOYRo6ETz6JtSgyXz6q5szJjp072bZtG9myZTMUTsSraQyTiCcJDw/nrbfe4rHHHnP7ySm9UtOm0LcvZM4MPj5QrRq+S5cyecoUkiVLxjvvvENUVJTplCJyDyqYRNxEly5d2LJlC1OnTiVdunSm40hCtG8PJ09CRAQsXAi5cpEpUyYmT57MypUrGThwoOmEInIPKphE3MCKFSvo06cPPXr0oGzZsqbjyMO6o3XwhRdeoE2bNnTo0IEdO3YYCiUi96MxTCIu7ty5cxQpUoSCBQuybNkyr7mlhseNYYrDjRs3KFOmDJGRkWzevJkUdwwSF5EkozFMIp6gefPmXLt2jQkTJnhNseSNkidPzpQpUzh48CBBQUGm44jIHfTpK+LCpk2bxrRp0xg+fDjZs2c3HUeSWOHChenbty+DBg1i+fLlpuOIyG3UJSfiKjZtggsXoGJFSJGC48ePU6RIEapUqcK0adNMp3M6b+uSuyUqKopXX32VPXv2sHPnTk1MKpL01CUn4hZOn4YyZRyPKlXgiSewV6zgww8/JEWKFAwfPtx0QnEiHx8fxo8fT1hYGE2aNOEh/qgVkUSkgknEtHbtHK1Lt5w9y8g6dVi6dCnjxo0jQ4YM5rIZEBISQkBAAIGBgaajGJMtWzZGjRpFaGgoEydONB1HRFCXnIh5GTPC2bMxTw8AxYH369Rh+KxZplIZ561dcrdr1KgRoaGh7Ny5k1y5cpmOI+Kp4tUlp4JJxLQCBWD/fgAigGeBs8D2bdtIXby4wWBmqWByHIOiRYuSO3dufvzxR10lKZI0NIZJxC20aBHzZV9gEzCpXDmvLpbEIW3atIwdO5aff/6ZUaNGmY4j4tVUMImY1qwZjB7NngIF6AZ8ERhIOV1SLtFeeuklmjRpQrt27Thy5IjpOCJeS11yIi4gMjKS8uXLc+nSJbZt26ZZnlGX3O0uX75M4cKFyZMnj7rmRBKfuuRE3MWgQYPYtGkT48aNU7Eksf39N2k+/phvT5/m559/ZmT9+qYTiXgltTCJGHbgwAGKFi1K06ZNdbf626iFCbBtKFECom/I2wSYCuzp3ZscX35pNJqIB9FVciKuLioqihdeeIG//vqLnTt3kjp1atORXIYKJmD9eihfPubpRaAgUDptWuZeuIBlxetzXkTuT11yIq5u1KhRrF69mrFjx6pYkv+6fDnW03TAMGD+pUvMnDnTSCQRb6UWJhFDjh49SuHChXn77bcZOXKk6TguRy1MQFgY5MgBZ87EWlwnXz7WX7rEvn37ePTRRw2FE/EYamEScVW2bdOsWTPSp09PcHCw6TjiqlKkgJkzIXv2f5fVrs2wRYu4fv067du3N5dNxMuoYBIxYObMmSxatIiQkBC3aT0ZPnw4uXLlIkWKFJQqVYo1a9bcd/0pU6ZQrFgxUqVKRdasWfnf//7H2dtuASPxVKkS/PGH436DR47A7Nlkz5uXvn37MnbsWFatWmU6oYhXUJeciJNdvHiRp59+mmeeeYZZbnKvuOnTp/Puu+8yfPhwKlSowKhRoxg7dix79+7lySef/M/6v/zyC5UqVWLgwIHUrFmT48eP07RpU/Lly8fs2bPjtU91yd1fVFQUlSpV4tSpU+zcuVPTUYgknLrkRFxRhw4duHLlCkOGDDEdJd4GDBjAhx9+yEcffUTBggUZNGgQOXLkYMSIEXddf8OGDTz11FN8/vnn5MqVi2effZYmTZqwefNmJyf3XD4+PowePZqjR4/So0cP03FEPJ4KJhEn2rBhAyNGjODrr7/miSeeMB0nXm7evMmWLVuoUqVKrOVVqlRh3bp1d92mfPny/PXXXyxatAjbtjl16hQzZ86kevXqzojsNQoWLEiHDh0IDg5m586dpuOIeDQVTCJOEh4eTpMmTShZsiSffvqp6TjxdubMGSIjI8mcOXOs5ZkzZ+bkyZN33aZ8+fJMmTKF+vXr4+/vT5YsWUifPj1Dhw69535u3LjBpUuXYj0kbl9++SX58uXj448/JjIy0nQcEY+lgknESQYPHszu3bsZNWoUvr6+puM8sDsnSbRt+54TJ+7du5fPP/+czp07s2XLFpYsWcKRI0do2rTpPb9/7969SZcuXcwjR44ciZrfUyVPnpyxY8eyadMmQkJCTMcR8Vga9C3iBH/88QeFChWicePGbnf7k5s3b5IqVSpmzJhBnTp1Ypa3aNGC7du33/UqrXfffZewsDBmzJgRs+yXX36hYsWK/P3332TNmvU/29y4cYMbN27EPL906RI5cuTQoO94at68ORMnTuS3334j++3TEIhIXDToW8SYiIiYL23bpnnz5mTIkIHu3bsbDJUw/v7+lCpViuXLl8davnz5csrfdtuO2127dg0fn9gfL7da1e71R1ry5MlJmzZtrIfEX8+ePUmVKhWtW7c2HUXEI6lgEklMq1dD6dLg5wd58sDUqYSGhrJo0SKGDh1KmjRpTCdMkNatWzN27FjGjRvHvn37aNWqFceOHYvpYgsKCuK9996LWb9mzZrMmjWLESNGcPjwYdauXcvnn39OmTJlyJYtm6m34dHSp09P//79+eGHH1i2bJnpOCIeR11yIonlxAnIlw+uXo1ZdBEo+NhjlHn2WebMmWMsWmIYPnw4wcHBnDhxgsKFCzNw4ECee+45ABo1asQff/zBypUrY9YfOnQoI0eO5MiRI6RPn57KlSvTt2/feHcXaR6mB2fbNi+88ALHjx9n165dmptJJH7i1SWngkkksQwcCHd0h7QExiZLxr7DhzWI+QGpYEqYvXv3UqxYMTp37kynTp1MxxFxBxrDJOJUtw1YBtiF487ynQsXVrEkThMQEEDbtm3p2bMnhw4dMh1HxGOohUkksRw4AAULQmQkNvA8cArYGRqKf926ZrO5IbUwJdzVq1cJCAigUKFCLFy48J7TP4gIoBYmESfLlw8mTYJMmZgGrAaGvP++iiVxutSpUzNkyBAWL14c73v3icj9qYVJJJFdPnuWpwMCKFeuHKFz55qO47bUwvRwbNumVq1a7Nixg3379pE6dWrTkURclVqYREzo0bcv5y9fZoAb3VxXPI9lWQx66y1OnzhBn/btTccRcXsqmEQS0W+//cbAgQPp0KEDOXPmNB1HvNW1a/DCC+R5+23aRUQQPHw4B93o/oUirkhdciKJxLZtqlSpwuHDh9mzZ4/mwHlI6pJ7CMHB8MUXAFwDCgJFgfm7dkHhwiaTibgidcmJONOsWbP48ccfGTx4sIolMevHH2O+TAUMABYACwcPNpVIxO2pYBJJBNeuXaNVq1bUqFGDGjVqmI4j3u6Oeb/qAi8CLRYsICwszEgkEXengkkkEfTu3ZtTp04xaNAg01FEoEULSJky5qkFDMmTh6NnzjBgwABzuUTcmAomkYd08OBBgoODad++PXny5DEdRwSKFoVffoEGDSAwENq0IWDDBlq0aEHPnj35888/TScUcTsa9C3ykGrUqMGuXbvYt28fqVKlMh3H7YWEhBASEkJkZCT79+/XoO9EdOnSJfLnz0/lypX5/vvvTccRcRW6+a5IUluwYAE1a9YkNDSUuprRO1HpKrmk8d133/HBBx/wyy+/UKFCBdNxRFyBCiaRpHTz5k0KFy7Mk08+yfLly3W/rkSmgilpREVFUaZMGQB+/fVXfHw0MkO8nqYVEElKISEhHDp0iIEDB6pYErfh4+PDoEGD2LJlCxMnTjQdR8RtqIVJJAHOnDlD3rx5adiwIcOHDzcdxyOphSlpvfXWW6xcuZL9+/eTJk0a03FETFILk0hS6dy5MwDdunUznEQkYfr27cvFixfp3bu36SgibkEFk8gD2r17N6NGjaJz5848/vjjpuOIJMiTTz5Ju3bt6N+/P4cPHzYdR8TlqUtO5AHcul/c0aNH2b17N/7+/qYjeSx1ySW9q1ev8vTTT1OmTBlCQ0NNxxExRV1yIoltwYIF/Pjjj/Tv31/Fkri91KlT07dvX2bNmsXPP/9sOo6IS1MLk8i9LFwI48dDVBS88w43q1encOHC5MyZk2XLlunKuCSmFibnsG2bChUqcPXqVbZu3Yqvr6/pSCLOFq8P82RJnULELQ0bBp999u/zWbMIqV6dQ4cOMWvWLBVL4jEsy2Lw4MGUKVOGsWPH0qRJE9ORRFySWphE7hQZCdmywenTMYv+AfIBDRs3ZvioUcaieRO1MDlXo0aNWLhwIQcOHCB9+vSm44g4k8YwiSTI5cuxiiWALtH/dmvWzPl5RJygV69eXL9+ne7du5uOIuKSVDCJ3Cl9eihUKObpbmAU0DlDBh4vWtRUKpEklS1bNr766iuGDh3K77//bjqOiMtRl5zI3axcCTVqYF+9ShXgqGWxe+5c/GvWNJ3Ma6hLzvnCwsIoWLAgRYoUYd68eabjiDiLuuREEuz55+HwYRY0bcqPQP/vvlOxJB4vRYoU9O3bl/nz52uaAZE7qIVJ5B5u3rypaQQMCAkJISQkhMjISPbv368WJiezo6KoUKIEYeHhbN69Gx8f/V0tHi9eH+4qmETuYeDAgbRt25YdO3ZQuHBh03G8jrrkDNizB+rVY/3vv1MeGJ8vH+//8gtkymQ6mUhSUsEkklBnzpwhb968NGzYkOHDh5uO45VUMBlQtCjs2gVAfWAtsL9OHVLNmmU0lkgS0xgmkYTq3r07tm3TrVs301FEnOPAgZhiCaAPjvnH+s+d65jtXsTLqWASucP+/fsZMWIEHTp04PHHHzcdR8Q5UqeG28bp5QI+B/pGRXHi1CljsURchQomkTsEBQWRLVs2Pv/8c9NRRJwnWzaoVSvWoq+AFClT0qlTJzOZRFyICiaR2/zyyy/MmjWLnj17kjJlStNxRJxr0iRo0sQxeWvWrKTv1IkuvXoxbtw4du7caTqdiFEa9C0SzbZtnnnmGcLDw9m0aZMupzZMg75dQ3h4eMz0GkuXLtX0GuKJNOhb5EHMmDGDjRs30q9fPxVLItH8/PwIDg5m+fLlLFmyxHQcEWPUwiQC3Lhxg4IFC1KoUCHmz59vOo6gFiZXYts2L7zwAv/88w87duwgWbJkpiOJJCa1MInEV0hICMeOHSM4ONh0FBGXY1kWAwYMYN++fXz77bem44gYoRYm8Xrnzp0jb9681K9fnxEjRpiOI9HUwuR63n//fZYsWcKBAwf0fyKeRC1MIvHRs2dPwsPD6dq1q+koIi6tZ8+eXL58mb59+5qOIuJ0KpjEqx0+fJihQ4fyxRdfkDlzZtNxRFzaE088QZs2bRgwYAB//vmn6TgiTqUuOfFqDRo0YM2aNRw4cIBUqVKZjiO3UZeca7p8+TL58uXj5ZdfZtKkSabjiCQGdcmJ3M+GDRuYPn06X3/9tYolkXhKkyYNPXr0YPLkyWzevNl0HBGnUQuTeCXbtqlYsSKXL19m69at+Pr6mo4kd1ALk+uKjIykePHiZMiQgZUrV2oyS3F3amESiXHhAmzdCleuADBnzhzWrl1Lv379VCy5mJCQEAICAggMDDQdRe7B19eXfv36sXr1aubNm2c6johTqIVJPN/XX0OvXnD9OqRJw80ePSg0bBh58uTRzMUuTC1Mru+VV17hjz/+YPfu3fj5+ZmOI5JQ8WphUsEknm3pUnj11ViLhgItfXzYvn07RYoUMZNL4qSCyfXt3LmT4sWLM2zYMJo1a2Y6jkhCqUtOhJkzYz29CHQDGhUpomJJ5CEVLVqURo0a0bVrVy5dumQ6jkiSUsEkni116lhPewPXge6vvGIkjoin6dGjB1euXNFkluLxVDCJZ/vf/yD6RqFHgUFA22TJyN68uclUIh4je/bsmsxSvIIKJvFsxYrBvHlQvDgdLYv0fn60W7AAnnzSdDIRj9G+fXvSpk1Lx44dTUcRSTIqmMTzVa3KlrFjmWzbdB82jEfUHSeSqNKkSUO3du2YNGkS21atMh1HJEnoKjnxeLZtU7lyZU6fPs2OHTtIFt1FJ65NV8m5ke7diejZkyI3b5Ldx4flPXtiffml6VQi8aWr5EQAFi5cyMqVKwkODlaxJJLYVqyALl1IdvMm3wAroqJYEhQEv/xiOplIolLBJB4tIiKCdu3aUblyZapVq2Y6jojnmTUr5svqwPNAWyAiNNRQIJGkoYJJPNrYsWP5/fff6devn+53JZIUbusutYB+wF7gu4MHTSUSSRIawyQe6/Lly+TNm5dXXnmFiRMnmo4jD0hjmNzE/v2Oq1HDwmIWvePry4oMGThw+DCPPPKIwXAi8aIxTOLdgoODuXTpEl9//bXpKB5h+PDh5MqVixQpUlCqVCnWrFlz3/Vv3LjBV199Rc6cOUmePDl58uRh3LhxTkorTpM/PyxZAhUqwCOPQMWK9Jw2jfOXLtGvXz/T6UQSjUbAikc6fvw4/fv3p1WrVjypOZce2vTp02nZsiXDhw+nQoUKjBo1iqpVq7J37957Ht8333yTU6dO8e2335I3b15Onz5NRESEk5OLU1SqFGuQd06g5ebNfPPNNzRp0oSsWbOayyaSSNQlJx7pgw8+YP78+Rw8eJB06dKZjuP2ypYtS8mSJRkxYkTMsoIFC1K7dm169+79n/WXLFlCgwYNOHz4MBkyZEjQPtUl594uXrxInjx5qFOnDmPGjDEdR+R+1CUn3mnHjh2MHz+erl27qlhKBDdv3mTLli1UqVIl1vIqVaqwbt26u24zb948SpcuTXBwMNmzZyd//vy0bduW69ev33M/N27c4NKlS7Ee4r7SpUtHly5dGDduHLt37zYdR+ShqWASj9OuXTvy5ctH48aNTUfxCGfOnCEyMpLMmTPHWp45c2ZOnjx5120OHz7ML7/8wu7du5k9ezaDBg1i5syZNL/PPfx69+5NunTpYh45cuRI1PchztekSRNy585N+/btTUcReWgqmMSjLF26lOXLl9O3b1/8/PxMx/Eod07LYNv2PadqiIqKwrIspkyZQpkyZahWrRoDBgxg/Pjx92xlCgoK4uLFizEP3cjV/fn7+9O3b18WL17Mjz/+aDqOyENRwSQeIzIykrZt21KxYkVee+0103E8RsaMGfH19f1Pa9Lp06f/0+p0S9asWcmePXusLtGCBQti2zZ//fXXXbdJnjw5adOmjfUQ91enTh0qVKhA27ZtiYyMNB1HJMFUMInHGD9+PLt379YklYnM39+fUqVKsXz58ljLly9fTvny5e+6TYUKFfj777+5cuVKzLL9+/fj4+PDE088kaR5xbVYlkW/fv3YsWMHkydPNh1HJMF0lZx4hKtXr5IvXz4qVarE1KlTTcfxONOnT+fdd99l5MiRPPPMM4wePZoxY8awZ88ecubMSVBQEMePH4+ZIPTKlSsULFiQcuXK0a1bN86cOcNHH31EpUqV4n3FlK6S8yz169dn7dq17N+/n1SpUpmOI3I7XSUn3qN///6cPXuWXr16mY7ikerXr8+gQYPo3r07xYsXZ/Xq1SxatIicOXMCcOLECY4dOxaz/iOPPMLy5cu5cOECpUuX5u2336ZmzZoMGTLE1FsQw3r37s3p06cZNGiQ6SgiCaIWJnF7J06cIF++fHzyySd88803puNIIlELk+dp3bo1Y8eO5eDBg2TKlMl0HJFb1MIk3qFz584kT56cDh06mI4iIvfRsWNHfH196datm+koIg9MBZO4l4gIWLsWtm0DYNeuXYwbN44uXbrw6KOPGg4nIveTIUMGOnbsyKhRo/jtt99MxxF5IOqSE/exaRPUrQu3LksPDKTqI49w8M8/2bNnD/7+/mbzSaJSl5xnunHjBk8//TRFixZl7ty5puOIQDy75HTzXXEPUVHQoMG/xRKwbNMmlgChoaEqlkTcRPLkyenduzdvvfUWq1atolKlSqYjicSLWpjEPezcCcWKxTyNBEoAaS2LNZGRmnfJA6mFyXPZtk25cuWIiopi48aN+PhodIgYpUHf4kHSp4/1dAKwC+ifIYOKJRE3c2syy82bNzNt2jTTcUTiRQWTuIcnn4RatQC4CnQEGgBlW7c2mUpEEqhixYrUrl2bDh06EBYWZjqOSJxUMIn7mDIFWrSgf5o0nLUsegUFQVCQ6VQikkB9+vThr7/+YujXX5uOIhInjWESt6JJKr2HxjB5uHPn4L33+HThQiYDh8qV47EffoAcOUwnE++jMUzieTRJpYiHaNYMFi6kCxAF9NiwAd55x3QqkXtSwSRu49YklZ07d9YklR4sJCSEgIAAAgMDTUeRpHLzJoSGAvA4EASEAAdXr4bjx00mE7kndcmJ26hatSoHDx7UJJVeQl1yHiw8HNKmhejB3teB/EA5YMbJk5A5s8l04n3UJSeeY9myZSxZsoS+ffuqWBJxd35+8O67MU9TAj2BmcC6Q4dMpRK5L7UwicuLjIykZMmSpEmThjVr1mjeJS+hFiYPd+0atGwJkyZBZCRRdetSau9eUj7yCGvXrtXPuTiTWpjEM0yYMIGdO3fSv39/fYiKeIpUqWD0aLh8Ga5exWfaNPoPGsT69esJjR7fJOJK1MIkLu3q1avky5eP5557TjMCexm1MHmn6tWr8/vvv7N37151v4uzqIVJ3F///v05e/YsvXv3Nh1FRJwgODiYI0eOMGLECNNRRGJRwSQu68SJEwQHB/P555+TK1cu03FExAkKFSrEhx9+SPfu3blw4YLpOCIxVDCJy+rSpYsmqRTxQt26dePGjRv06tXLdBSRGCqYxCXt3r2bb7/9VpNUinihrFmz0q5dOwYPHswff/xhOo4IoEHf4qI0SaVo0Ld30wUf4kQa9C3uaenSpSxZsoQ+ffqoWBLxUqlTp6ZXr15Mnz6ddevWmY4johYmcS0REREUK1aMjBkzsnLlSs275MXUwiRRUVEEBgbi6+vLhg0b8PHR3/iSJNTCJO5n9OjR7Nu3j4EDB6pYEvFyPj4+DBo0iE2bNjFlyhTTccTLqYVJXMb58+fJly8ftWrVYty4cabjiGFqYZJb3nzzTdauXcv+/ftJnTq16TjiedTCJO6lR48ehIWF0bNnT9NRRMSF9O3bl7NnzxIcHGw6ingxFUxixvHj0KmT447l337L/t27GTp0KB06dCBr1qym04mIC8mVKxetW7cmODiYY8eOmY4jXkpdcuJ8hw9D2bJw5kzMolqZMrEzZUr27dtHypQpDYYT00JCQggJCSEyMpL9+/erS04AuHz5Mvny5ePFF1/UeCZJbOqSExf1zTexiqUfgfmnTxPcqJGKJaF58+bs3buXTZs2mY4iLiRNmjT06tWL77//nvXr15uOI15IBZM43549MV9GAK2ACsAbjz9uKpGIuIH333+f4sWL06pVK6KiokzHES+jgkmcLzAw5stvgd3AQMC6bbmIyJ18fX0ZOHAgGzduZOrUqabjiJdRwSTO17Yt5MrFRaAj8B4Q+P77UKaM4WAi4uqeL1+euqVK8UWTJlzt1QvOnTMdSbyECiZxvqxZYft2vn7pJa75+dFr0iT47jvTqUTE1UVFQY0aBG/Zwj9Xr9Lvq6+gZEk4dcp0MvECKpjEiIOnTzN41Sq+7NSJ7O+8A5rVW0TismQJLF9OHqAl0Bf46+hRGDLEbC7xCiqYxIj27duTJUsW2rRpYzqKiLiLHTtivvwKSAMEAezcaSiQeBMVTOJ0P//8M7Nnz6Zv376kSpXKdBwRcRclSsR8mRb4GpgM/Jopk6lE4kU0caU4VWRkJKVKlSJlypSsW7dON9iVe9K95OQ/oscwsXgxAJFASX9/Uhctytpff9XniSRUvE6cZEmdQuR2Y8aMYceOHWzYsEEfbiLyYHx8YN48mDkT1q7FN39+Bj31FJVr1WL69Ok0aNDAdELxYGphEqc5d+4c+fLlo1atWnynq+IkDmphkviqU6cOW7ZsYd++faROndp0HHE/ujWKuJZOnToRERFBnz59TEcREQ/Sr18/Tp06Rd++fU1HEQ+mgkmcYseOHYwcOZIuXbqQOXNm03FExIPkyZOHtm3bEhwczJEjR0zHEQ+lLjlJcrZt8/zzz3P69Gl27NiBv7+/6UjiBtQlJw/i6tWrFChQgMDAQGbPnm06jrgXdcmJa5g+fTqrV69myJAhKpZEJEmkTp2afv36MWfOHJYtW2Y6jnggtTBJktJffZJQamGSB3V7a/bOnTvx8/MzHUncg1qYxLxevXpx5swZBgwYYDqKuImQkBACAgIIDAw0HUXcjGVZDBkyhP379zN06FDTccTDqIVJksyhQ4cICAjgiy++oHv37qbjiJtRC5MkVPPmzZk0aRL79+8nS5YspuOI64tXC5MKJkkytWrVYvv27fz222+6BYo8MBVMklCa800ekLrkxJzFixczf/58+vfvr2JJRJwqQ4YM9OzZk/Hjx7Nx40bTccRDqIVJEt3NmzcpUqQI2bNnZ8WKFboFiiSIWpjkYURGRlK6dGn8/PzYsGEDPj5qH5B7UguTmDF48GAOHTrE4MGDVSyJiBG+vr4MHTqUTZs28e2335qOIx5ABZM8nKNHoXZtSJ4ccuTg727d6N69O82aNaNIkSKm04mIF3v22Wd5//33+eKLL/jnn39MxxE3py45SbioKAgIgN9/j1n0FrAiTRp+P3qURx991Fw2cXvqkpPEcPr0aQoUKEDdunXV0iT3oi45SWIrV8YqlpYD04B+jz+uYklEXEKmTJno3bs348aNY926dabjiBtTwSQJd+1azJdhQHOgEvBuihSmEomI/MfHH39MYGAgn3z8MRG3/ZEn8iBUMEnCVa4M0S1J3wBHgOGA9cYbJlOJiMTiGx7OiMceY9fevQx7+mkoXBh27DAdS9yMCiZJuFSpYPZsDmXJQk+gDRDQsCEEBZlOJiLyry5dKLVkCc2ATsDxPXvgtdcgMtJ0MnEjKpjkodjPPcenxYqRJWtWOh04AFOmOK6YE48zfPhwcuXKRYoUKShVqhRr1qyJ13Zr164lWbJkFC9ePGkDitzLtGkAfA2kwvHHHUePwvr1BkOJu1HBJA8lNDSUJUuXMmTkSFLnzWs6jiSR6dOn07JlS7766iu2bdtGxYoVqVq1KseOHbvvdhcvXuS9997jxRdfdFJSkbuI/iMuPdAPmI7jIhX9cScPQtMKSIJdvnyZggULUqpUKebOnWs6jiShsmXLUrJkSUaMGBGzrGDBgtSuXZvevXvfc7sGDRqQL18+fH19mTNnDtu3b4/3PjWtgCSaPn1ihgrYwPPACX9/dl26RHIVTaJpBSSpde3alXPnzjF48GDTUSQJ3bx5ky1btlClSpVYy6tUqXLfy7S/++47Dh06RJcuXeK1nxs3bnDp0qVYD5FE0b49dOwI6dNj+fgw/LnnOBIVxTfffGM6mbgRFUySIDt37mTw4MF07tyZp556ynQcSUJnzpwhMjKSzJkzx1qeOXNmTp48eddtDhw4wJdffsmUKVNIlixZvPbTu3dv0qVLF/PIkSPHQ2cXAcDHB3r0gHPnICyMQqtW0bp1a3r27Mnhw4dNpxM3oYJJHlhUVBSffPIJBQoUoHXr1qbjiJPceV9A27bveq/AyMhIGjZsSLdu3cifP3+8v39QUBAXL16Mefz5558PnVkkFssCPz8AOnXqxOOPP85nn33GQwxNES8Svz/9RG7z3XffsW7dOlauXIm/v7/pOJLEMmbMiK+v739ak06fPv2fVidwjG3bvHkz27Zt49NPPwUcRbZt2yRLloxly5ZRuXLl/2yXPHlyjScRp3nkkUcYPHgwdevWZe7cudSuXdt0JHFxGvQtD+TMmTMUKFCA6tWrM3HiRNNxxEnKli1LqVKlGD58eMyygIAAXnvttf8M+o6KimLv3r2xlg0fPpyffvqJmTNnkitXLlKnTh3nPjXoW5KabdvUqFGDXbt2sXfvXh555BHTkcSMeA36VguTPJAvv/ySqKgo+vXrZzqKOFHr1q159913KV26NM888wyjR4/m2LFjNG3aFHB0px0/fpyJEyfi4+ND4cKFY22fKVMmUqRI8Z/lIiZZlsXQoUMpXLgwXbp0oX///qYjiQtTwSTxtmrVKr799ltGjBhBpkyZTMcRJ6pfvz5nz56le/funDhxgsKFC7No0SJy5swJwIkTJ+Kck0nEFeXOnZuuXbsSFBREw4YNKVWqlOlI4qLUJSfxEhYWRtGiRcmUKROrV6/Gx0fXC0jSUpecOEt4eDiBgYH4+Pjw66+/xvvKTvEYmodJEk/37t05evQoY8aMUbEkIh7Fz8+PMWPGsH37ds0rJ/ek33wSpx07dhAcHEzHjh0pWLCg6TgiIokuMDCQzz//nM6dO3PkyBHTccQFqUtO7isiIoJy5cpx48YNtmzZomkExGnUJSfOdvnyZQoVKkShQoVYtGjRXecZE4+kLjl5eIMHD2br1q2MHTtWxZKIeLQ0adIQEhLCkiVLmDZtmuk44mLUwiT3dOjQIYoUKUKTJk0YOHCg6TjiZdTCJKa8+eabrFy5kr1795IxY0bTcSTpxauFSQWT/GvRIhgyBM6cwa5enZdXr+bgkSPs3r1bE7qJ06lgElNOnjxJQEAA1apVY/LkyabjSNLTxJXyAGbPhnr1ILqAnrBlCyuAJUuWqFgSEa+SJUsWBg8ezHvvvUeDBg2oUaOG6UjiAtTCJA7ly8P69QCcAgoCNYCJO3dCkSImk4mXUguTmGTbNtWrV2fnzp3smTSJdHnzQo4cpmNJ0tCgb3kAx4/HfPk5jqbHgQB//WUokIiIOZZlMeqtt7j099+0q1wZcuaEt96CGzdMRxNDVDCJw0svATAP+AEYDDyWOjVUqGAylYiIGZcukaN5c4JtmzHACtuGadPgm29MJxNDVDCJQ8+eXHj6aZoB1YEG/v4wejSoK0ScLCQkhICAAAIDA01HEW+2fDlcvkxj4HngY+AqQGioyVRikAomcciShVZly3I5VSpGDBmC9ddf0LCh6VTihZo3b87evXvZtGmT6SjizdKkARy/JMcCJ4Gvblsu3kcFkwCwYMECxk+YwKBhw8jx2Wfw+OOmI4mImPPii5AvHwB5gK+BIcDaypVNphKDdJWccO7cOQoXLkyJEiVYsGCBbgcgLkFXyYlxR49C69awdCmRWbJQ0bb5x8eH7du3kzp1atPpJPHoKjmJnxYtWnD9+nVGjx6tYklE5JacOR1jlq5cwffgQcYvXszx48f58ssvTScTA1Qwebk5c+YwefJkhgwZQvbs2U3HERFxWfnz56dv374MGzaMFStWmI4jTqYuOS925swZChUqRLly5ZgzZ45al8SlqEtOXFFUVBQvvfQShw4dYufOnaRLl850JHl46pKT+/v000+JiIhg1KhRKpZEROLBx8eHcePGcf78eVq3bm06jjiRCiYvNWPGDKZPn86wYcPIkiWL6TgiIm7jqaeeYuDAgYwbN4758+ebjiNOoi45L/T3339TpEgRXnjhBWbMmKHWJXFJ6pITV2bbNjVq1GDLli3s2bOHxx57zHQkSTh1ycl/2bbNBx98gL+/PyNHjlSxJCKSAJZlMWbMGG7evEnz5s1NxxEnUMHkZUaMGMHSpUsZN24cGTNmNB1HRMRtZcuWjZCQEKZPn87UqVNNx5Ekpi45L/L7779TokQJ/ve//xESEmI6jsh9qUtO3IFt2zRs2JDFixezY8cOcubMaTqSPLh4dbWoYPIS4eHhVKhQgYsXL7J161bNUisuTwWTuIsLFy5QrFgxcubMyc8//4yvr6/pSPJgNIZJ/tWjRw+2bt3KpEmTVCyJiCSi9OnTM3nyZNauXUvv3r1Nx5Ekksx0AEkCtg3Ll8Mvv0Du3KzKlo2ePXvStWtXypQpYzqdiIjHqVixIh06dKBr1668/PLLlC1b1nQkSWTqkvNE774LkycDcBYoniwZecqUYcXq1WoqFrehLjlxN+Hh4VSsWJF//vmH7Rs3kiYyEjJnNh1L4qYuOa/0yy8xxZINfARci4hgcmCgiiVxCyEhIQQEBBAYGGg6isgD8fPzY8qUKZz+6y8+z5YNsmSBwoVh3TrT0SQRqGDyNBs3xnw5EpgDfAs8sW+foUAiD6Z58+bs3buXTZs2mY4i8sDyrFnDsJs3GR8ezg8Ae/ZA9epw+bLpaPKQVDB5mqefBmA30BpoBtS+bbmIiCShSZN4D6gPNAGOAVy4ALqFittTweRpqlblevnyNADyAv0AHn8cWrUym0tExEtYwAggDfA2EAGguyq4PRVMnsbHhzaFCnHIz4+pVaqQMigItm6Fp54ynUxExPO99x4AjwJTgfVApxQpoEYNk6kkEahg8jCzZ89mxJgxDBwyhMJLl0KvXvDEE6ZjiYh4h/ffh+7dIX16KgC9smShT1gYi9asMZ1MHpKmFfAgR44coWTJkrzwwguEhobqxrri1jStgLi1mzfh0iWiMmSgVq1arF+/nu3bt5MjRw7TyeS/dGsUbxIWFkaFChU4f/48W7Zs4dFHHzUdSeShqGAST3H27FlKlCjBE088wapVq/Dz8zMdSWLTPEzepEWLFuzZs4fQ0FAVSyIiLuSxxx7jhx9+YNOmTXTo0MF0HEkgFUweYOLEiYwePZphw4ZRokQJ03FEROQO5cqVIzg4mH79+jFv3jzTcSQB1CXn5nbt2kXZsmWpX78+48aN07gl8RjqkhNPY9s2devWZeXKlWzbto2ndPWyq9AYJk936dIlSpcuTcqUKVm/fj2pUqUyHUkk0ahgEk90/vx5SpYsSaZMmVizZg3+/v6mI4nGMHk227b54IMPOHXqFKGhoSqWRETcwKOPPsoPP/zAtm3baNu2rek48gBUMLmpQYMGERoayvjx48mbN6/pOCIiEk+BgYEMGjSIoUOHMnHiRNNxJJ7UJeeG1q5dy/PPP0/Lli355ptvTMcRSRLqkhNPZts2H330EVOmTGHNmjUEBgaajuTNNIbJE508eZJSpUqRJ08efvrpJ5IlS2Y6kkiSUMEknu7GjRtUqlSJv/76iy1btpA5c2bTkbyVxjB5hIgIWLoU5s3jxvnz1K1bF9u2mT59uoolERE3ljx5cmbNmkVkZCSvv/46N2/eNB1J7kMFkyvbvx/y54dXX8V+7TU+yZKFrVu2MGfOHLJmzWo6nUiSCAkJISAgQF0U4hWyZcvGrFmz2LhxIy2KFIE0aSBrVsf96KKiTMeT26hLzpVVqQLLlwMwGGgJTHz8cd49dQo035J4OHXJiTcZ+/TTfPz774wCGt9a2KMHdOxoMJXX0BgmtxYZCX5+YNssB14FWgH9AH77DQoUMBpPJKmpYBKvsX8/FChAM2AssBIoD5A9O/z1l8lk3kJjmNyary88/jgHgfrAy0BfgGTJIGNGo9FERCQRXb0KwCCgHFAPOA5w5YqxSPJfKphc2KXmzakFZASmAr4AjRrBY4+ZjCUiIompWDHIkwd/YAaQDKgLhNWubTSWxKaCyUVFRkby9qZNHE+ZknlFi/Jo8eLQpw+MGGE6moiIJCYfH5g5E/LmJTMwG9jp48P/Ll8mSgO/XYauS3dRbdq0YdGiRcyfP5+nq1UzHUdERJJS8eKOsUy7dlH6kUeYvG0bb7zxBrk7daJnz56m0wkqmFzS0KFDGTx4MCEhIVRTsSQi4h0sC4oWBaBe7twEBwfTrl07cuXKxUcffWQ4nKhgcjHz58+nZcuWtG7dmmbNmpmOIyIihrRp04bDhw/TtGlTnnzySapUqWI6klfTtAIuZOvWrVSsWJFXXnmFGTNm4OvrazqSiDGaVkAEIiIieO2111izZg1r166lSJEipiN5Is3D5E7+/PNPypYtyxNPPMHKlStJlSqV6UgiRqlgEnG4cuUKzz33HP/88w8bN24kW7ZspiN5Gs3D5C7Onz9PtWrV8Pf3Z/78+SqWREQkxiOPPMKCBQsAqFGjBlc0P5MRKpgMu3btGjVr1uTvv/9m8eLFulu1uKzhw4eTK1cuUqRIQalSpVizZs091501axYvv/wyjz/+OGnTpuWZZ55h6dKlTkwr4lmyZcvGokWLOHjwIA0aNCAiIsJ0JK+jgsmg8PBw3nzzTbZv386iRYsoWLCg6UgidzV9+nRatmzJV199xbZt26hYsSJVq1bl2LFjd11/9erVvPzyyyxatIgtW7bwwgsvULNmTbZt2+bk5CKeo0iRIsycOZOlS5fSuHFjHmJIjSSAxjAZEhUVRaNGjZg2bRoLFizQ1Q/i0sqWLUvJkiUZcdvEqQULFqR27dr07t07Xt+jUKFC1K9fn86dO8drfY1hErm7KVOm8M4779CuXTuCg4NNx/EE8RrDpGkFnGnDBti2DbtoUdqGhjJ58mSmTp2qYklc2s2bN9myZQtffvllrOVVqlRh3bp18foeUVFRXL58mQwZMiRFRBGv8vbbb3P27FlatGhBxgsXaL93L+zbB2XKQO/ejkkwJdGpYHKGqCho2BCmTwccN9EdCAwbMoT69esbjSYSlzNnzhAZGfmf8XWZM2fm5MmT8foe/fv35+rVq7z55pv3XOfGjRvcuHEj5vmlS5cSFljEC3z++ef8s2cPX4weTUbgA4AlSxx/mP/+O2TKZDih59EYJmeYPTumWBoLBAFdgOYa4C1uxLJit1rbtv2fZXczdepUunbtyvTp08l0nw/x3r17ky5duphHjhw5HjqziCfr7udHE+BjYNathRcuwJQpxjJ5MhVMzrByJeA4oZsAzXEUTPz0k7FIIvGVMWNGfH19/9OadPr06Tiv6pw+fToffvghP/zwAy+99NJ91w0KCuLixYsxjz///POhs4t4MuviRUKAN4AGwMJbL5w/byyTJ1PB5Aw5c7IAxwn9JjCE6BFmTz1lMJRI/Pj7+1OqVCmWL18ea/ny5cspX778PbebOnUqjRo14vvvv6d69epx7id58uSkTZs21kNE7qN2bXyBSUB1oB6wPHq5JD4VTE6w9MknqQfUBCYSfdCzZIEPPjCaSyS+WrduzdixYxk3bhz79u2jVatWHDt2jKZNmwKO1qH33nsvZv2pU6fy3nvv0b9/f8qVK8fJkyc5efIkFy9eNPUWRDxPvXrQqhV+yZIxDXjRx4fX/PxYdfmy6WQeSdMKJLEVK1ZQo0YNXn72WWbmy4f/rl1QrBi0awc5c5qOJxJvw4cPJzg4mBMnTlC4cGEGDhzIc889B0CjRo34448/WBnd/fz888+zatWq/3yP999/n/Hjx8drf5pWQCSe/v4bDhwgLF8+ar7/PuvXr2fZsmX3bQGWWHQvOdNWrVpF1apVqVSpEnPmzCF58uSmI4m4DRVMIg/u6tWrVKtWje3bt7NixQpKly5tOpI70L3kTFq5ciXVq1enQoUKzJo1S8WSiIgkudSpU7NgwQICAgJ4+eWX+fXXX01H8hgqmJLA4sWLqVq1KhUqVGDu3LmkTJnSdCQREfESadKkYenSpQQEBPDSSy+xdu1a05E8ggqmRDZr1ixee+01XnnlFebNm0eqVKlMRxIRES+TNm1ali5dSqlSpXjllVdixhdKwqlgSkTff/89b775JvXq1WPGjBnqhhMREWMeeeQRFi5cSPny5alWrRo//vij6UhuTQVTIhk7dizvvPMO7733HpMnT8bPz890JBER8XKpUqVi3rx5vPDCC9SoUYMFCxaYjuS2VDAlgsGDB/Pxxx/TrFkzxo4di6+vr+lIIiIiAKRIkYJZs2ZRrVo1ateuzfgxY6B9e8ia1fFo3x5uu4+j3J2mFUiIiRNh2jTw86N3unR0mDSJ9u3b06dPn3jdW0tE4qZpBUQSV2RkJM2aNWP06NH0Adpz2/X0jRvDqFHmwpmleZiSxFdfQa9eRAFtgYFA1+rV6Tx/voolkUSkgkkk8dlXrtAtfXq6RUbSAhhAdFdT8uRw9iykTm02oBnx+uWdLKlTeJTLl2HgQG4AjYDpwDCg+bZtRmOJiIjEh3XzJl0jI8kCNANOA+MB/xs3HN1y3lkwxYvGMD2IU6e4dP061YDZwAygOTimpb9502g0ERGROGXIAJUq0RTH77BQHDfuvVyhguM1uScVTA/gLz8/nvPzYyuOO0LXu/VC2bKO5kwRERFX9913UKQI9YBlwK8+Prxw8SKnT582ncylqWCKp02bNlHmmWc4nz49a/z9qXjrhUcfhaFDTUYT8SghISEEBAQQGBhoOoqIZ8qVC3buhG3bqLR9O6u3buX4mTM888wz7N2713Q6l6VB3/EwY8YM3nvvPYoXL86cOXPIbFkwZw74+0OdOpAunemIIh5Hg75FnOePP/6gZs2aHDt2jOnTp/Pqq6+ajuRMuvnuw7Jtm6+//po333yTOnXq8PPPP5M5c2bIlMlxCWajRiqWRETE7T311FOsW7eO5557jurVqzNo0CAeokHFI6lguoewsDDeffddOnXqRLdu3ZgyZQopUqQwHUtERCRJpEmThjlz5tCmTRtatWpFkyZNuKkLmmJoWoG7+PPPP3n99dfZuXMn06ZNo379+qYjiYiIJDlfX1+Cg4MJCAigcePG7N+/nx9++IFMmTKZjmacWpjusGLFCkqWLMnJkydZs2aNiiUREfE6jRo14qeffmLfvn2ULFmSdevWwZYt8PLLkD49PPMMLF9uOqZTadD30aNw+DBRxYrRd9QoOnbsyEsvvcSUKVPImDGj6XQiXkuDvkXMO378OPXr12fjxo30S5aMz8PC/h0h7ecHmzdD0aImIyYGDfq+r6go+PhjyJ2bi5UrUzdTJjp06ECHDh1YtGiRiiUREfF62bNn5+eff+az8uVpGRZGA+DyrRfDw2HMGIPpnMt7C6bx42HsWHZERVEaWBkZyTygx7vv4uvraziciIiIa/Dz82PA888zA1gMBAJ7br148aKxXM7mtQVT1Jw5DALKAKmBzUBNgPnzDaYSERFxQbVr8zqO35V+OH53TgTs2rVNpnIqryyYTp48SbUtW2iF4+aDG4C8t1589FFjuURERFxSiRIwaBD5U6ViA/C6ZfE+0GD6dM6dO2c6nVN4XcG0cOFCihYtyvawMBb7+DAQiJldKXNmeOMNg+lERERcVIsWcPw4qdesYcKJE0ybNo1ly5ZRtGhRfvzxR9PpkpzXFEzXr1/n008/pUaNGpQpU4ad+/bx6sKFjksjM2aE2rVh5UpIk8Z0VBEREdeUPj08+yxkzkz9+vXZtWsXBQoU4OWXX6ZVq1aEhYWZTphkvGJagc2bN9OoUSMOHjxIv379aN68OZYVr6sIRcQQTSsg4h6ioqIYPHgwQUFB5M2blylTplCsWDHTsR6EphW4du0abdu2pWzZsvj5+bF582Y+/fRTFUsiIiKJxMfHh1atWrF582Z8fX0JDAykT58+REREOFZYtAhq1YIXXoBhwyAy0mzgBPLYFqYVK1bQuHFj/v77b7p27Urr1q3x8/MzHUtE4kktTCLu58aNG3Tq1In+/ftTtGhRxtarR6lOnWKv9OGHMHasmYB3F69WFM8omG7cgP79Yf58zqdJQxvgu+XLef755xk9ejT58uUznVBEHpAKJhH3tWnTJj7++GN27dhBS6A7jil8APD1hWPHIFs2Y/nu4EUFU5062HPmEAp8ClwH+jVtyochIfj4eHSvo4jHCQkJISQkhMjISPbv36+CScRNhYeHMyB1arqGh5MFGAm8cuvFjRuhTBlz4WLzkoLpt984UrAgrYC5QG0gBMhWpozjP0RE3JJamEQ8QI0aHFy4kCbAT0B94Jv06clx4gSkSBHHxk7j+YO+r169SqcePSiIY/bRGcAsIBvAn3+ajCYiIiL9+5M3e3Z+BCYAK4EC167R45tvuH79utlsD8gtCybbtpk2bRpPP/0034SG0s7Pj9+B17mtTHzpJXMBRUREBAoUgAMHsKZP572xY9l/4ACftmhBjx49CAgIIDQ0lIfo6XIqtyuY1q9fz3PPPcdbb71FYGAge/fupcf48aT29/93pXz5oFcvcyFFRETEIWVKePNN+PBD0ubNS3BwMLt376ZQoUK8/vrrvPjii+zYseO/2129CqtXw5Ejzs98F25TMO3du5fatWtTvnx5Ll++zLJly5g1axa5c+eGhg0dI+4nTHDcPHfvXnjiCdORRURE5C7y58/PggULWLRoEcePH6d48eK89dZbHDhwwLHC9OmQPTtUqgR58jh+z4eHG83smoO+z5+Hs2chTx7+/OsvunTpwoQJE8iZMydff/01DRo00NVvIh5Og75FvEN4eDgTJkygW7dunDhxgg/eeovO06fzxJ0FUr9+0KZNUkRww6vkIiLgs8/g2285HR5O3/TpCbl2jbTp0tGpUyeaNGmC/+1dbyLisVQwiXiXsLAwRowYQa/Onbl85QrNgS+Bx2+t8OyzsGZNUuzaDa+SGziQkyNH0iY8nKeAMRcuEOTry6Fdu/jss89ULImIiHioFClS0KpVKw6NHEkHYAyQG/gK+AccN/41yGUKpoMHD/JJcDC5gLFAW+APoMv166TRfEoiIiJeIe0bb9D5qac4AnwCDAZyAi18fPjT4JRBxgumjRs3Uq9ePfLnz8+sCxfoiKNQ6g5kuLWSWpZERES8g78//PQTj73xBsEZMnC0cGG+eP11Jq1ZQ+7cufnggw/47bff7r396dOwe3ei3+TXSMEUFRXF/Pnzee655yhXrhy7d+9m5MiRHB08mK+AR29fOWdOePFFEzFFRETEhFy54Icf4OxZHtu1iy4zZnD06FH69OnDkiVLKFiwIFWrVmXZsmX/zuMUHg4ffOC4R12RIpA7N/z8c6JFSvqC6dQpx81xgfPnzzNo0CAKFixIrVq1iIiIYPbs2ezbt4/GjRuTolkz6N0bMmZ0bFu5MixZAn5+SR5TREREXFeaNGlo06YNR44cYcKECZw8eZJXXnmFwoULM3r0aK727g3fffdvy9KxY1C3Lly5kij7T7qr5NauhSZNYM8etqROzYgCBfh+3z4iIiKoV68en376KRUqVLjHd7YdlaK64kS8lq6SE5H7sW2bNWvWMGjQIObOncsjwHtRUTQFCt2+4qxZUKfO/b6VwWkFLl6EnDk5d/EiVYFfgRxA03r1+DAkhMyZMyd0nyLiJVQwiUh8HT16lDHPPsvYv/7iFFARWAKkAli8GF599X6bG5xWYP58uHiRR4HSwFzgCNDhwgUVSyIiIpKocubMydc9enAM+AFH7ZEKHGOhEuneskk6hskCQoBagC+AFa8iTkS8WEhICAEBAQQGBpqOIiLupFEj/AcN4o3s2Rng4wNVq8KyZZAsWaJ8+6Tpkrt0yXF124ULsZdPnAjvvpvQ/YmIF1GXnIg4icEuubRpYdEiKFrU8Tx9evj6axVLIiIi4paS/l5yZ844Cihd8SYiD0AtTCLiJPFqYUqcjr37uTWnkoiIiIibMn5rFBERERFX9zBdciIiScayrLTARSCdbduXTOcREe+mgklEXJJlWRaQBrhs64NKRAxTwSQiIiISB41hEhEREYmDCiYRERGROKhgEhEREYmDCiYRERGROKhgEhEREYmDCiYRERGROKhgEhEREYmDCiYRERGROKhgEhEREYmDCiYRERGROCRLyEa33eNJRERExN3Fec/KBBVMOIqliwncVkRERMSVpAMu3W+FBN189wFbmNIAfwFPAJcfYDe/AmUeMFpCt3PmvnQ8/pXQY5GQfSV0G2du58nHw5k/Kwndn6ceD1d/XwndTsfjX87+WUnodq58PJKmhSn6m963ErvFUVvFhInXNtHbRT3I+g+znZP3detLrz8eCT0WCdlXQrdx5naefDyc+bPyEPvzyOPh6u/rIfZ160uvPx7O/llJ6HbucDzux5UHfYc4cTtn7iuhdDwefl/Ofl86Hg+/r4TS8Xj4fbn6MUwoHY/E2ZenHo97SlCX3APtwLLS4hjvlC4xKz13pePxLx2L2HQ8YtPxiE3HIzYdj3/pWMSWVMfDGS1MN4Bu0f+KjsftdCxi0/GITccjNh2P2HQ8/qVjEVuSHI8kb2ESERERcXeuPIZJRERExCWoYBIRERGJgwomERERkTioYBIRERGJQ5IVTJZlPWVZ1reWZR2xLOu6ZVmHLMvqZlmW/x3rPWlZ1nzLsq5alnXGsqwhd67jCSzL+sqyrHWWZV2zLOvCPdax7/Jo6uSoThHP4+EV58bdWJb1x13OhT6mczmLZVnNoj87wizL2mJZVkXTmUywLKvrXc6Dk6ZzOYtlWc9Ffwb8Hf3ea9/xuhV9jP6O/j2z0rKsQobiJrl4HI/xdzlfNhiKm6QsywqyLGuTZVmXLcs6bVnWHMuyCtyxTqKeH0nZwvR09PdvAhQCWgFNgV63VrAsyxdYCKQGngUaAPWA/kmYyxR/YAYwIo71/gdkve0xIYlzmXLf4+Fl58a9dCb2ufC12TjOYVlWfWAQ0BMoAawBFluW9aTJXAbtIfZ5UMRsHKdKDewAPr3H6+2B1tGvBwIngeWWZXnqzeHjOh4AS4h9vlRzQi4TKuGY0LIc8DKOO5cssywr9W3rJO75Ydu20x5AO+Dwbc+rApFAttuWNQDCgLTOzObEY9AIuHCP12ygtumMrnA8vPHcuOP9/wG0NJ3D0HvfCIy4Y9k+oLfpbAaORVdgu+kcrvC48/MRsIATwBe3LUsOXACamM7r7OMRvWw8MMd0NkPH4/HoY/JcUp0fzh7DlA44d9vzZ4Ddtm3/fduypTjeVClnBnMhw6K7nzZZltXUsixvHWemcwO+sCzrrGVZ26O7MD2+OzL6PZYClt3x0jKgvPMTuYR80V0KRyzLmmZZVm7TgVxELiALt50rtm3fAFbhvecKwPPRXVT7LcsaY1lWJtOBnCRd9L+3aoxEPz8SdPPdhLAsKw/wGdDmtsVZgFO3r2fb9nnLsm5Gv+ZtOgErgOvAizi6nzLiJV0xd/D2c2MwsBU4j+Mu3b1xfAB8ZDKUE2QEfLnj/z76uTf8v99pI/AesB/IDHQE1lmWVci27bNGk5l363y427mS08lZXMViHEMdjuL4vOgB/GRZVqnoYsEjWZZlAQOAX2zb3h29ONHPjwduvbjHIMQ7H6Xv2CYbjn7VGbZtj73jW95tqnHrHstdSkKOxf3Ytv21bdvrbdvebtt2fxxjWNol3TtIXIl9PHDjc+NuHuT42LY90LbtVbZt74z+mWkKfGhZ1mNm34XT3Pl/7Lb/7w/Dtu3Ftm2H2ra9y7btH4Hq0S+9bzKXi9G5Es227em2bS+0bXu3bdvzcQxtyM+/542nGgYUBd66y2uJdn4kpIVpGDAtjnX+uPVFdLH0M7AeaHzHeieBsrcvsCzrUcCP/1aFruiBjkUCbADSWpaV2bZtbzse7n5u3M3DHJ9bV7rkBTy5ZeEMjrFrd7YmZcJ9/98TjW3bVy3L2gXkM53FBdy6WjALjrEqt+hciWbb9gnLso7iweeLZVlDgVo4xi79ddtLiX5+PHDBZNv2GRwfanGyLCs7jmJpC/A/27aj7lhlPfCVZVlZbdu+9Yaq4Lhh3pYHzeZsD3IsEqgEjkHOF5JwH4kmkY+HW58bd/OQx6dE9L8n7ruWm7Nt+6ZlWVtwXPUy+7aXXgbmmknlOizLSg4UxHHloLc7guOX4svANogZA1cJ+MJgLpcR3SKdAw/83IjuhhsK1AGet237yB2rJPr5kWRjmKJbllYCx4C2wOOO9we2bd+q/JYBe4FJlmW1AzIA/YAxtm1fSqpsJkRfEp0BeBLwtSyrePRLB23bvmJZVk0clfB6HGOYXsBxWfVoT+x7jut44EXnxp0sy3oGx6WyPwMXcVwOOxCYZ9v2MZPZnGQAjv/3zfzbMv0kMNJoKgMsy+oHzMfxOZoJxximtHjudCOxWJb1CI5W1VtyRX9WnLNt+5hlWYOADpZlHQAOAB2Aa8D3zs7qDPc7HtGPrkAojgLpKRzT+Jwh9h8fniIEaAi8Bly2LOtWq/RF27av27ZtJ/r5kYSX+DXC0U/4n8cd6z0JLIh+E2dxVIzJTV+imATHY/w9jsfz0a+/iqMKvgxcBXYBLYBkprObOB7edG7c5diUxNEFdwFH8fwbjg/CVKazOfEYNMPRPXmrRfE505kMHYdpwN/ATeA4jl+GAaZzOfH9P3+Pz4nx0a9b0T8bJ3C0xq8CCpvObeJ4AClxXEl8Ovp8ORq9PIfp3El0LO5aXwCNblsnUc8PK/qbioiIiMg9eOscPyIiIiLxpoJJREREJA4qmERERETioIJJREREJA4qmERERETioIJJREREJA4qmERERETioIJJREREJA4qmERERETioIJJREREJA4qmERERETioIJJREREJA7/B+2FefvlJ090AAAAAElFTkSuQmCC\n", "text/plain": [ "Graphics object consisting of 42 graphics primitives" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Plot the series coefficients of [y^N]F(x,y)/asm(n=N) compared to the limit curve\n", "N = 100\n", "SN = ((x+1+2/x)^N).expand()\n", "mm = m[0]*N\n", "\n", "pt = plot([])\n", "for k in range(-20,21):\n", " pt += point([k,SN.coefficient(x,mm+k)/asm.subs(n=N)], size=20, color='red')\n", "\n", "var('L')\n", "pt += plot(norm.subs(ep=m[0]+L/N,n=N),L,-20,20, color='black')\n", "pt" ] } ], "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 }