{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercises: An Invitation to Analytic Combinatorics in Several Variables\n", "\n", "Created by Stephen Melczer\n", "\n", "A quick Sage tutorial can be found [here](https://melczer.ca/files/SageIntro.html) (try an interactive version in your browser [here](https://mybinder.org/v2/git/https%3A%2F%2Fgit.uwaterloo.ca%2Fsmelczer%2Fintro-to-sage/HEAD?filepath=A%20Brief%20Introduction%20to%20Sage.ipynb)).\n", "\n", "See [https://melczer.ca/textbook](https://melczer.ca/textbook) for further Sage notebooks solving problems in analytic combinatorics in several varibles. \n", "\n", "In particular, [this notebook](https://melczer.ca/files/TextbookCode/Chapter5/Example5-SmoothASM.ipynb) (also available as a [static HTML page](https://melczer.ca/files/TextbookCode/Chapter5/Example5-SmoothASM.html)) gives an algorithm to compute asymptotic terms. Don't use it to solve these exercises, but use it to further check your answers! (Or to compute asymptotics for other problems!)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Helper function to compute the Hessian matrix M from the function H(vars)\n", "# in the direction R at the point CP, specified by a list of substitutions.\n", "# Copied from melczer.ca/textbook/\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.subs(CP)\n", "\n", "\n", "# Helper function to compute leading asymptotics of the R-diagonal of G(vars)/H(vars)\n", "# determined by the Main Asymptotic Theorem of Smooth ACSV at the point CP, specified\n", "# by a list of substitutions. We take det(M) as an input that can be computed by the\n", "# above function.\n", "var('n')\n", "def leadingASM(G,H,detM,R,vars,CP):\n", " dd = len(R)\n", " lcoeff = -G/vars[-1]/H.diff(vars[-1])\n", " exp = 1/mul([vars[k]^R[k] for k in range(dd)])^n\n", " \n", " ASM = exp * (2*pi*n*R[-1])^((1-dd)/2) / sqrt(detM) * lcoeff\n", " return ASM.subs(CP)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Question 1: Delannoy Numbers\n", "The *Delannoy number* $d_{a,b}$ is the number of paths from the origin $(0,0)$ to the point $(a,b)$ using only the steps $\\textsf{N}=(0,1)$, $\\textsf{E} = (1,0)$, and $\\textsf{NE}=(1,1)$.\n", "\n", "**(a)** Prove the recurrence \n", "$$ d_{a,b} = \n", "\\begin{cases} \n", "1 &: \\text{ if $a=0$ or $b=0$} \\\\\n", "d_{a-1,b} + d_{a,b-1} + d_{a-1,b-1} &:\\text{ otherwise}\n", "\\end{cases}\n", "$$\n", "Conclude that \n", "$$ D(x,y) = \\sum_{a,b\\geq0}d_{a,b}x^ay^b = \\frac{1}{1-x-y-xy}. $$\n", "\n", "**(b)** Use the Main Theorem of Smooth ACSV to find asymptotics of $d_{n,n}$ as the $(1,1)$-diagonal of $D(x,y)$. What are the critical points in the $(1,1)$ direction? Which are minimal?\n", "\n", "**(c)** Use the Main Theorem of Smooth ACSV to find asymptotics of the $(r,s)$-diagonal of $D(x,y)$ for any $r,s>0$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Answer to 1a\n", "\n", "If $a=b=0$ then we say $d_{0,0}=1$ by convention (there is one way to do nothing). Otherwise, if $a=0$ then there is exactly one walk to $(a,b)$, consisting of $\\textsf{N}$ repeated $b$ times. Similarly, if $b=0$ then there is exactly one walk to $(a,b)$, consisting of $\\textsf{E}$ repeated $a$ times. Thus, the initial conditions hold.\n", "\n", "If $a,b>0$ then a walk ending at $(a,b)$ is \n", "- a walk ending at $(a-1,b)$ followed by an $\\textsf{E}$ step, or\n", "- a walk ending at $(a,b-1)$ followed by an $\\textsf{N}$ step, or\n", "- a walk ending at $(a-1,b-1)$ followed by a $\\textsf{NE}$ step.\n", "\n", "The possibilities are counted by $d_{a-1,b}, d_{a,b-1}, $ and $d_{a-1,b-1}$, respectively, giving the stated recurrence. The recurrence then implies\n", "$$ (1-x-y-xy) \\sum_{a,b\\geq0}d_{a,b}x^ay^b = 1, $$\n", "proving $D(x,y)=1/(1-x-y-xy)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Answers to 1b and 1c\n", "\n", "Part (b) is a special case of part (c), so we simply solve part (c)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "We consider the (r,s)-diagonal of\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|F(x,y)|\\phantom{\\verb!x!}\\verb|=| -\\frac{1}{x y + x + y - 1}\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|F(x,y)|\\phantom{\\verb!x!}\\verb|=| -\\frac{1}{x y + x + y - 1}$$" ], "text/plain": [ "'F(x,y) = ' -1/(x*y + x + y - 1)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Introduce variables and function\n", "var('x y r s t n')\n", "H = 1 - x - y - x*y\n", "print('We consider the (r,s)-diagonal of')\n", "show('F(x,y) = ', 1/H)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The critical points in direction (r,s) are\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = -\\frac{s + \\sqrt{r^{2} + s^{2}}}{r}, y = -\\frac{r + \\sqrt{r^{2} + s^{2}}}{s}\\right]\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = -\\frac{s + \\sqrt{r^{2} + s^{2}}}{r}, y = -\\frac{r + \\sqrt{r^{2} + s^{2}}}{s}\\right]$$" ], "text/plain": [ "[x == -(s + sqrt(r^2 + s^2))/r, y == -(r + sqrt(r^2 + s^2))/s]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = -\\frac{s - \\sqrt{r^{2} + s^{2}}}{r}, y = -\\frac{r - \\sqrt{r^{2} + s^{2}}}{s}\\right]\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = -\\frac{s - \\sqrt{r^{2} + s^{2}}}{r}, y = -\\frac{r - \\sqrt{r^{2} + s^{2}}}{s}\\right]$$" ], "text/plain": [ "[x == -(s - sqrt(r^2 + s^2))/r, y == -(r - sqrt(r^2 + s^2))/s]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Solve the critical point equations in direction (r,s)\n", "rts = solve([H,s*x*diff(H,x) - r*y*diff(H,y)],[x,y])\n", "\n", "# Print results\n", "print(\"The critical points in direction (r,s) are\")\n", "for k in rts:\n", " show(k)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The unique minimal critical point is\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = -\\frac{s - \\sqrt{r^{2} + s^{2}}}{r}, y = -\\frac{r - \\sqrt{r^{2} + s^{2}}}{s}\\right]\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = -\\frac{s - \\sqrt{r^{2} + s^{2}}}{r}, y = -\\frac{r - \\sqrt{r^{2} + s^{2}}}{s}\\right]$$" ], "text/plain": [ "[x == -(s - sqrt(r^2 + s^2))/r, y == -(r - sqrt(r^2 + s^2))/s]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# When r,s > 0 only one of these points has positive coordinates, which lie in (0,1)\n", "# It is also coordinate-wise smaller than the other critical point\n", "CP = [x == -(s - sqrt(r^2 + s^2))/r, y == -(r - sqrt(r^2 + s^2))/s]\n", "\n", "# It is minimal as |x| = |1-y|/|1+y| decreases as y increases in (0,1) when H(x,y) = 0\n", "print('The unique minimal critical point is')\n", "show(CP)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The matrix M is the 1 x 1 matrix (i.e., constant)\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{r}\n", "\\frac{{\\left(r^{2} + s^{2} + \\sqrt{r^{2} + s^{2}} r - \\sqrt{r^{2} + s^{2}} s\\right)} r}{{\\left(r - s + \\sqrt{r^{2} + s^{2}}\\right)} s^{2}}\n", "\\end{array}\\right)\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{r}\n", "\\frac{{\\left(r^{2} + s^{2} + \\sqrt{r^{2} + s^{2}} r - \\sqrt{r^{2} + s^{2}} s\\right)} r}{{\\left(r - s + \\sqrt{r^{2} + s^{2}}\\right)} s^{2}}\n", "\\end{array}\\right)$$" ], "text/plain": [ "[(r^2 + s^2 + sqrt(r^2 + s^2)*r - sqrt(r^2 + s^2)*s)*r/((r - s + sqrt(r^2 + s^2))*s^2)]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Asymptotics depends on the 1x1 Hessian matrix M\n", "M = getHes(H,[r,s],[x,y],CP)\n", "detM = M.determinant().factor()\n", "print(\"The matrix M is the 1 x 1 matrix (i.e., constant)\")\n", "show(M.factor())" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "In this case, M is also the second derivative of the function\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\log\\left(-\\frac{s e^{\\left(i \\, t\\right)} - \\sqrt{r^{2} + s^{2}} e^{\\left(i \\, t\\right)} + r}{s e^{\\left(i \\, t\\right)} - \\sqrt{r^{2} + s^{2}} e^{\\left(i \\, t\\right)} - r}\\right)\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\log\\left(-\\frac{s e^{\\left(i \\, t\\right)} - \\sqrt{r^{2} + s^{2}} e^{\\left(i \\, t\\right)} + r}{s e^{\\left(i \\, t\\right)} - \\sqrt{r^{2} + s^{2}} e^{\\left(i \\, t\\right)} - r}\\right)$$" ], "text/plain": [ "log(-(s*e^(I*t) - sqrt(r^2 + s^2)*e^(I*t) + r)/(s*e^(I*t) - sqrt(r^2 + s^2)*e^(I*t) - r))" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "These two expressions equal: True\n" ] } ], "source": [ "# OPTIONAL NOTE: We could have derived M by parametrizing y = g(x) and taking the\n", "# \"Hessian\" (which is the second derivative in one dimension) of ϕ(t)\n", "g = solve(H,y)[0].rhs()\n", "phi = log(g.subs(x=CP[0].rhs()*exp(I*t)).simplify_full())\n", "print('In this case, M is also the second derivative of the function')\n", "show(phi)\n", "eq = bool(1==(detM/phi.diff(t,t).subs(t==0)).simplify_full())\n", "print('These two expressions equal:', eq)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The dominant asymptotic behaviour in the direction (r,s) is\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\sqrt{2} s}{2 \\, \\sqrt{\\pi n s} \\left(\\left(-\\frac{s - \\sqrt{r^{2} + s^{2}}}{r}\\right)^{r} \\left(-\\frac{r - \\sqrt{r^{2} + s^{2}}}{s}\\right)^{s}\\right)^{n} {\\left(r - \\sqrt{r^{2} + s^{2}}\\right)} {\\left(\\frac{s - \\sqrt{r^{2} + s^{2}}}{r} - 1\\right)} \\sqrt{\\frac{{\\left(r^{2} + s^{2} + \\sqrt{r^{2} + s^{2}} r - \\sqrt{r^{2} + s^{2}} s\\right)} r}{{\\left(r - s + \\sqrt{r^{2} + s^{2}}\\right)} s^{2}}}}\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\sqrt{2} s}{2 \\, \\sqrt{\\pi n s} \\left(\\left(-\\frac{s - \\sqrt{r^{2} + s^{2}}}{r}\\right)^{r} \\left(-\\frac{r - \\sqrt{r^{2} + s^{2}}}{s}\\right)^{s}\\right)^{n} {\\left(r - \\sqrt{r^{2} + s^{2}}\\right)} {\\left(\\frac{s - \\sqrt{r^{2} + s^{2}}}{r} - 1\\right)} \\sqrt{\\frac{{\\left(r^{2} + s^{2} + \\sqrt{r^{2} + s^{2}} r - \\sqrt{r^{2} + s^{2}} s\\right)} r}{{\\left(r - s + \\sqrt{r^{2} + s^{2}}\\right)} s^{2}}}}$$" ], "text/plain": [ "1/2*sqrt(2)*s/(sqrt(pi*n*s)*((-(s - sqrt(r^2 + s^2))/r)^r*(-(r - sqrt(r^2 + s^2))/s)^s)^n*(r - sqrt(r^2 + s^2))*((s - sqrt(r^2 + s^2))/r - 1)*sqrt((r^2 + s^2 + sqrt(r^2 + s^2)*r - sqrt(r^2 + s^2)*s)*r/((r - s + sqrt(r^2 + s^2))*s^2)))" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "For instance, when r = s = 1 we have\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{2^{\\frac{3}{4}}}{4 \\, \\sqrt{\\pi n} {\\left({\\left(\\sqrt{2} - 1\\right)}^{2}\\right)}^{n} {\\left(\\sqrt{2} - 1\\right)}}\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{2^{\\frac{3}{4}}}{4 \\, \\sqrt{\\pi n} {\\left({\\left(\\sqrt{2} - 1\\right)}^{2}\\right)}^{n} {\\left(\\sqrt{2} - 1\\right)}}$$" ], "text/plain": [ "1/4*2^(3/4)/(sqrt(pi*n)*((sqrt(2) - 1)^2)^n*(sqrt(2) - 1))" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# The partial derivative of H with respect to y doesn't vanish where H is zero,\n", "# thus we get the following dominant asymptotics\n", "detM = getHes(H,[r,s],[x,y],CP).determinant().factor()\n", "ASM = leadingASM(1,H,detM,[r,s],[x,y],CP)\n", "print(\"The dominant asymptotic behaviour in the direction (r,s) is\")\n", "show(ASM)\n", "print(\"For instance, when r = s = 1 we have\")\n", "show(ASM.subs(r==1,s==1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can numerically check this asymptotic formula with actual series terms." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Checking the 50th term on the (r,s)-diagonal with our asymptotic formula,\n", "In direction (1,1) we get ratio 0.997656071379095\n", "In direction (2,1) we get ratio 0.998120119323750\n", "In direction (2,3) we get ratio 0.998996193385674\n" ] } ], "source": [ "# First, define the ring of formal power series (more efficient for computations)\n", "P. = QQ[['X,Y']] \n", "\n", "# Function to check ratio of actual coefficients to asymptotic approximation\n", "# for x^(N*R) * y^(N*S) where R, S, and N are positive integers\n", "def asmCheck(R,S,N):\n", " N2 = (R+S)*N\n", " ser = 1/(1-X-Y-X*Y + O(X^(N2+1)))\n", " cfs = ser.coefficients()\n", " return (cfs[X^(R*N)*Y^(S*N)]/ASM.subs(r==R,s==S,n==N)).n()\n", "\n", "print(\"Checking the 50th term on the (r,s)-diagonal with our asymptotic formula,\")\n", "print(\"In direction (1,1) we get ratio\", asmCheck(1,1,50))\n", "print(\"In direction (2,1) we get ratio\", asmCheck(2,1,50))\n", "print(\"In direction (2,3) we get ratio\", asmCheck(2,3,50))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Question 2: Apéry Asymptotics\n", "\n", "Recall from lecture that a key step in Apéry's proof of the irrationality of $\\zeta(3)$ is determining the exponential growth of the sequence that can be encoded as the main diagonal\n", "of \n", "$$ F(x,y,z,t) = \\frac{1}{1 - t(1+x)(1+y)(1+z)(1+y+z+yz+xyz)}. $$\n", "Use the Main Theorem of Smooth ACSV to find dominant asymptotics of this sequence." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Answer to 2" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The critical point equations are\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[-{\\left(x y z + y z + y + z + 1\\right)} t {\\left(x + 1\\right)} {\\left(y + 1\\right)} {\\left(z + 1\\right)} + 1, -{\\left(t {\\left(x + 1\\right)} {\\left(y + 1\\right)} y {\\left(z + 1\\right)} z + {\\left(x y z + y z + y + z + 1\\right)} t {\\left(y + 1\\right)} {\\left(z + 1\\right)}\\right)} x + {\\left({\\left(x z + z + 1\\right)} t {\\left(x + 1\\right)} {\\left(y + 1\\right)} {\\left(z + 1\\right)} + {\\left(x y z + y z + y + z + 1\\right)} t {\\left(x + 1\\right)} {\\left(z + 1\\right)}\\right)} y, -{\\left(t {\\left(x + 1\\right)} {\\left(y + 1\\right)} y {\\left(z + 1\\right)} z + {\\left(x y z + y z + y + z + 1\\right)} t {\\left(y + 1\\right)} {\\left(z + 1\\right)}\\right)} x + {\\left({\\left(x y + y + 1\\right)} t {\\left(x + 1\\right)} {\\left(y + 1\\right)} {\\left(z + 1\\right)} + {\\left(x y z + y z + y + z + 1\\right)} t {\\left(x + 1\\right)} {\\left(y + 1\\right)}\\right)} z, {\\left(x y z + y z + y + z + 1\\right)} t {\\left(x + 1\\right)} {\\left(y + 1\\right)} {\\left(z + 1\\right)} - {\\left(t {\\left(x + 1\\right)} {\\left(y + 1\\right)} y {\\left(z + 1\\right)} z + {\\left(x y z + y z + y + z + 1\\right)} t {\\left(y + 1\\right)} {\\left(z + 1\\right)}\\right)} x\\right]\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[-{\\left(x y z + y z + y + z + 1\\right)} t {\\left(x + 1\\right)} {\\left(y + 1\\right)} {\\left(z + 1\\right)} + 1, -{\\left(t {\\left(x + 1\\right)} {\\left(y + 1\\right)} y {\\left(z + 1\\right)} z + {\\left(x y z + y z + y + z + 1\\right)} t {\\left(y + 1\\right)} {\\left(z + 1\\right)}\\right)} x + {\\left({\\left(x z + z + 1\\right)} t {\\left(x + 1\\right)} {\\left(y + 1\\right)} {\\left(z + 1\\right)} + {\\left(x y z + y z + y + z + 1\\right)} t {\\left(x + 1\\right)} {\\left(z + 1\\right)}\\right)} y, -{\\left(t {\\left(x + 1\\right)} {\\left(y + 1\\right)} y {\\left(z + 1\\right)} z + {\\left(x y z + y z + y + z + 1\\right)} t {\\left(y + 1\\right)} {\\left(z + 1\\right)}\\right)} x + {\\left({\\left(x y + y + 1\\right)} t {\\left(x + 1\\right)} {\\left(y + 1\\right)} {\\left(z + 1\\right)} + {\\left(x y z + y z + y + z + 1\\right)} t {\\left(x + 1\\right)} {\\left(y + 1\\right)}\\right)} z, {\\left(x y z + y z + y + z + 1\\right)} t {\\left(x + 1\\right)} {\\left(y + 1\\right)} {\\left(z + 1\\right)} - {\\left(t {\\left(x + 1\\right)} {\\left(y + 1\\right)} y {\\left(z + 1\\right)} z + {\\left(x y z + y z + y + z + 1\\right)} t {\\left(y + 1\\right)} {\\left(z + 1\\right)}\\right)} x\\right]$$" ], "text/plain": [ "[-(x*y*z + y*z + y + z + 1)*t*(x + 1)*(y + 1)*(z + 1) + 1,\n", " -(t*(x + 1)*(y + 1)*y*(z + 1)*z + (x*y*z + y*z + y + z + 1)*t*(y + 1)*(z + 1))*x + ((x*z + z + 1)*t*(x + 1)*(y + 1)*(z + 1) + (x*y*z + y*z + y + z + 1)*t*(x + 1)*(z + 1))*y,\n", " -(t*(x + 1)*(y + 1)*y*(z + 1)*z + (x*y*z + y*z + y + z + 1)*t*(y + 1)*(z + 1))*x + ((x*y + y + 1)*t*(x + 1)*(y + 1)*(z + 1) + (x*y*z + y*z + y + z + 1)*t*(x + 1)*(y + 1))*z,\n", " (x*y*z + y*z + y + z + 1)*t*(x + 1)*(y + 1)*(z + 1) - (t*(x + 1)*(y + 1)*y*(z + 1)*z + (x*y*z + y*z + y + z + 1)*t*(y + 1)*(z + 1))*x]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Introduce variables and function\n", "var('x,y,z,t')\n", "H = 1 - t*(1+x)*(1+y)*(1+z)*(1+y+z+y*z+x*y*z)\n", "R = [1,1,1,1]\n", "vari = [x,y,z,t]\n", "\n", "# Get the smooth critical point equations\n", "d = len(vari)\n", "criteqs = [H] + [R[j]*vari[0]*diff(H,vari[0]) - R[0]*vari[j]*diff(H,vari[j]) for j in range(1,d)]\n", "print(\"The critical point equations are\")\n", "show(criteqs)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The critical points are\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[t = -58 \\, \\sqrt{2} - 82, x = -\\sqrt{2} + 1, y = -\\frac{1}{2} \\, \\sqrt{2}, z = -\\frac{1}{2} \\, \\sqrt{2}\\right]\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[t = -58 \\, \\sqrt{2} - 82, x = -\\sqrt{2} + 1, y = -\\frac{1}{2} \\, \\sqrt{2}, z = -\\frac{1}{2} \\, \\sqrt{2}\\right]$$" ], "text/plain": [ "[t == -58*sqrt(2) - 82,\n", " x == -sqrt(2) + 1,\n", " y == -1/2*sqrt(2),\n", " z == -1/2*sqrt(2)]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[t = 58 \\, \\sqrt{2} - 82, x = \\sqrt{2} + 1, y = \\frac{1}{2} \\, \\sqrt{2}, z = \\frac{1}{2} \\, \\sqrt{2}\\right]\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[t = 58 \\, \\sqrt{2} - 82, x = \\sqrt{2} + 1, y = \\frac{1}{2} \\, \\sqrt{2}, z = \\frac{1}{2} \\, \\sqrt{2}\\right]$$" ], "text/plain": [ "[t == 58*sqrt(2) - 82, x == sqrt(2) + 1, y == 1/2*sqrt(2), z == 1/2*sqrt(2)]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Of these, the minimal critical point is\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[t = 58 \\, \\sqrt{2} - 82, x = \\sqrt{2} + 1, y = \\frac{1}{2} \\, \\sqrt{2}, z = \\frac{1}{2} \\, \\sqrt{2}\\right]\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[t = 58 \\, \\sqrt{2} - 82, x = \\sqrt{2} + 1, y = \\frac{1}{2} \\, \\sqrt{2}, z = \\frac{1}{2} \\, \\sqrt{2}\\right]$$" ], "text/plain": [ "[t == 58*sqrt(2) - 82, x == sqrt(2) + 1, y == 1/2*sqrt(2), z == 1/2*sqrt(2)]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# There are two critical points, the one with positive coordinates is minimal\n", "[CP1,CP2] = solve(criteqs,vari)\n", "\n", "if t.subs(CP1) > 0:\n", " CP = CP1\n", "else:\n", " CP = CP2\n", "\n", "print(\"The critical points are\")\n", "show(CP1); show(CP2)\n", "print(\"Of these, the minimal critical point is\")\n", "show(CP)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dominant asymptotic behaviour of the main diagonal is\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\sqrt{2}}{8 \\, \\left(\\pi n\\right)^{\\frac{3}{2}} \\left({\\left(29 \\, \\sqrt{2} - 41\\right)} {\\left(\\sqrt{2} + 1\\right)}\\right)^{n} {\\left(29 \\, \\sqrt{2} - 41\\right)} {\\left(3 \\, \\sqrt{2} + 4\\right)} \\sqrt{\\frac{408 \\, \\sqrt{2} + 577}{{\\left(3 \\, \\sqrt{2} + 4\\right)}^{3}}}}\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\sqrt{2}}{8 \\, \\left(\\pi n\\right)^{\\frac{3}{2}} \\left({\\left(29 \\, \\sqrt{2} - 41\\right)} {\\left(\\sqrt{2} + 1\\right)}\\right)^{n} {\\left(29 \\, \\sqrt{2} - 41\\right)} {\\left(3 \\, \\sqrt{2} + 4\\right)} \\sqrt{\\frac{408 \\, \\sqrt{2} + 577}{{\\left(3 \\, \\sqrt{2} + 4\\right)}^{3}}}}$$" ], "text/plain": [ "1/8*sqrt(2)/((pi*n)^(3/2)*((29*sqrt(2) - 41)*(sqrt(2) + 1))^n*(29*sqrt(2) - 41)*(3*sqrt(2) + 4)*sqrt((408*sqrt(2) + 577)/(3*sqrt(2) + 4)^3))" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Note: it's possible to algebraically simplify this formula somewhat\n" ] } ], "source": [ "print(\"Dominant asymptotic behaviour of the main diagonal is\")\n", "detM = getHes(H,[1,1,1,1],[x,y,z,t],CP).determinant().factor()\n", "ASM = leadingASM(1,H,detM,[1,1,1,1],[x,y,z,t],CP)\n", "show(ASM)\n", "print(\"Note: it's possible to algebraically simplify this formula somewhat\")" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "When n = 40 the ratio of our formula and actual coefficient is 0.989591444404954\n" ] } ], "source": [ "# Numerically check the formula (this can take a decent amount of time for large N)\n", "P. = QQ['X,Y,Z'] \n", "N = 40\n", "ser = ((1+X)*(1+Y)*(1+Z)*(1+Y+Z+Y*Z+X*Y*Z))^N\n", "rat = (ser[X^N*Y^N*Z^N]/ASM.subs(n==N)).n()\n", "print(\"When n = {} the ratio of our formula and actual coefficient is\".format(N),rat)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Question 3: Pathological Directions\n", "**(a)** Find asymptotics of the $(r,s)$-diagonal of $F(x,y) = \\frac{1}{1-x-xy}$\n", "for any $0\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = \\frac{r - s}{r}, y = \\frac{s}{r - s}\\right]\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[x = \\frac{r - s}{r}, y = \\frac{s}{r - s}\\right]$$" ], "text/plain": [ "[x == (r - s)/r, y == s/(r - s)]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Solve the critical point equations in direction (r,s) and print results\n", "[CP] = solve([H,s*x*diff(H,x) - r*y*diff(H,y)],[x,y])\n", "print(\"The critical point in direction (r,s) is\")\n", "show(CP)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "If r > s then this point is minimal and the dominant asymptotic behaviour is\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\sqrt{2} r}{2 \\, \\sqrt{\\pi n s} \\left(\\left(\\frac{r - s}{r}\\right)^{r} \\left(\\frac{s}{r - s}\\right)^{s}\\right)^{n} s \\sqrt{\\frac{{\\left(r - s\\right)} r}{s^{2}}}}\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\sqrt{2} r}{2 \\, \\sqrt{\\pi n s} \\left(\\left(\\frac{r - s}{r}\\right)^{r} \\left(\\frac{s}{r - s}\\right)^{s}\\right)^{n} s \\sqrt{\\frac{{\\left(r - s\\right)} r}{s^{2}}}}$$" ], "text/plain": [ "1/2*sqrt(2)*r/(sqrt(pi*n*s)*(((r - s)/r)^r*(s/(r - s))^s)^n*s*sqrt((r - s)*r/s^2))" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Compute asymptotic contribution of this point\n", "detM = getHes(H,[r,s],[x,y],CP).determinant().factor()\n", "ASM = leadingASM(1,H,detM,[r,s],[x,y],CP)\n", "print('If r > s then this point is minimal and the dominant asymptotic behaviour is')\n", "show(ASM)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If $r < s$ then this point is **not** minimal, and there are **no minimal critical points**. Note that \n", "$$F(x,y) = 1/(1-x-xy) = \\sum_{n \\geq 0}x^n(1+y)^n$$ \n", "so $[x^ay^b]F(x,y)=0$ if $a = QQ[['X,Y']] \n", "def asmCheck(R,S,N):\n", " N2 = (R+S)*N\n", " ser = 1/(1-X-X*Y + O(X^(N2+1)))\n", " cfs = ser.coefficients()\n", " return (cfs[X^(R*N)*Y^(S*N)]/ASM.subs(r==R,s==S,n==N)).n()\n", "\n", "print(\"Checking the 50th term on the (r,s)-diagonal with our asymptotic formula,\")\n", "print(\"In direction (2,1) we get ratio\", asmCheck(2,1,50))\n", "print(\"In direction (3,2) we get ratio\", asmCheck(3,2,50))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Question 4: A Composition Limit Theorem\n", "\n", "An *integer composition* of size $n\\in\\mathbb{N}$ is an ordered tuple of positive integers (of any length) that sum to $n$. A composition can be viewed as an integer partition where the order of the summands matters. Let $c_{k,n}$ denote the number of compositions of size $n$ that contain $k$ ones. \n", "\n", "**(a)** If you know the symbolic method, species theory, or similar enumerative constructions, prove that \n", "$$ C(u,x) = \\sum_{n,k\\geq0}c_{k,n}u^kx^n = \\frac{1-x}{1-2x-(u-1)x(1-x)}. $$\n", "\n", "**(b)** Prove that the distribution for the number of ones in a composition of size $n$ satisfies a local central limit theorem as $n\\rightarrow\\infty$. More precisely, find a constant $t>0$ and normal density $\\nu_n(s)$ such that \n", "$$ \\sup_{s \\in \\mathbb{Z}} |t^nc_{s,n} - \\nu_n(s)| \\rightarrow 0 $$\n", "as $n\\rightarrow\\infty$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Answer to 4a\n", "\n", "A composition can be viewed as a sequence of positive integers. Separating and marking the integer 1, the combinatorial class $\\mathcal{C}$ has the specification\n", "$$ \\mathcal{C} = \\text{SEQ}\\left(Z\\times U \\; + \\; Z^2 \\times \\text{SEQ}(Z) \\right) $$\n", "giving the generating function\n", "$$ C(u,x) = \\frac{1}{1-xu - \\frac{x^2}{1-x}} , $$\n", "which simplies to the stated expression. (See *Analytic Combinatorics* by Flajolet and Sedgewick for information on the symbolic method)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "var('u x')\n", "(1/(1-x*u-x^2/(1-x)) - (1-x)/(1-2*x-(u-1)*x*(1-x))).simplify_full()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Answer to 4b" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "# Introduce variables and function\n", "var('u x s')\n", "G = 1 - x\n", "H = 1 - 2*x - (u-1)*x*(1-x)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "H(1,x) = 0 implies [u == 1, x == (1/2)]\n" ] } ], "source": [ "# We want a critical point where u = 1\n", "# There is only one such point on V\n", "CP = [u==1] + solve(H.subs(u==1),x)\n", "print(\"H(1,x) = 0 implies\", CP)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAGFCAYAAAAPa6wiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAuNElEQVR4nO3dd5xU1f3/8dehuHTF/k1sUQQUCxZQmiJEjQUFY2JBY/km6lcSvyYxUdREf0mMxnyTqImmmdjBgiIKiqgEVBBRsAY1RsBeopGi0rm/P85sdlkXd3d2Zu+dmdfz8ZgHODN35vPZu9F3zjn33JAkCZIkSVq/VmkXIEmSlHUGJkmSpAYYmCRJkhpgYJIkSWqAgUmSJKkBBiZJkqQGGJgkSZIaYGCSJElqgIFJkiSpAQYmSZKkBjQ5MIUQ9gsh3BtCeDuEkIQQhtd5PYQQLs69viyEMC2E0KsJnx9CCF1CCKGptUmSJBVDmzyO6Qg8C1wH3FnP6z8EvgecDPwDuBB4MITQI0mSpY34/M7AYmBDYEkTa/PGeJIkaX3yHowJzbn5bgghAUYkSXJ37p8D8DZwRZIkv8g9VwW8B5ybJMkfG/GZXcgFpiRJDEySJKlQ8g5MhV7D9CVgS2BK9RNJkqwApgP96zsghFCVm4LrkgtLnQtckyRJUrMUOjBtmfvzvTrPv1frtbpGE0eUqh9vFrgmSZKkZinWVXJ1p8ZCPc9Vu5S4Xqn6sVWRapIkScpLoQPTu7k/644mbc5nR52AOGWXJMmS6gfQmIXhzfLqq8X+BkmSVE4KHZgWEEPTgdVPhBA2APYHZhb4u/IybhzsuCO88UbalUiSpFKRzz5MnUIIvUMIvXNPfSn3z9sk8ZK7K4DzQwgjQgi7ANcDnwJjClRzsxx0EFRVwdixaVciSZJKRZO3FQghDAb+Vs9LNyRJcnJua4GLgNOBrsATwKgkSV5o5OcXfVuBY4+FefPgueea+OmSJKmUpbMPUzG0RGCaOBGGDYNnn4XddmtyiZIkqTRlZh+mknDwwbDJJnDzzWlXIkmSSkFFBqa2beGYY2DMGFizJu1qJElS1lVkYAI44QR46y2YPj3tSiRJUtZVbGDad1/Yfnun5SRJUsMqNjCFEEeZxo2DZcvSrkaSJGVZxQYmgJEjYenSeNWcJEnS+lR0YOreHfr2dVpOkiR9vooOTBCn5e67D95/P+1KJElSVlV8YDruuLie6ZZb0q5EkiRlVcUHpk03hSOPhL/+FTK26bkkScqIzASmEMKoEMI8YHZLf/cpp8ALL8CcOS39zZIkqRRU5L3k6lq9GrbdNo40XXNNPp8gSZJKgPeSa442beCkk2DsWPdkkiRJn2Vgyjn5ZFi0CO6+O+VCJElS5jglV8ugQdC+PUyZ0pxPkSRJGeWUXCGccgo89BC8/nralUiSpCwxMNXyta9Bhw5www1pVyJJkrLEwFRL587w9a/DddfB2rVpVyNJkrLCwFTHqafCggUwdWralUiSpKxw0XfdD0hgl11gp51g3LjmfpokScoQF30XSghwxhkwYQK8807a1UiSpCwwMNXjxBOhbdt4fzlJkiQDUz022giOOw7+9CdYsybtaiRJUtoMTOtx+ulxP6bJk9OuRJIkpc1F3+v7oAT22gu22gruuadQnypJklLkou9Cq178PWmSO39LklTpDEyf47jjoGNHuPbatCuRJElpMjB9js6d4YQTYmBatSrtaiRJUloyE5hCCKNCCPOA2WnXUtvpp8f9mCZMSLsSSZKUFhd9N8J++8U1TdOnF+PTJUlSC3HRdzGddRY88gg880zalUiSpDQYmBph+HDYemu46qq0K5EkSWkwMDVCmzYwahSMGQP/+lfa1UiSpJZmYGqkb34zrmP685/TrkSSJLU0A1MjbbJJ3GLgmmvcYkCSpEpjYGqCs86Ct96Cu+5KuxJJktSS3FagiYYMgRUrYMaMYn+TJEkqMLcVaClnnQUzZ8JTT6VdiSRJaimOMDXRmjWw446w777xqjlJklQyHGFqKa1bw3e/C7ffDq+9lnY1kiSpJRiY8nDqqbDhhvCb36RdiSRJagkGpjx07AhnngnXXgsffZR2NZIkqdgMTHn69rdh9Wr4wx/SrkSSJBWbgSlPW2wB3/hGvL/cihVpVyNJkorJwNQM3/8+vPce3Hxz2pVIkqRicluBZho+HF5+Gf7+d2hl/JQkKctKf1uBEMKoEMI8YHbatTTFD34AL70E992XdiWSJKlYHGEqgP794/5Mjz6axrdLkqRGKv0RplJ23nnw2GPwyCNpVyJJkorBEaYCWLsW9tgDttwSHnggjQokSVIjOMKUplat4IILYMoUmF1SK7AkSVJjOMJUIGvWQK9e0KMHTJiQVhWSJOlzOMKUttat4fzz4Z574Nln065GkiQVkiNMBbRqVRxh6tMHbrstzUokSVI9HGHKgrZt4xVzd9wR92aSJEnlwRGmAluxAnbYAYYOhRtuSLsaSZJUiyNMWVFVFXf/vuUWmD8/7WokSVIhGJiK4Fvfgs02g5/+NO1KJElSIRiYiqBDBxg9Gm68Md6YV5IklTbXMBXJ8uWw444waBCMGZN2NZIkiSytYQohtAkh/CyEsCCEsCyEMD+E8OMQQkWNZrVrBxdeCLfeCi+8kHY1kiSpOQo+whRCuAD4LnAS8Hdgb+A64MIkSa5sxPFlMcIEsHIl9OwZ7zN3551pVyNJUsXLzggT0A+YkCTJpCRJFiZJMg6YQgxOFWWDDeDHP4a77oK5c9OuRpIk5asYgekxYGgIoTtACGF3YCBwX31vDiFUhRC6VD+AzkWoKTUnnADdu8fgJEmSSlMxAtMvgLHASyGEVcDTwBVJkoxdz/tHE6fgqh9vFqGm1LRpAxdfDJMmwaxZaVcjSZLyUYw1TMcCvwR+QFzD1Bu4AvhekiSf2fs6hFAFVNV6qjMxNJX8GqZqa9fCbrvBFlvAQw9ByHsGVZIkNUPe/wUuRmB6A7gsSZKraz13IXBCkiQ9G3F82Sz6rm3CBBg+HCZPhoMPTrsaSZIqUqYWfXcA1tZ5bk2RvqtkHHEEDBgA554bR5wkSVLpKEaIuRe4IIRwWAhhuxDCCOB7wPgifFfJCAEuvxyefTbeZ06SJJWOYkzJdQZ+CowANgfeJi4C/0mSJCsbcXxZTslVO+oomDMn3jKlXbu0q5EkqaJkZw1Tc5V7YHr5ZejVC37xC/j+99OuRpKkimJgyslWM+txxhlw++3w6qvQtWva1UiSVDEytehbDbjoIlixAi67LO1KJElSYxiYUvBf/xWn4668El5/Pe1qJElSQ5ySS8nSpbDjjjBkCIwZk3Y1kiRVBKfkSk3nzvDzn8PYsTBzZtrVSJKkz+MIU4rWroU+faBVK3jiifinJEkqGkeYSlGrVnDFFfDUU3DTTWlXI0mS1scRpgw45hh49FH4xz+gU6e0q5EkqWw5wlTKLr8cPvoILr007UokSVJ9DEwZsO22cM458KtfwYIFaVcjSZLqckouIz7+GHr0gH79YNy4tKuRJKkslf6UXAhhVAhhHjA77VrS0KlTnJq7806YMiXtaiRJUm2OMGVIksSNLN96C55/Hqqq0q5IkqSyUvojTIIQ4Oqr4zqmX/4y7WokSVI1A1PG7LwzfO97cMklLgCXJCkrnJLLoI8/hp12gt694d57065GkqSy4ZRcOenUCa68EiZOhHvuSbsaSZLkCFNGJQkceii8+CLMmwcdOqRdkSRJJc8RpnITAvz2t/Duu/DTn6ZdjSRJlc3AlGHdusEFF8Qr5p59Nu1qJEmqXE7JZdzKlbDnntC+PcyaBa1bp12RJEklyym5crXBBnDttTBnDlx1VdrVSJJUmQxMJWDffeE734ELL3RvJkmS0uCUXIlYuhR22QV69oTJk+OicEmS1CROyZW7zp3hD3+IN+a9+ea0q5EkqbI4wlRiRo6MI0wvvgibb552NZIklRRHmCrFFVdAq1Zwxhlxc0tJklR8BqYSs9lmcWpu/HgYMybtaiRJqgyZmZILIYwCRhFDXA+ckvtcxx8P998Pf/87fOELaVcjSVJJyHtKLjOBqZprmBrn3/+GXr3ippYTJ3rVnCRJjeAapkqz8cbw5z/DfffBddelXY0kSeXNEaYSd8opcOed8MILsM02aVcjSVKmOSWXk61mWsDixTUbWk6Z4tScJEmfwym5SrXhhvCXv8BDD8Fvf5t2NZIklScDUxk46KB4r7kf/hCefz7taiRJKj9OyZWJ5cuhT5/49yefhHbt0q1HkqQMckqu0rVrFzeyfOUVOPfctKuRJKm8GJjKyK67wuWXw1VXxU0tJUlSYTglV2aSBA49FObOjeuZvEGvJEn/4ZScohDg+utjcDr1VG/QK0lSIRiYytAWW8TQNGkS/OY3aVcjSVLpMzCVqUMPhXPOiQvAZ81KuxpJkkqba5jK2KpVsP/+8NZb8PTT8f5zkiRVMNcw6bPatoVbb4WPP4aTTnI9kyRJ+TIwlblttoEbb4SJE+FXv0q7GkmSSpOBqQIcdli8bcp558HMmWlXI0lS6cnMGqYQwihgFDHE9cA1TAW1ahUMHgxvvBH3aNp007QrkiSpxeW9hikzgamai76L5403YM89oXfvuBN4mzZpVyRJUoty0bcatvXWcNttMHUqXHBB2tVIklQ6DEwVZsiQeL+5yy+HO+5IuxpJkkqDU3IVKEng+OPh3nvjppa77JJ2RZIktQjXMOVkq5kM++QT6NcPli2DJ5+EjTZKuyJJkorONUxqmo4dYfx4+OADOPFEWLs27YokScouA1MF22EHGDMm3qT3oovSrkaSpOwyMFW4Qw6BSy+Fn/0shidJkvRZrmESSQInnxy3HJg+HfbZJ+2KJEkqChd952SrmRKyYkXccuDVV+Mi8K23TrsiSZIKzkXfap6qqrgIvKoKhg2Djz9OuyJJkrLDwKT/2HzzuDfTq6965ZwkSbUZmLSO3XaLi78nTPD2KZIkVStKYAohfDGEcHMI4cMQwqchhGdCCHsV47tUeMOGwS9/CZddBn/8Y9rVSJKUvoLfrz6E0BWYAfwNOAR4H9gBWFTo71LxfO97sHAhnHkmfPGLcPjhaVckSVJ6Cn6VXAjhMmBAkiSD8jzeq+QyYs0aOPpomDIFpk2DPn3SrkiSpGbJzrYCIYR5wAPAVsD+wFvANUmS/Hk9768Cqmo91Rl4EwNTJixbFrcbmD8fHn8ctt8+7YokScpbprYV2B74H+AV4GDgD8BVIYRvrOf9o4kjStWPN4tQk/LUvj3ccw906RJ3Bf/ww7QrkiSp5RVjhGkl8FSSJP1rPXcV0CdJkn71vN8RphLwz39Cv37QvTs89FAMUpIklZhMjTC9A8yr89yLwDb1vTlJkhVJkiypfgBLi1CTmqlbN5g4EZ5+Go45BlatSrsiSZJaTjEC0wygR53nugOvFeG71IL22QfuugsmT4ZTT3VjS0lS5ShGYPoNsG8I4fwQQrcQwvHAacDVRfgutbCvfAVuugluuQXOPjveuFeSpHJX8H2YkiR5MoQwArgU+DGwADg7SZJbCv1dSscxx8CiRXDGGbDJJnDRRWlXJElScRU8MAEkSTIRmFiMz1Y2nH46/PvfcP750LUrnHVW2hVJklQ8RQlMqgznnRdD0//+bwxNJ56YdkWSJBWHgUl5CwEuvxw++ghOOQU6dYIRI9KuSpKkwivKzXdVOUKIN+j96lfj2qaJTsRKksqQgUnN1ro13HwzDBsWg9PkyWlXJElSYRmYVBBt28LYsXDwwTB8eNwNXJKkcmFgUsFssAHccUe8We8RR8C0aWlXJElSYRiYVFBVVXE38IED4fDD4bHH0q5IkqTmMzCp4Nq1g7vvhr59487gjzySdkWSJDWPgUlF0aED3Hsv7LtvDE0PP5x2RZIk5c/ApKLp2DGGpv33h8MOg/vvT7siSZLyY2BSUbVvH6fnqq+eu+eetCuSJKnpMhOYQgijQgjzgNlp16LCqqqKV88dcUTcp2ncuLQrkiSpaUKSJGnXsI4QQhdgMbBhkiRLmnh4tprROlavhpNOgltvhRtvhJEj065IklRhQr4Hei85tZg2bWJQqqqKN+r96CP49rfTrkqSpIYZmNSiWreGa6+Frl3hO9+BDz6Aiy6K96STJCmrDExqca1awf/9H2y2GYweDR9+CFdeGZ+XJCmLDExKRQhw3nmw8cZwxhkxNF1/fby9iiRJWWNgUqpOOy2GppEjYdGieAVdhw5pVyVJ0rqcBFHqjj4aJk2Kt1A58MA42iRJUpYYmJQJX/4yTJ0K//gH9O8P8+enXZEkSTUMTMqMvn3h8cdh7dp4D7rZbmEqScoIA5MypVu3GJq6dYPBg72ViiQpGwxMypxNN4WHH4ZDD4URI+Dqq9OuSJJU6QxMyqT27eH22+Hss+Nu4OecE6fqJElKg9sKKLNatYJf/Qq23TYGp4UL4YYboGPHtCuTJFUab76rkjBhQtyrqXv3uK5pq63SrkiSVILyvhGXU3IqCUceCTNmxD2a+vSBJ55IuyJJUiUxMKlk7L573Gpg++1h//3hllvSrkiSVCkMTCopW2wRN7g85hg44QQ4/3wXg0uSii8zi75DCKOAURji1ICqqnij3l12gXPPhXnz4MYboUuXtCuTJJUrF32rpN17b1wM/sUvwvjx0LNn2hVJkjLMRd+qTMOGwZNPQgjx1irjx6ddkSSpHBmYVPJ69IhXzR18MBx1FIweDWvWpF2VJKmcGJhUFjp3jjuDX355fBxySNyCQJKkQnANk8rOww/DscfGHcHvugv23DPtiiRJGeEaJqna0KEwZw5sthn07w9/+hNk7P8XSJJKjIFJZWmbbeDRR+HUU+H00+H442FJU8crJUnKMTCpbLVrB9dcA7fdBpMmwV57wdNPp12VJKkUGZhU9r7+dZg7Ny4M79cPfv97p+gkSU1jYFJF6NYNZs6Eb34TzjwzLgpfvDjtqiRJpcLApIrRrh387ndwxx0weXKconvqqbSrkiSVAgOTKs7RR8e1TF27xim6Sy91o0tJ0uczMKkibb89zJgBP/gBXHABDBkCr7+edlWSpKwyMKlibbAB/PznMG0aLFwIu+0GY8emXZUkKYsMTKp4++0Hzz4Lhx4a92s68UQXhEuS1mVgkoCNNoIxY+Dmm+Gee2D33eGxx9KuSpKUFZkJTCGEUSGEecDstGtR5Ro5Mo42bb11HHk65xxYtiztqiRJafPmu1I91qyB3/wGLrwQttsObrgB9tkn7aokSc3kzXelQmrdOo4uPf00dOkSb+I7ejSsWJF2ZZKkNBiYpM+x005xh/BLLoFf/9rNLiWpUhmYpAa0aQPnnQdz5kBVFey7L/zoR442SVIlMTBJjbTLLjBrFlx0EVx2Gey5Zxx9kiSVPwOT1ARt28bRpblzoVMnGDgQRo2CJU29PEGSVFIMTFIedt01ji5dcUW8gm7nneP+TZKk8mRgkvLUujWcdRb8/e9xo8sjj4SvfQ3efTftyiRJhWZgkppp221h4sS4U/j06fHKumuvhYxtcSZJagYDk1QAIcBxx8GLL8Lw4fCtb8HgwfDCC2lXJkkqBAOTVECbbALXXQcPPRSn5nr3jhtgLl2admWSpOYwMElFMHQoPPcc/OQncM010LMn3H6703SSVKoMTFKRVFXB+efDvHnQty8ccwwcdBC8/HLalUmSmsrAJBXZdtvB+PFxYfj8+XFLgvPPh08+SbsySVJjFT0whRBGhxCSEMIVxf4uKcsOOywuAj///Hhfup4945V1TtNJUvYVNTCFEPoApwHPFfN7pFLRvj1cfHHNNN3IkXG3cG/oK0nZVrTAFELoBNwCfAv46HPeVxVC6FL9ADoXqyYpK7bfHu68E6ZOjVfQ9ekDp5wC77yTdmWSpPoUc4TpamBSkiQPNfC+0cDiWo83i1iTlCkHHBDvS/f738O990L37vCLX8CKFWlXJkmqrSiBKYRwLLAnMQw15FJgw1qPrYpRk5RVbdrAGWfAK6/Af/83XHhhvDfdnXe6vkmSsqLggSmEsDVwJXBCkiTLG3p/kiQrkiRZUv0A3OJPFalr13gz3+eegx494OijYcAAmDEj7cokScUYYdoL2ByYE0JYHUJYDewPnJX759ZF+E6pbOy0E9x3Hzz4ICxfHheFH3WU+zdJUpqKEZgeBnYFetd6PEVcAN47SZI1RfhOqex8+cvx6rmbb4Y5c6BXLzjzTHjvvbQrk6TKE5IWWCQRQpgGPJMkydmNeG8X4uLvDXNTdE3hig+VpeXL4Xe/g0sugdWr4Qc/gO9/Hzp2TLsySSopId8D3elbKgHt2sWb+L76Kpx2WgxO3brB1VfDypVpVydJ5a9FRpiawhEmqWELF8KPfgS33ALbbgsXXQQnnBCvuJMkrZcjTFIl2W47uOkmeP552HPPuOnlrrvCHXfA2rVpVydJ5cfAJJWwXr3ifk1PPhlHmr7+ddhrr3iVXcYGjyWppBmYpDKw994weTI88gh07hxv9DtwIEyblnZlklQeDExSGRk0CKZPj+FpxYp465UDDoC//c0RJ0lqDgOTVGZCgIMPjtN048fD4sUwZAjsvz88/LDBSZLyYWCSylQIMHx43PTy3nth2bK4GebAgTBlisFJkprCwCSVuRDg8MNh9uy4GHzNmjgC1a8f3H+/wUmSGsPAJFWIEOCQQ+Dxx+GBB6B1azj0UOjbN45AGZwkaf0MTFKFCQEOOggeewweegjat4cjjoDdd4cxY+KtVyRJ6zIwSRUqBBg6NG5FMH06bLUVjBwJO+4Yb7ny6adpVyhJ2ZGZwBRCGBVCmAfMTrsWqdLst19c3/TMM3Ft01lnxd3EL7kEPvoo7eokKX3eS07SZ8yfD7/6Ffz1r/H+dKefDt/9Lnzxi2lXJknNkve95AxMktbrvffgqqtqpui+8Q045xzo2TPtyiQpLwamnGw1I5WJJUvgT3+CX/8a3nknXl333e/GNVAh73/9SFKLy/vfWJlZwyQpu7p0iSNLCxbADTfA22/DgQfGK+uuuy7ehkWSypmBSVKjVVXFabm5c2HqVNh2Wzj11PjnT38K//pX2hVKUnE4JSepWV5+Ga68Eq6/Pm5+eeKJcPbZsPPOaVcmSZ/hGqacbDUjVZAPP4zrnH7727jO6Stfge98J/7ZyrFsSdngGiZJ6dpkExg9GhYuhBtvhPffh8MOg+7d42Jx93OSVMocYZJUFEkCTzwBv/sd3H573M/phBPg29+G3XZLuzpJFcopuZxsNSMJiPs5/elP8Ic/xCvsBg2K03XDh0PbtmlXJ6mCOCUnKbu22AJ+9KM4XXf77XHvpq9/Pd5+5Sc/iSFKkrLMESZJqXjuubiD+M03x32cjjgi3oLlwANdJC6paJySy8lWM5IatHhxDE1//CM8/zx86UvwrW/BKafAllumXZ2kMmNgyslWM5IaLUlg1qwYnG67DVavjmucTj8dhgxx1ElSQRiYcrLVjKS8fPQR3HRTDE/z5sEOO8Bpp8HJJ8Pmm6ddnaQSZmDKyVYzkpolSWDGjBic7rgD1q6FYcPi7VgOPjhuVSBJTWBgyslWM5IK5sMP41qnv/41Lhj/whfgpJPiWqcdd0y7OkklovQDUwhhFDCKuNVBDwxMkuqRJPD00zE43XILLFoU93U69VQ4+mjo1CntCiVlWOkHpmqOMElqrOXL4e674S9/gYcfho4d4ZhjYnjq1y/u9yRJtRiYcrLVjKQWs3Ah3HADXHcdvPYa9OgRp+tGjoSttkq7OkkZYWDKyVYzklrc2rXwt7/FKbu77oqbYg4ZAieeCEcdBZ07p12hpBQZmHKy1YykVC1ZAnfeCTfeCNOmQYcOMGJEDE9f/jK0bp12hZJamIEpJ1vNSMqM11+Pi8RvvBFeeinuIj5yZAxPu++ednWSWoiBKSdbzUjKnCSBOXPixphjxsAHH8Cuu8I3vgHHHx+3K5BUtgxMOdlqRlKmrVoFDzwQw9OECfGfBw+GY4+Fr34VNt447QolFZiBKSdbzUgqGYsWxfVOY8fGReOtW8fdxI89Fo480v2dpDJhYMrJVjOSStK778ZbsYwdC48/Du3bw+GHw3HHwSGHQLt2aVcoKU8GppxsNSOp5C1cCLfdBrfeCs88A126xCvtjjsOhg71fnZSiTEw5WSrGUll5cUXY3AaOxZeeQU23RS+9rU4bTdggNsUSCXAwJSTrWYklaXq+9mNHRtHn954I25TcNRR8X52gwY58iRllIEpJ1vNSCp7a9fCE0/AuHHx8frrsNlmNeFp8GDDk5QhBqacbDUjqaIkCTz5ZFwwPm5cXP+0ySZxzdPRR8dbtLRtm3aVUkUzMOVkqxlJFat6g8xx42KAmj8funaF4cPjuqehQ2GDDdKuUqo4BqacbDUjScTw9MwzNeHplVdgww3jVgUjRsT9ntznSWoRBqacbDUjSXUkCTz/fAxPd98d/15VBQceGMPTsGFxDZSkoij9wBRCGAWMAloBPTAwSaoAr74ab8syfjzMmAEhxC0KRoyI03df+lLaFUplpfQDUzVHmCRVqvffh3vvjeHpoYdgxQrYffcYnIYPj38Pef/rXhIGpv/IVjOSlKelS+ONgcePh0mTYPFi2G67mvA0YIDbFUh5MDDlZKsZSSqAlSth+vQYniZMgLffjtsVHHpoXPN00EFxEbmkBhmYcrLVjCQV2Nq18NRTccH4xIlx0XibNrD//vGqu8MPh27d0q5SyiwDU062mpGkInvttRicJk6EqVPjaFTPnjXhyak7aR0GppxsNSNJLejjj+Ni8YkT47qnd9+FjTaCQw6J4emQQ+LmmVIFMzDlZKsZSUrJ2rUwd2686m7ixPj31q3jiNOwYTFA9ejhVXeqOAamnGw1I0kZ8dZbcdRp4sQ4CrVsGeywQxx1OuSQeJPgDh3SrlIqOgNTTraakaQMWrYsrneaNAnuvz/eJLiqKi4crw5Q3bs7+qSyZGDKyVYzkpRxSQIvvxyD0+TJcfuCFSviDuNf+UoMT0OGQMeOaVcqFYSBKSdbzUhSifnkE5g2LQao+++H+fNhgw1g0KCa0aeddnL0SSXLwJSTrWYkqYQlCfzznzXhado0WL4cttmmZvRp6FDo3DntSqVGMzDlZKsZSSojy5bVjD5NngyvvBL3eNp337jb+IEHwt57u++TMi07gSmEMBo4CugJLANmAucmSfJyI483MElSCXj1VZgyJT6mToUlS+K+T0OHxvB00EFxLZSUIZkKTJOBW4EngTbAJcCuwM5JknzSiOMNTJJUYlavhtmz4cEHY4B64glYsyZuXVA9+nTAATFQSSnKTmD6zBeEsBnwPrB/kiSP1PN6FVBV66nOwJsYmCSpZC1eDH/7WwxPDz4Y10K1bg19+9YEqH32cfpOLS7Tgakb8Aqwa5IkL9Tz+sXARfUcamCSpDKxYEHN6NPDD8OiRdClSxx1qg5Q3bp59Z2KLpuBKYQQgAlA1yRJBq3nPY4wSVIFWbMGnnqqJkA9/nic0ttuO/jyl+O+TwccAFtumXalKkOZDUxXA4cBA5MkebORx7iGSZIqyNKl8eq7Bx+Mo0/z5sXne/WK4Wno0LgLueufVADZC0whhN8Cw4H9kiRZ0ITjDEySVMHefTdedTd1agxQCxdCq1aw1141AWrAAO99p7xkJzDlpuF+C4wABidJ8koTjzcwSZL+Y8GCGJyqQ9R778Xdx/v1qwlQfftC27ZpV6oSkKnAdA1wPHAkUHvvpcVJkixrxPEGJklSvZIkTtlVjz5NmxavyOvYEfbbLwaoIUOgd+84KiXVkanAtL4PPCVJkusbcbyBSZLUKGvWwNy5NQHqscfijuQbbwyDB8fRpwMOgJ49vQJPQJYCU3MZmCRJ+VqxIm6aWT2FN2tWvAJviy3iwvHBg+Of3kC4YhmYcrLVjCQpVR9/DDNnxqm7adPgySdjgNp883UD1M47G6AqhIEpJ1vNSJIy5ZNP1g1Qs2fHALXZZp8NUK6BKksGppxsNSNJyrRPPokbZ06bBtOnx+m8Vatg003XDVC9ehmgyoSBKSdbzUiSSsqnn64boGbNigFqk03WDVC77GKAKlEGppxsNSNJKmmffhpDU+0AtXJlvApv//3jY7/9YLfd4s2FlXkGppxsNSNJKivLlsXQNH16DFGzZsUr87p0ibuPDxoUA9Tee0NVVYMfp5ZnYMrJVjOSpLK2fHm88u7RR+GRR+KC8qVLoV072GefGKAGDYq7knfunHa1wsD0H9lqRpJUUVavhmefrQlQjz4KH3wQp+v22COOPg0aBAMHxoXlanEGppxsNSNJqmhJAi+9tG6Aev31+NrOO9dM4Q0aBFtvnW6tFcLAlJOtZiRJquO112Jwqg5RL70Un99uu5opvP32g+7d3UyzCEo/MIUQRgGjgFZADwxMkqQK8P778R541QHqmWdg7dq4G/nAgTUjULvtBm3apF1tySv9wFTNESZJUiVbsiQuHq8OULNnx60MOnWKi8cHDIhBap994nNqEgNTTraakSSpmaqvxJsxI45EzZgBixbFheS9e9cEqAED4AtfSLvazDMw5WSrGUmSCmztWnjxxZrw9NhjsGBBfO1LX6oJTwMHwk47uSN5HQamnGw1I0lSC3j77RieqgPUM8/AmjXQtSv0718ToPr0iXtEVTADU062mpEkKQUffxxvJFw9CvX44/G5DTaAvfaqGYUaMKDi9oMyMOVkqxlJkjJg9Wp4/vmaAPXoo3FUCqBHjxigqkNUt25lvZ2BgSknW81IkpRBSRL3g6q9kPyFF+Lzm29eM4U3cGDcobxt27QrLhgDU062mpEkqUQsWhSn7qoD1BNPxCv02reHvn1jiOrfP25tsPHGaVebNwNTTraakSSpRK1cCXPn1iwmnzEjbrIJ8eq7/v1rFpSX0K7kBqacbDUjSVKZSBKYPz9uqjlz5rrTeJtsUrOpZv/+sPfe0KFD2hXXy8CUk61mJEkqY4sXx6m7GTNiiJo1K16N16YN7LnnuqNQGdlU08CUk61mJEmqIGvWxFGn6gA1c2bNpprbblsTnvr3h113TeXeeAamnGw1I0lShXv77biYvHoab+5cWLUKOnaM98OrDlD77gsbbVT0cgxMOdlqRpIkrWPZMpgzZ91RqA8+iIvGe/VadxRqhx0KvpjcwJSTrWYkSdLnShJ45ZV1F5PPmxdf22yzdQPUXns1+9YuBqacbDUjSZKa7KOP4gLy6lGoJ56ATz+tubVL9WLy/v1hyy2b9NEGppxsNSNJkppt9Wp47rl1p/Fefz2+tv32MTj98peNCk+lH5hCCKOAUUAroAcGJkmStB5vvlkTnmbNgqlTG7X3U+kHpmqOMEmSpCLJOzC1KmQVkiRJ5cjAJEmS1AADkyRJUgMMTJIkSQ1o+bu4FFdh9wOVJEkim1fJBaAzsDTJWnGSJKkiZS4wSZIkZY1rmCRJkhpgYJIkSWqAgUmSJKkBBiZJkqQGGJgkSZIaYGCSJElqgIFJkiSpAQYmSZKkBhiYJEmSGmBgkiRJakDZ3Hy31j3oJEmS1ieve9WWTWAihqXFaRchSZIybUNgSVMPKpub7zZhhKkz8CawFbC0CV8xG+jbxLKaekyxv6Oles/izyqf3luij3yOKZfey+X3PZ9jyqX3LP6ssvr7ns8x5dJ71n7fK3uEKdd8g4kx5iog/sAanTBDCGub8v58jin2d7RU71n8WeXTe0v0kc8x5dJ7ufy+53NMufSexZ9VVn/f8zmmXHrP6u97U7nou/GuboFjWuI78pHFPrLYd0sdUy69l8vvez7HlEvvWf1ZNZX/Wy/uMVn9fW+SspmSa6wQQhfiWqcNi5lEs8je7b2Seq/UvsHesfeK6r2l+q7EEaYVwP/L/Vlp7N3eK0ml9g32bu+VpUX6rrgRJkmSpKaqxBEmSZKkJjEwSZIkNcDAJEmS1AADkyRJUgMMTJIkSQ0oi8AUQjgzhLAghLA8hDAnhDDoc977XyGEMSGEl0MIa0MIV6znfV8NIcwLIazI/TmiaA00Q6F7DyGcHEJI6nm0K2ojTdTEvo8KITwYQvhXCGFJCOHxEMLB9byvHM95g72XyjmHJvc+MIQwI4TwYQhhWQjhpRDCd+t5Xzme9wZ7L5Xz3pS+6xw3IISwOoTwTD2vld05r3Ncvb2XyjmHJv++D15PXz3rvK955z1JkpJ+AMcAK4FvAjsBVwAfA9us5/3bAVcC3wCeBq6o5z39gNXAaKBn7s9VwD5p99sCvZ9M3ABsy9qPtHttZt9XAD8E+gA7Aj/PHb9HBZzzxvSe+XOeZ+97AMcBvXK/+ycAnwCnVcB5b0zvmT/vTe271nEbAq8CDwDP1HmtLM95I3vP/DnPp3dgMJAA3ev01rqQ5z31H0wBfrBPAL+v89yLwKWNOHYa9YeG24D76zw3GRibdr8t0PvJwKK0eytW37Xe/3fgx5V0zj+n98yf8wL2fhdwU4We97q9Z/6859s3cCvwU+BiPhsayvqcN9B75s95Pr1TE5g2+pzPbPZ5L+kpuRDCBsBewJQ6L00B+jfjo/vV85kPNPMzC6qIvQN0CiG8FkJ4M4QwMYSwRzM/r2AK0XcIoRXx7tb/rvV0RZzz9fQOGT7nULDe98i9d3qtpyvlvNfXO2T4vOfbdwjhFGAH4s7P9Snbc96I3iHD5xya/fv+dAjhnRDCwyGEA+q81uzzXtKBCdgUaA28V+f594jDcfnasgifWWjF6v0l4v8LOYI4pL8cmBFC2LEZn1lIhej7+0BH4PZaz1XKOa+v96yfc2hG77n/MKwAngKuTpLk2lovl/V5b6D3rJ/3Jvedq/0yYGSSJKvX87llec4b2XvWzznk9/v+DnAa8FXgKOBl4OEQwn613tPs896msW/MuLr3dwn1PJeFzyyGgtaZJMksYNZ/PiyEGcBc4DvAWfl+bhHk1XcI4TjiUPWRSZK8X4jPTEFBey+hcw759T4I6ATsC1wWQvhnkiRjm/mZaSho7yV03hvVdwihNTAGuChJkn8U4jMzoKC9l9A5hyacoyRJXiaGpGqPhxC2Bs4BHsnnM+tT6oHpA2ANn02Im/PZJNkU7xbhMwutWL2vI0mStSGEJ4kLhrMg775DCMcAfwG+liTJQ3VeLutz3kDv68jgOYdm9J4kyYLcX58PIWxBDI3Vgamsz3sDvdd9b9bOe1P77gzsDewRQvhd7rlWQAghrAYOSpJkKuV5zhvb+zoyeM6hcP9tm0W82KFas897SU/JJUmyEpgDHFjnpQOBmc346Mfr+cyDmvmZBVXE3tcRQghAb+KQZ+ry7Ts3unI9cHySJJPqeUvZnvNG9F73/Zk651DQ3/cAVNX657I97/Wo2/u6L2bsvOfR9xJgV2IP1Y8/EEceehMXEkN5nvPG9r6OrJ1zKOjv+x6s21fzz3vaq+ELsJq++vLDU4mXH/6GePnhtrnXLwVurHNM79zjKeCW3N93rvV6f+Llh+cSLz88l2xfdlrI3i8CDga2z73211zvfdPuN9++iXP1q4AzWfeS0w3L/Zw3svfMn/M8ex8FDCP+v+cdgVOIl1T/rALOe2N6z/x5b2rf9Rx/MZ+9Uqwsz3kje8/8Oc+nd+BsYHjud71X7vUEOKqQ5z31H0yBfrhnAguBFcRkul+t164HptV5f1LPY2Gd9xxNXCC3kng541HF7iMLved+MV/Lfd77xKsI+qXdZ3P6Jm6hUF/f15f7OW9M76VyzvPo/TvAC8T9hxYT12r8D9CqAs57g72XynlvSt/1HHsxdUJDuZ7zxvReKue8qb0T95r7J7CMeAXwo8ChhT7vIfchkiRJWo+SXsMkSZLUEgxMkiRJDTAwSZIkNcDAJEmS1AADkyRJUgMMTJIkSQ0wMEmSJDXAwCRJktQAA5MkSVIDDEySJEkNMDBJkiQ14P8DZrc81vbEYSoAAAAASUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 1 graphics primitive" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# This point is minimal as H(u,x) = 0 and x in (0,1) implies u increases as x decreases\n", "solve(H,u)[0].rhs().plot(x,0.1,1/2)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The minimal point [u == 1, x == (1/2)] is critical in the direction [1/4, 1]\n" ] } ], "source": [ "# Find the direction for which CP is a critical point, scaled to have final coordinate 1\n", "m = [k.subs(CP) for k in [u*diff(H,u)/x/diff(H,x),1]]\n", "print(\"The minimal point\", CP, \"is critical in the direction\", m)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "As n goes to infinity 1/2^n * c(s,n) approaches\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|ν(s)|\\phantom{\\verb!x!}\\verb|=|\\phantom{\\verb!x!}\\verb|\t| \\frac{\\sqrt{5} \\sqrt{2} e^{\\left(-\\frac{{\\left(n - 4 \\, s\\right)}^{2}}{10 \\, n}\\right)}}{5 \\, \\sqrt{\\pi n}}\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|ν(s)|\\phantom{\\verb!x!}\\verb|=|\\phantom{\\verb!x!}\\verb|\t| \\frac{\\sqrt{5} \\sqrt{2} e^{\\left(-\\frac{{\\left(n - 4 \\, s\\right)}^{2}}{10 \\, n}\\right)}}{5 \\, \\sqrt{\\pi n}}$$" ], "text/plain": [ "'ν(s) = \\t ' 1/5*sqrt(5)*sqrt(2)*e^(-1/10*(n - 4*s)^2/n)/sqrt(pi*n)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Get density for limit theorem\n", "detM = getHes(H,m,[u,x],CP).determinant()\n", "EXP = (s - m[0]*n)^2/detM/2/n\n", "C = -(G/x/diff(H,x)).subs(CP)\n", "nu = C * exp(-EXP) * (2*pi*n)^(-1/2)/sqrt(detM)\n", "print(\"As n goes to infinity {}^n * c(s,n) approaches\".format(x.subs(CP)))\n", "show(\"ν(s) = \\t \", nu)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The following plot shows c(s,n)/2^n vs computed density when n = 200\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAGFCAYAAAAPa6wiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAA9hAAAPYQGoP6dpAABBU0lEQVR4nO3deXxU1f3/8dchG2vCDklIIGxCBAU1KloVrUvVXzdtq9a2tmpb29SltFqxi35bK37r2tb41dbWpdatLl2lIFatC2pYVDAgQQIJWdhJWANJzu+PM0NCyJB15s6d+34+HvMYcufMvZ/JnZm8Offcc421FhERERGJrJfXBYiIiIjEOwUmERERkXYoMImIiIi0Q4FJREREpB0KTCIiIiLtUGASERERaYcCk4iIiEg7FJhERERE2hGowGScdGOM8boWERER8Y9krwvoQe1OWV5bW0tGRga1tbWxqEdERETiR7c6SwLVwyQiIiLSFQpMIiIiIu1QYBIRERFphwKTiIiISDu6FJiMMd81xpQZY/YaYxYbY05pp/1poXZ7jTFrjDFXtXr8m8aY140x20K3BcaY41u1ucUYY1vdarpSvwTHSy/Beec08sUvwpo1XlcjIiJ+1enAZIy5CLgX+CUwHXgdmGuMyY3QPg94MdRuOnAb8BtjzIUtms0EngROB2YA5cB8Y0x2q9V9CGS2uE3tbP0SHGVl8MT5f+Yv89O5/dlx3PuJZ70uSUREfMpY2+7Z+Ac/wZh3gCXW2u+0WLYC+Ku1dnYb7f8X+Iy1dnKLZQ8AR1trZ0TYRhKwDfietfax0LJbgM9Za6dFKK3dF1JXV3dgWoH09PT2motPNTbCb38L/PGPXLfsCp7iIvqym3OZS+MHJfSeOsHrEkVEJPZiN62AMSYVOBaY3+qh+cBJEZ42o43284DjjDEpEZ7TF0gBtrZaPsEYUxU6HPiUMWZsx6uXoJg9G370/Xq+tOwn/JkvcwlPchFPsyk5k5UX3MTcuV5XKCIiftPZQ3JDgSRgQ6vlG4CREZ4zMkL75ND62nI7UAksaLHsHeBrwDnAN0PrfcsYMyRSsfX19dTV1R10k8T30kvwVf7ESGr4OT+joMAw4/Q+zG74BdNWP8vPz3tboUlERDqlq2fJtT78ZdpY1l77tpZjjLkBuAS4wFq798AKrJ1rrX3OWrvMWrsAOD/00GWRNjpnzhwyMjIO3HJycg5ToiSKY45u5AZ+xQt8nlUcwY03QkMDPM5XWMUEvs0DzJvndZUiIuInnQ1Mm4FGDu1NGs6hvUhhNRHaNwBbWi40xvwQuAk421r7weEKsdbuApYBEQekzJ49m9ra2gO3ioqKw61SEsR9F7/BREp5vWAW998PF1wARx8NTSTxF77IZ/g7R+fv97pMERHxkU5dS85au88Ysxg4C3ihxUNnAX+L8LSFwKdbLTsbWGStPfBXyxhzPfAT4Bxr7aL2ajHGpAGTcWfftSktLY20tLT2ViUJps/L/4QRI7j37RkH/kvwq1+5+7VvXsDgpbfxjXH/BT7pWY0iIuIvXTkkdzdwpTHmcmPMZGPMPUAu8ACAMWaOMeaxFu0fAEYbY+4Otb8cuAK4M9wgdBjuVuByYK0xZmTo1r9FmztD8znlGWNOAJ4F0oFHu/AaJJH9619w/vnQq/nt3aePO3Pu94uPgdxceP55DwsUERG/6XRgstY+DVwH/Ax4DzgVOM9auy7UJBMXoMLty4DzcHMtvQf8FLjGWvtci9V+F0jFhaDqFrcftmgzCjdX00fA88A+4MQW2xWh5q01sGKFC0xtMcYdo3vhBWhqim1xIiLiW52ehymOaR6mAGtqgksugeHP/Ja7+AF//+MWvvCNAW03fvVVOP10WLwYjjkmpnWKiIhnYjcPk0i8mjcPnnkGzmY+r3MKV34/QlgC9h9zAk1Jyax6dGEMKxQRET9TYJKE0NAAhiZO5k1e4zQaGiK3O/eCPixqnM67v1nIFVfEtk4REfEnBSZJCJ/6FFx+4goGs42F5mTuuKPtdosXw8svw0JmMIOF/PGPsGVL221FRETCFJgkIaSkwO8uexOblMQflp/Ad77TdruMDHe/kBmMYw2jUjfSp0/s6hQREX9SYJKE0WvJIsyRR5Kb3z9im0mT3JxMi1PcdZ+fuHohffvGqkIREfErBSZJHEuWdOist+uvh4/25GIzMzkl5e0YFCYiIn6nwCSJYd8+WLasw9ME9EoymBkzYKHOlBMRkfYFIjAVFRWRn59PQUGB16VItJSUuNB07LEdf86MGVBcDI2N0atLREQSQiACU2FhISUlJRQXF3tdikTLkiVuFu+jj+74c6ZPh927YfXq6NUlIiIJIRCBSQJg8WI3ortfv44/Z+pUd798eXRqEhGRhKHAJIlh6VLXY9QZw4fDsGFu7JOIiMhhKDCJrzU0wD13W/Ys/pCy/lM7v4KpU9XDJCIi7VJgEl+79lq44wfV9NlXxw8emkynh6lNnaoeJhERaZcCk/javHkwiZUALG+azCuvdO75q1Kn0FS6msVv7IlCdSIikigUmMTXpk2DyaxgHymsYWynTpJ7+mn42h1T6WWb+M7MFbz8ctTKFBERn1NgEl976CG4cPIKqvpN4IHfJ3POOR1/7uOPw4fkAzC5cRlPPhmlIkVExPeSvS5ApDsGDoQzMldA/mSuvLJzz83NhZ0MYA15TGE5e0dHpUQREUkACkzifytWwBVXdPppt90GNTWw7sVJfHLYKo68IQq1iYhIQlBgEn+rrYXqapg8udNPzciA554DrhkPL78MaT1fnoiIJAaNYRJ/W7HC3XchMB0wfjx8/DE0NfVMTSIiknAUmMTfVqxw15A74oiur2P8eKivh8rKnqtLREQSSiACU1FREfn5+RQUFHhdivS0FStg9Gjo27fr6xg/3t3rIrwiIhJBIAJTYWEhJSUlFHd6GmiJeytWdO9wHMCYMZCUpMAkIiIRBSIwSQJbtap7h+MAUlNdL1Vpac/UJCIiCUeBSfyrsRHKypoPqXXH+PHqYRIRkYgUmMS/Kith/34YO7b761JgEhGRw1BgEv9as8bd92Rgsrb76xIRkYSjwCT+tWaNm1JgdA9c02TCBNizx02CKSIi0ooCk/jS2rXw9pNr2DUoG5vWu/srDI+D0sBvERFpgwKT+M769VBQAGsWrGHR1rFcfXUPrDQvz/VWaRyTiIi0QYFJfGf+fNi8GcayhjWM5Yknur/O8g1pbOmfy4IHVlNV1f31iYhIYlFgEt/Jy3P34/iYjxl34Oeu2r0bTj0Vlu4Yz7ZFq5k5E/bt63aZIiKSQBSYxHdOPx1+/Ys6hrGZlIljefLJ7q1v9WpYtw5KmcB4VlNaCuXlPVOriIgkBgUm8aVrPl0GwM2PjmXixO6ta/RoGDIEVjOeCZQycoQlK6sHihQRkYQRiMCki+8moB6cgykjA156CYaeMJ7+7OKVJ2u6dS1fERFJPMYmzkR97b6Quro6MjIyqK2tJT09PRY1SbTcdRfcfDPs2OHObusJy5bBUUfBW2/BjBk9s04REYkX3fpjEYgeJklAa9a43qWeCksAubnuXgOYRESkFQUm8aePP4Zx43p2nRkZkJ6uwCQiIodQYBJ/WrOGbs8n0JbRo90pcyIiIi0oMIn/WOt6gcaM6fl15+aqh0lERA6hwCT+s2kT1Nc3jznqSaNHKzCJiMghFJjEf8KBJhqBKTdXh+REROQQCkziP+HAlJPT8+vOzYXt26GurufXLSIivqXAJP5TXg69e8PQoT2/7tGj3X1FRc+vW0REfEuBSfynosL1BPXkHExh4cN8OiwnIiItKDCJ/5SXR2f8EkBmJiQna+C3iIgcRIFJ/Ke8PDrjlwCSkmDUKPUwiYjIQQIRmHTx3QQTzR4m0FxMIiJyiEAEpsLCQkpKSiguLva6FOmu+nqoqVFgEhGRmApEYJLEYC387f8qAahKjmJg0uVRRESkFQUm8Y0f/ADu+b7r+fnc1TmsXRud7SzelEvT+kr+8mRDdDYgIiK+o8AkvvHEE5CLC0zL63KYN6/nt/HYY/Dj3+XSyzbxgy9X8fDDPb8NERHxHwUm8Y2xYyGHCjYxlD30JS+v57cxdy6sw01eOZp1zJ3b89sQERH/UWAS33j8cTgxq5wNabncfjucfXbPb2PqVKjATVmQSzlTp/b8NkRExH+SvS5ApKPGjoWxR5VDWg5TfhSdbdxwA+zY0Z/aOwfz1ZPKOXN2dLYjIiL+oh4m8Zcoz8GUnAxz5kDG1NF8avI6kvVfChERoYuByRjzXWNMmTFmrzFmsTHmlHbanxZqt9cYs8YYc1Wrx79pjHndGLMtdFtgjDm+u9uVBGNt9CetDNNcTCIi0kKnA5Mx5iLgXuCXwHTgdWCuMabNv2LGmDzgxVC76cBtwG+MMRe2aDYTeBI4HZgBlAPzjTHZXd2uJKDaWti5MzaBadQoWL8++tsRERFf6EoP0yzgD9bah6y1K6y11wEVwHcitL8KKLfWXhdq/xDwR+CH4QbW2kuttfdba9+z1q4Evhmq7ZPd2K4kmnCPT7SuI9dSdjZUVkZ/OyIi4gudCkzGmFTgWGB+q4fmAydFeNqMNtrPA44zxqREeE5fIAXY2o3tSqKpqnL32dmHb9cTsrNh61bYsyf62xIRkbjX2SGtQ4EkYEOr5RuAkRGeMzJC++TQ+qrbeM7tQCWwoBvbpb6+nvr6+gM/19XVRWoqfhAOTCMj7vKeEw5lVVUwblz0tyciInGtq2fJ2VY/mzaWtde+reUYY24ALgEusNbu7c5258yZQ0ZGxoFbTiwO5Uj0VFXB0KGQmhr9bY0a5e41jklEROh8YNoMNHJor85wDu39CauJ0L4B2NJyoTHmh8BNwNnW2g+6uV1mz55NbW3tgVtFRUWkpuIH1dWQlRWbbYV7mDSOSURE6GRgstbuAxYDZ7V66CzgrQhPW9hG+7OBRdba/eEFxpjrgZ8Cn7LWLuqB7ZKWlkZ6evpBN/GxqqrYBab+/SE9XYFJRESArs30fTfwJ2PMIlwY+haQCzwAYIyZA2Rba78Wav8A8D1jzN3A73GDwK/AHXYj9JwbgF8AXwbWGmPCPUk7rbU7O7JdCYDqasjPj932NLWAiIiEdDowWWufNsYMAX4GZALLgfOstetCTTJxQSbcvswYcx5wD1AIVAHXWGufa7Ha7wKpwLOtNvc/wC0d3K4kuqoqOPPM2G1PUwuIiEhIly78YK29H7g/wmNfb2PZa8Axh1nfmO5uVxJcU5PrYcrMjN02s7Nh5crYbU9EROKWriUn/rBlCzQ0xG4ME7jApENyIiKCApP4RXgOplgGplGjXK9WY2PstikiInFJgUn8oTo0v2msD8k1NsLGjbHbpoiIxCUFJvGHWM7yHRaei0mH5UREAk+BSfyhqgqGDYvNLN9h4dm+daaciEjgKTCJP8T6DDmAoUOxKSlsfl+BSUQk6BSYJO6tWgWlr1VRbWI44Bv46c29WLc/i9/fsp6rr47ppkVEJM4EIjAVFRWRn59PQUGB16VIJ61eDQUFsPXDKua+n8mNN8Zmuxs2wK23QiXZZFPJfffBRx/FZtsiIhJ/AhGYCgsLKSkpobi42OtSpJP+9S+oq4NMqqkiiyefjO321zOKbNwhOWNiu20REYkfgQhM4l9jx4Kh6UBgysuLzXZHjICbb27uYbr2Wpg4MTbbFhGR+KPAJHHt05+GO2/cQgoNDD8qk0cfjd22b7kFLv9JNhP7rufee2zsNiwiInFHgUni3qyL3RxMt/wui9GjY7vtgUdm02v3LtixI7YbFhGRuKLAJPEvPGllrKcVaLnN8EzjIiISSApMEv/CYSWWs3yHhQNTOLSJiEggKTBJ/PNilu8w9TCJiAgKTOIHXszyHTZgAPTvr8AkIhJwCkwS/6qqICu2s3wfJCtLh+RERAJOgUnin9eBKTNTPUwiIgGnwCTxz8tDcqDAJCIiCkwS55qaXFjRITkREfGQApPEt82boaFBPUwiIuKpQASmoqIi8vPzKSgo8LoU6axwUPF6DNOOHbBzp3c1iIiIpwIRmAoLCykpKaG4uNjrUqSzwofCvD4kB+plEhEJsEAEJvGxcEgZMcK7GjR5pYhI4CkwSXzzcpbvMAUmEZHAU2CS+FZV5e2Ab4D0dOjbV2fKiYgEmAKTxDevpxQAMEZnyomIBJwCk8Q3r2f5DsvMVA+TiEiAKTBJfIuHQ3LgQpt6mEREAkuBSeJXUxPU1MRFD9OujEx2llazbp3XlYiIiBcUmCR+hWf59jgwffgh3Pl4Jg0VVeTnw6uvelqOiIh4QIFJ4lf4EJjHh+Tuvx9W78liILU07d7DXXd5Wo6IiHhAgUnikrXw36fcIOsNSd72MPXvD9W40JZJNenpnpYjIiIeUGCSuFRYCI/e7gLTCZ8ZwcaN3tVy440wbKoLTDNyq7jtNu9qERERbwQiMOniu/7z6KOuN2cjw1hXncqCBd7VMmgQPPma6+X6853VjB7tXS0iIuKNQAQmXXzXf3JyIIsqqsg68LOnBg6EtDRNLSAiElCBCEziP888A5MHVrOtdyb/+79wyikeFxSe7VuTV4qIBFKy1wWItOWoo4CJVTBlCqff4HU1IZq8UkQksNTDJPErXi6LEqbryYmIBJYCk8Sn8Czf8XBZlLCsLB2SExEJKAUmiU9xMsv3QdTDJCISWApMEp/CPTnx1MOUmQlbt0J9vdeViIhIjCkwSXwK9+TEUw9TuBb1MomIBI4Ck8SncA/TyJHe1tFSuLdLgUlEJHAUmCQ+VVfDsGGQkuJ1Jc0UmEREAkuBSeJTvE0pADBkiAtwOlNORCRwFJgkPlVVxdeAb2ie7Vs9TCIigaPAJPGpujr+ephAgUlEJKACEZiKiorIz8+noKDA61Kko+LxkBxo8koRkYAKRGAqLCykpKSE4uJir0uRjojHWb7D1MMkIhJIgQhM4jPxOMt3mAKTiEggKTBJ/Akf8orHwJSVBZs2wb59XlciIiIxpMAk8SfcgxOvh+QANmzwtg4REYkpBSaJP/E4y3dYODBp4LeISKAoMEn8qaqKv1m+Q9bsdYcJf/Hdat54w+NiREQkZhSYJP7E6RxMTU1w5sVD2U8yVUuqOfdcjf8WEQkKBSaJP3E6B9PWrVC2rhcbGEEm1ezcCatWeV2ViIjEQpcCkzHmu8aYMmPMXmPMYmPMKe20Py3Ubq8xZo0x5qpWjx9pjHnOGLPWGGONMde1sY5bQo+1vNV0pX6Jc9XVcTnge8gQmD4dqskkk2pGjICpU72uSkREYqHTgckYcxFwL/BLYDrwOjDXGJMboX0e8GKo3XTgNuA3xpgLWzTrC6wBbgQOF4I+BDJb3PTnKhHFaQ+TMfDSS9BvfCYnjXFjmAYP9roqERGJheQuPGcW8Adr7UOhn68zxpwDfAeY3Ub7q4Bya+11oZ9XGGOOA34IPAdgrS0GigGMMbcfZtsN1lr1KiWwyoomMqtrWLc3kzyvi2nDkCEw5IxMKC6G8V5XIyIisdKpHiZjTCpwLDC/1UPzgZMiPG1GG+3nAccZYzp7GtQEY0xV6HDgU8aYsYdrXF9fT11d3UE3iV9r1sCZ0zbTq7GB79+ZxYMPel1RBJrtW0QkcDp7SG4okAS0nrVvAxBp0pyREdonh9bXUe8AXwPOAb4ZWu9bxpghkZ4wZ84cMjIyDtxycnI6sTmJtaeegrStbn6jKrK47z6PC4okMxM2boTGRq8rERGRGOnqWXK21c+mjWXttW9reeQVWDvXWvuctXaZtXYBcH7oocsiPWf27NnU1tYeuFVUVHR0c+KBYcMgE9dzU00mw4d7XFAkmZlujoGNG72uREREYqSzY5g2A40c2ps0nEN7kcJqIrRvALZ0cvsHWGt3GWOWARMitUlLSyMtLa2rm5AYu/xy6PVwFSyE4VNH8sADXlcUQfgMvjg9m09ERHpep3qYrLX7gMXAWa0eOgt4K8LTFrbR/mxgkbV2f2e235IxJg2YDGgwSYJISoIrzq2C4cNZ/EEKEyJGYY+1DEwiIhIIXTlL7m7gT8aYRbgw9C0gF3gAwBgzB8i21n4t1P4B4HvGmLuB3+MGgV8BXBJeYWgweX7ox1Qg2xgzDdhprV0danMn8A+gHNdD9RMgHXi0C69B4pUfem1GjHBzDCgwiYgERqcDk7X26dBA65/h5kJaDpxnrV0XapKJC1Dh9mXGmPOAe4BCoAq4xlr7XIvVZgFLW/z8w9DtNWBmaNko4EncQPFNwNvAiS22K4kgTudgOkhKihtwpcAkIhIYXelhwlp7P3B/hMe+3say14BjDrO+tTQPBI/U5uJOFSn+VFXlj+mzNbWAiEig6FpyEl/i9MK7h1BgEhEJFAUmiR9NTVBTo8AkIiJxR4FJ4sfmzdDQEP+DvkGBSUQkYBSYJH5UuVm+fdPDVFMDtsNzr4qIiI8pMEn88Ftg2rcPtm71uhIREYkBBSaJH9XVbn6jESO8rqR94cOG4ZAnIiIJLRCBqaioiPz8fAoKCrwuRQ6nqsrNb5SS4nUl7dNs3yIigRKIwFRYWEhJSQnFxcVelyKH44dZvsMUmEREAiUQgUl8wg+zfIf17g0DByowiYgEhAKTxA8/BSbQ1AIiIgGiwCTxw0+H5ID9QzOpWVpNWZnXlYiISLQpMEl8aGryz2VRgDVr4G/FWZS+Xs3kyfDvf3tdkYiIRJMCk8SHTZugsdE3PUy/+x2s2ZtJJtXU18Ptt3tdkYiIRJMCk8SH8Fggn/Qw9e8P1bjABJb+/b2uSEREokmBSeKDn2b5Bq69FjKOyKQfu5mSu4O77vK6IhERiSYFJokPfprlGxgwAG550B0+XDa/miOO8LggERGJKgUmiQ9+muU7TJNXiogEhgKTxAe/zcEECkwiIgGiwCRxYWdpNdv7ZtLY6HUlnTBgAPTrp8AkIhIACkziuTvvhJKXq3j2rSzOOw8aGryuqBMyM5sHrIuISMIKRGAqKioiPz+fgoICr0uRVpqa4Mc/hkyqqSaT+fNhwQKvq+qErCz1MImIBEAgAlNhYSElJSUUFxd7XYq0YgykJjeRSTVVuDFMvXt7XFRnZGdDZaXXVYiISJQFIjBJ/DIGHrljE8k0UkUW3/gGzJzpdVWdkJWlQ3IiIgGQ7HUBIhee5A5pPfLvTAad43ExnZWV5XqYrHXpT0REEpJ6mMR7oR6aQUf6bFoBcIfkdu+GujqvKxERkShSYBLvVVX5apbvg4TnjtJhORGRhKbAJN6rrvbfLN9h2dnuXgO/RUQSmgKTeK+ysjl4+E14tm/1MImIJDQFJvFeZaX/LosS1qcPDBqkwCQikuAUmMR7fu5hAs3FJCISAApM4j2/BybNxSQikvAUmMRb+/bBxo3+D0zqYRIRSWgKTOKt8HXY/ByYsrPVwyQikuACEZh08d04Fu6Z8XNgCl+At6nJ60pERCRKAhGYdPHdOBbumfFxYNozOBsaGrilcBOrVnldjYiIREMgApPEscpK6N3bnZrvU7PudFMi/O2BKk4+GTZs8LggERHpcQpM4q3wHEw+vXBtfT38fbELTNlUsnkzLFrkcVEiItLjFJjEWz6fUiAtDQZOHEEjvciiitRUOOIIr6sSEZGepsAk3vJ5YAL4x9xkanuP4OQxVbzwAowf73VFIiLS05K9LkACrrISjjvO6yq6ZexY4MgsLpteCed5XY2IiESDepjEO9YmRA8ToLmYREQSnAKTeGf7dtizJzECky6PIiKS0BSYxDsJMAfTAboAr4hIQlNgEu+EA0ZWlrd19ISsLNi0yV0bT0REEo4Ck3gn0QITNF8bT0REEooCk3iioQH+WlTJFjOUM85N83/OCB9W1DgmEZGEpMAknnjgAaheXEmFzeaVV+Caa7yuqJvCPUwKTCIiCSkQgamoqIj8/HwKCgq8LkVC1q93lxKpxPXM+H689ODBbtpv378QERFpSyACU2FhISUlJRQXF3tdioRcfDHkmObA9I1veFxQdxmjqQVERBKYZvoWT0ybBg1DKtl//Gd45XqYOdPrinqAApOISMJSYBJv7N9P8paNHP/5bJjpdTE9RHMxiYgkrEAckpM4VFPjLo2SCFMKhKmHSUQkYSkwiTfCPTGJMMt3WFaWephERBKUApN4IxEDU3Y27NjhbiIiklAUmMQblZXuNPwhQ7yupOdotm8RkYSlwCTeqKx0AcMYryvpOeHApMNyIiIJp0uByRjzXWNMmTFmrzFmsTHmlHbanxZqt9cYs8YYc1Wrx480xjxnjFlrjLHGmOt6YrsSxyorE+twHBx4PZveq2TXLo9rERGRHtXpwGSMuQi4F/glMB14HZhrjMmN0D4PeDHUbjpwG/AbY8yFLZr1BdYANwI1PbFdiXMJGJjqGvtRmzSIu2dVMGoUvPmm1xWJiEhP6UoP0yzgD9bah6y1K6y11wEVwHcitL8KKLfWXhdq/xDwR+CH4QbW2mJr7fXW2qeA+h7arsSzqqqEC0y/+x2sbcwhhwq2b4frr/e6IhER6SmdCkzGmFTgWGB+q4fmAydFeNqMNtrPA44zxqREcbvU19dTV1d30E3igLXNY5gSSGMjVOACE0BDg8cFiYhIj+lsD9NQIAnY0Gr5BmBkhOeMjNA+ObS+aG2XOXPmkJGRceCWk5PTwc1JVNXVwa5dCdfDdOWVsGvQKEaxnr594dZbva5IRER6SlfPkrOtfjZtLGuvfVvLe3S7s2fPpra29sCtoqKik5uTqEjEOZhwMyRccG0OUwdWsHYtnH221xWJiEhP6ey15DYDjRzaqzOcQ3t/wmoitG8AtkRxu6SlpZGWltbBTUjMJGhgAkjJy4HtWxjWbzfuXAYREUkEnephstbuAxYDZ7V66CzgrQhPW9hG+7OBRdba/VHcrsSrcGBKsDFMAIQP+65f720dIiLSo7pySO5u4EpjzOXGmMnGmHuAXOABAGPMHGPMYy3aPwCMNsbcHWp/OXAFcGe4gTEm1RgzzRgzDUgFskM/j+/odsVHKivd8avevb2upOeFA5MO/4qIJJTOHpLDWvu0MWYI8DMgE1gOnGetXRdqkokLMuH2ZcaY84B7gEKgCrjGWvtci9VmAUtb/PzD0O01YGYHtyt+UV4OuQk6fdaoUe5egUlEJKF0OjABWGvvB+6P8NjX21j2GnDMYda3luaB4F3arvhIIgem3r1h2DAFJhGRBKNryUnsVVQkbmAC18ukMUwiIglFgUliy1pYty6xA1NOjnqYREQSjAKTxFZtLezcycbeCkwiIuIfCkwSU3/6ZTkAn706h8sv97iYaFFgEhFJOApMEjN79sBf7nKBqZxcHn4Y3nnH46KiISfH9aTt2OF1JSIi0kMUmCSmcihnP8nUhCZtt529OI4faC4mEZGEE4jAVFRURH5+PgUFBV6XEmh9+sBXZ1awnlE0kcRXvgInnuh1VVGgwCQiknCMTZz/4rf7Qurq6sjIyKC2tpb09PRY1CStXXop9R+vp/xPrzFhgtfFRMm+fW4+pt//Hq64wutqRETEaXe+x8MJRA+TxJHyctIm5CZuWAJITYURI9TDJCKSQBSYJLbKy5sPWSUynSknIpJQFJgkdhob3YV3E3nSyrCcHDa/V8ETT8D27V4XIyIi3aXAJLFTXe1CUwAC0+trc9i0pIJLL4UTTlBoEhHxOwUmiZ1yNwdTEALTvz7IIYcKwLJqFcyf73VFIiLSHQpMEjsBCky1A0fTn10MZivgxoCLiIh/KTBJ7JSXQ0YGBGBKh6tuHwPAlH5r+clP4LTTvK1HRES6J9nrAiRAgnKGHHD05/LgSnjtkTL4wrFelyMiIt2kHiaJnfLyQByOA2DwYBgwAMrKvK5ERER6gAKTxE5FRXACkzEwZgysXet1JSIi0gMUmCR2gtTDBJCXpx4mEZEEocAksbFzJ2zdqsAkIiK+FIjAVFRURH5+PgUFBV6XElzhKQUCMugbaD4klzgXuBYRCaxABKbCwkJKSkooLi72upTgCve05OV5W0cs5eXB3r2wYYPXlYiISDcFIjBJHCgrg5QUyMryupLYGTPG3Wvgt4iI7ykwSWyUlbkAkZTkdSWxE+5N0zgmERHfU2CS2CgrC9bhOHAzmg8erMAkIpIAFJgkJnZ8UMbbG/N4+mmvK4kxzcUkIpIQFJgk6l78l6Xp4zW88F4eF18Mv/2t1xXFkKYWEBFJCApMEnUvP7uNDOoowx2S+8c/PC4olhSYREQSggKTRF3BUBcYwoFp8mQvq4mxMWPcHFSNjV5XIiIi3aDAJFH3pQIXmDKmjeVb34I5czwuKJby8mD/fm6/pooHH4SmJq8LEhGRrkj2ugBJfL3WlcGAASxYMhiM19XE1qr9eUwE/nX/Wt4gh5IS+PWvva5KREQ6Sz1MEn3hKQVMwNISMLdkNAB5uF62QI3fEhFJIApMEn1BnIMpZOK0vtQw4kBgCtT4LRGRBBKIwKSL73pszZrABqZzzwU7Jo+CoWv50pfg4Ye9rkhERLrC2MS5knq7L6Suro6MjAxqa2tJT0+PRU3S1AR9+sCdd8LVV3tdjTcuuQSqq+HVV72uREQkyLo1LiQQPUzioepq2LcvsD1MgHvtmu1bRMTXFJgkusKTNgY5MI0ZAxUVsH+/15WIiEgXKTBJdIUD05gxnpbhqbw8d2iyosLrSkREpIsUmCS6SkshMxP69fO6Eu+Ee9d0WE5ExLcUmCS6Skth4kSvq/BWbi706gUff+x1JSIi0kUKTBJdq1YpMKWmwujRLjyKiIgvKTBJ9FjrAtOECV5X4r0JExSYRER8TIFJomfDBti5Uz1M4ALT6tVeVyEiIl2kwCTRs2qVu1cPU3NgamryuhIREekCBSaJntJSd8HdceO8rsR748fD3r1QWel1JSIi0gUKTBI1DStWsWPIaH7/WBo7d3pdjcfCvWwaxyQi4ksKTBIV1sI7j63irc0T+da3YOZMd4WUwMrLwyYl8eD1pZx8MjzyiNcFiYhIZwQiMBUVFZGfn09BQYHXpQRGZSVkbCplFW7A9+LFsGyZx0V5KSWFypQx1C1ZzVtvweWXw5tvel2UiIh0VCACU2FhISUlJRQXF3tdSmAMymhiPKspxR2KSkmBkSM9LspjKxsnMAF3SM5aKCnxuCAREemwQAQmib1+WyvoTT27syeSlwePPw7Z2V5X5a2mcc2BqV8/d5hSRET8IdnrAiRBhaYUeOi1CaCT5AD45FUTsT94kB99v5GLL03SbAsiIj6iwCTRUVrqjsONHu11JXEj6chJ0LiP269aq6kWRER8RofkJDpWrYKxYyFZmfyASZPc/YoV3tYhIiKdpsAk0aGL7h4qO9sNXlq50utKRESkkxSYJDpKSxWYWjPG9TIpMImI+I4Ck/S8/fuhrEzXkGvL5Mk6JCci4kNdCkzGmO8aY8qMMXuNMYuNMae00/60ULu9xpg1xpir2mhzoTGmxBhTH7r/fKvHbzHG2Fa3mq7UL1FWWgqNjXDEEV5XEn8mTXKByVqvKxERkU7odGAyxlwE3Av8EpgOvA7MNcbkRmifB7wYajcduA34jTHmwhZtZgBPA38Cjg7dP2OMOaHV6j4EMlvcpna2fomB5cvd/ZQp3tYRjyZNgm3bYPNmrysREZFO6EoP0yzgD9bah6y1K6y11wEVwHcitL8KKLfWXhdq/xDwR+CHLdpcB7xkrZ1jrV1prZ0DvBxa3lKDtbamxW1TF+qXaFu+3E3rPXSo15XEn/CZchrHJCLiK50KTMaYVOBYYH6rh+YDJ0V42ow22s8DjjPGpLTTpvU6JxhjqkKHA58yxow9XL319fXU1dUddJMYWL5cvUuRjB8PSUm6LoqIiM90todpKJAEbGi1fAMQ6UphIyO0Tw6t73BtWq7zHeBrwDnAN0OPvWWMGRKp2Dlz5pCRkXHglpOTE6mp9KQPP1RgiiQtDSZMwC7/kLIy2L7d64JERKQjunqWXOsRq6aNZe21b738sOu01s611j5nrV1mrV0AnB966LJIG509eza1tbUHbhUVFYcpUXrEnj2wejUceaTXlcStpvwpvP/n5YwdC5mZ8Le/eV2RiIi0p7OBaTPQyKG9ScM5tIcorCZC+wZgSzttIq0Ta+0uYBkQ8dz1tLQ00tPTD7pJlK1cCU1N7J2gHqZIljOFrG1uYPzevXDNNR4XJCIi7epUYLLW7gMWA2e1eugs4K0IT1vYRvuzgUXW2v3ttIm0TowxacBkoLr9yiVW6t5yQWD4zHyOPhqqqjwuKA5ty57CcDYxjI2Am4FBRETiW1cOyd0NXGmMudwYM9kYcw+QCzwAYIyZY4x5rEX7B4DRxpi7Q+0vB64A7mzR5tfA2caYHxljJhljfgSciZu+gNB67wzN55QXmm7gWSAdeLQLr0Gi5N2Hl7OW0ewgnQ8+gJtv9rqi+HPiFe5w5ZF8SEoK3HlnO08QERHPdfrKqNbap0MDrX+GmwtpOXCetXZdqEkmLkCF25cZY84D7gEKgSrgGmvtcy3avGWMuRi4FfgF8DFwkbX2nRabHgU8iRsovgl4GzixxXYlDgzbsJzlNB+O04mJh0o7cjw2NZU/z1qOueZ0MjO9rkhERNpjbOLMONzuC6mrqyMjI4Pa2lqNZ4qS+szR3L/1Embtu51+/WD+fDgp0oQTQTZtGpxwAjz4oNeViIgEhWm/SWSd7mESiaiujrSaci6790gmTXQzC2gmhwimTIEPPvC6ChER6SAFJuk5ockYB586hXOne1xLvJs2DV54wY34TkryuhoREWlHV+dhEjnU8uXQq1fz5T8ksunTYfdud6FiERGJewpM0nM++MBd+qNPH68riX/TQ11w773naRkiItIxCkzSc5YsgWOO8boKfxg8GHJzYelSrysREZEOUGCSntHU5HpLFJg6bto0BSYREZ9QYJKeUVoKu3YpMHXG9OnYpUv59b2WWbNg4UKvCxIRkUgCEZiKiorIz8+noKDA61IS15Il7n66To/rsOnTMZs386vvV3HPPTBzpoY0iYjEq0AEpsLCQkpKSiguLva6lMS1dCmMHu3G5kjHhMLldNxhuX374D//8bIgERGJJBCBSWJAA747LyeHupTBBwITwNSpHtYjIiIRKTBJ9zU1YYuLeS/lOJ580vWUSAcYQ+9PHMeF2e9w6qnuKilnneV1USIi0hbN9C3d1rishKS6Or7/zAxefQY++Ul3DbleiuPtSj11BtM+uI/XXrVgunWZIxERiSL9SZNu2/i3hTTSi2LcoPqXX4aPP/a4KL848UTYsgVWr/a6EhEROQwFJum2gSsWsswcxS76A5CWprHfHXbCCe7+7be9rUNERA5LgUm6rc97C+l/5gyys2HUKHj8cRgyxOuqfGLQIHftPU3CJCIS1zSGSbpn61ZYuZLxN93E+q96XYxPzZihwCQiEufUwyTd88477n7GDG/r8LMTT3QXLt61y+tKREQkAgUm6Z6FC2HoUBg3zutK/GvGDGhq4rGrixk61B2he/ddr4sSEZGWFJikexYudD0kOiW+6/Lzaeg7gBUPL2TLFvjoI7joIq+LEhGRlhSYpOsaG90hOR2O656kJLaMPZ4TaT5TrqbGw3pEROQQgQhMuvhulJSUwI4dCkw9IP2cGZzS600MTQB8+9seFyQiIgcx1lqva+gp7b6Quro6MjIyqK2tJT09PRY1Jbb77oNZs2DbNujXz+tq/O3VV+H00/nH/yyh17HTOf98rwsSEUk43Ro7EogeJomS//yHnVNO5A9P9ePDD70uxudmzIC+ffl075cUlkRE4pACk3RNYyP75r/CXe9/kiuvhGOPhf/+1+uifCwtDU49FRYs8LoSERFpgwKTdM3SpaTu2s6CpjMAqK+HRx7xtiTfO+sseP112LvX60pERKQVBSbpmv/8h/rkvrzDCQcWjRzpYT2J4MwzXVh64w3q692ZcokzxFBExN8UmKRrXn4Zc9qpnH52KgMHwmc/Czfd5HVRPjd1KowYQdWfFpCdDZmZcMopmgBcRCQeKDBJ5+3eDf/9L6nnncW8ee4kub/+Ffr397ownzMGzjyTHc+9xJYtbtGbb8KDD3pbloiIKDBJV7z6qjt0dN55XleSeM48kwm7ljKEzQcW1dd7WI+IiAAKTNIV//wnjBkDRxzhdSWJ56yz6IXlU8kvAzB+PFxxhcc1iYiIApN0UlOTO/72+c/r+nHRkJ0N06bxu/P/ysKF8N57MHy410WJiIgCk3TOwoVQXQ0XXuh1JYnrwgvp+/I/OXHaXk2gLiISJxSYpHOee469gzM5+YczmDkTFi3yuqAE9IUvwM6dMH++15WIiEhIIAKTLr7bQ6yl4S/P82jt53nr7V689pob961ByT1s0iQ48kh49lkeeQQmT4YTToAlS7wuTEQkuHTxXem4RYugoIAzeJlXOOPA4spKyMrysK5EdMstNN59L/12bqTepgJueNP69R7XJSLiX7r4rsTIX/6CHTyEqnGnHlh0wgma4TsqLryQpB21nGGbry1XVQX79nlYk4hIgCkwScc0NMBjj2EuvojX3kzm5z+H22+Hl16CXnoX9bwpU2g8Ip9v9338wKILLoDUVA9rEhEJMB2Sk47529/gc5+DpUth2jSvqwmGO+7A/vSn3P+TanpnDuKyyyA52euiRER8S4fkJAYeegiOO05hKZa++lVMQwOFg5/kiitcWNq/3+uiRESCSYFJ2ldZCS++SMU5V/L441BW5nVBATFyJJx/Ptx/Px+vtkya5A7JnX22LsgrIhJrCkzSvocfpiGlN0fNuYSvfhWOOsrNQC0xMGsWfPghj3/l33z0kVv00kvwm994W5aISNAoMMnh7d0LRUXMG3op25vcuK+dO+Hhhz2uKyhOPRUKCvj0yjsOWrx9uzfliIgElQKTHN4f/wgbN7Jg+g0HLdb1zWLEGLj+eo6pfYXjkxcDMHQoXHmlx3WJiASMzpKTyPbvhwkT4KSTqLrzCS68EN5/342hefJJ6NPH6wIDorERJkygbvLxvHbVU0yZAm+/7X7/n/mMpnUQEemgbp0lp8AkkT3yCHzjG7BsGUyZ4nU1wXbffXDttdQvW8UpXx9HcbFbfMkl8MQT3pYmIuITCkwhCkw9ae9emDyZTaOm8fX0F+jXz01UOXas14UF1O7dMHEimyaexPBXnjnooa1bYdAgj+oSEfGPbgUmTYMnbfv1r7Hr13NWzb95f69btHQplJZ6W1Zg9e0Lv/wlw77+dU7mTd7kZAD694d+/TyuTUQkAALRw1RUVERRURGNjY2sWrVKPUztWbUKpk1j3ae+zZgX7jnoobo6GDDAo7qCrqkJjj+eLVX1TNm7iF590jj5ZKithTPOgBtucGPERUSkTTokF6JDcj2hsRFOOQU2bWLTS++RX9CPzZvdQ8cfD++84215gffBB27G9e9/n5uS/pc5c5of+r//g6uu8q40EZE4p0Ny0oPuvhv79ts8d+3r/OdX/bjtNjfmu39/uP56r4sTjjoKfvELmD2b+oJPA5848NCSJd6VJSKS6NTDJM1efx3OPJP/TruG0951EyUa42aW/uQnPa5NmjU2wqmnUruiignb3mETblKs4cMhPd3NAn7uuR7XKCISf3TxXekBK1a4SX1OPpkbG249sNhaWLDAw7rkUElJ8MQTZKTtZfmY/8cVF+8iKQk2boTVq+GLX4QdO7wuUkQksSgwibu47rnnsmdwNt/Lep6d+9MOeviYYzyqSyIbPRpefJHhm1dw1/ovkdRYf+ChXbvc5KI/+YnrjBIRke7TIbmgW7kSzj2Xhn2NTNn+Jh/tzgEgKwvy8+Gzn4Xvfc/jGiWy+fOxn/kMxWmf4My659nBwe/rH/8Yrr4aRozwqD4RkfihQ3LSBdbCQw9hjz2WrXv6cOURbxwISwDV1TBvnsJS3Dv7bMy8eRRQTPmIAi6dfPDI71/+EkaOdMFJRES6ToEpiDZswH7hi/DNb/LfnK+Qs6GYR1/JPajJccfpGmW+cdppmHffZWBWPx5ddSI/5lZ6s+egJrfd5i7ae/75sGULNDR4VKuIiE/pkFyQ1NSw5cY7GPjU/7GrsTffaPg9z3PhQU2OPBJOPBFuvdX1TIiP1NfDz35G0113sz1tBLN3/5Q/cym76H9Qs+TQZCLXXOP299SpUFDgQb0iIrGliStDAh2YrHUDfI2BtWth2DA3L8/2jzYwaOl/SPnbsxRs+Ae7bR/u5Tru5Tq2c+gFyP75T9cLIT62ejX89Kc0PfU0u+nLX/gif+ZS3uAT1NP7kObGNIfjz37W3U+aBBMnwvbtbhbx9evduLaBA90ZeMOHw759kJoas1clItJdCkwhEV9IQwMUFsK//lVHZWUG+fm17NmTztix7tpoWVmwZ4+77Me4ce7KIMOHuwCybRuMH+/aDRniDlNt2uT+mJSWQkYG9O7txvxMnAgff+yu7ZWeDhUVbtnatZCW5p6/dq1bX0WFW9fIkVBWBnl5UFPjrn4xapT7mzdmjNvWvn2Qk+O2l5MDmze7ekeNcrMBZGXB5poGBuyo4oSBHzF820qm9FrBjKY3OIplACxhOo/xNR7lskOCUu/erlfpoos0U3RCKS9n812PkvbkwwzYVMYeerOQGSxlOh9wFCuYzHpGsYERNJEUcTWpqe49mJzsbnv3uvfh+vXuor/9+rn3aX6+ey8PHOiW1dS44FVW5n7OyHDPmTgRysvdeocMcf8OfyaSktwA9bIyd6Hnlp+Jjz92Jwdu3gz797t/l5ZCbq4Ldnv2uM9RaSlkZ7uzBevq3LpXrXKftf37D/5MDx3qAuOmTTBhgvvcDRzoaqupccvWrHETt/bv704oHT8e1q1zn+nBg92/J0xwryM52f1nZe1a911SVeV+h1lZrv68PDf9Q0OD+x2uXu3q37bN/V7z8lyto0bBzp0unLasv77eXQpnwgS3bNgw9/vZsqX5O2nQIFfHhg3N30kDBrjLEVZVueeWlbmfMzJc3eH6U1LcOsvKXP2Vle73k5npfg/hfdLY2Fz/6NFu+/X1zfXn5Ljf/a5dHPiezcx0bbZvd9v76CP3Pduy/lWr3O80PE3GEUc0f8/26ePqnzjR1dK376Hfs6mpbp+Gv2fXr2/+no1U/5gxzfWPGeNqyM11v+eW9WdluWtg79jR/HdixAi3LzvynkpLi/x3YsKEQ99TbX0mxo1z62hqcu/x8Htq0yb33s7N7fpnIvyeall/+D3Vlc9EuP7kZLefw/V35DMxerS7qHjLz0ROjvvdt/xMZGa6NrW17uSWWbPa/VZUYDLGmNra2qbWy+vr66mvr+ehh+BX/7ObL/InHuNWLuPnpOFOnXe/PRv6d/O9aWNZVx+L3L7tbbd+Xi+a6MMe+rCHvuyhD7tD/97NILYxnI0MZuuBAWn1pLCa8SzhGF5lJq9xGpto+zSp1FT43e/g859v77csvtXUxOq/l7Dl+dfo9fZbDN7wIeMoO/BwA73YwAg2M5Sd9GcHA9jBAHbSn3rSaCSJ/aTQQPKB+waSaWwRsuyBd3PH7lt+b9kI6+gJdWTwDBf12PpEJH79+98wY0bkxzMyMjKAHbaLwSdRAlM6UOt1HSIiIhLXMqy1dV15YqIEpsP2MC1eDF/4AjQ0VAPHAyVAdhQrKgCKo7h+tw1jimm9+3r3dl2UqamuG3vHDjjnHNclP2oUXHedG+/Ut28HtlBQQHFxtF9H9LdTV1dHTk4OFRUVUR27FovfVzS20dQEf/+7e6/U1rqZ3YuLCxg0qJht29whnooK996qr3fvn/BhupZ69XLraskYDnmPtng1RPdzUgfkABVANMcsxuLzHqvtJMI2YrXfITF+X4mxjZEj4cUX6zjmmMjf9d3tYUqIi++29+JHjYJFi+Cvf4VbboF77hlASko606dDcbE7Jrt9uztmeuyxbtno0e4PQk2NO8W+uNgdA05Kcsf6jzsOFi92x/v793fHYgsK3EDrOXOSeOihdD76yC17/30XUEaNguXL3XOXL3fjBcaNg/feg+nT3TFZa91YkMWL4eij3XHfffvcNVcXL4YpU9zx3uuuS2LevHSWLnXjRJqa3HHok05yx48zM11t3ZGUlBSTwfGx2k56enpUtxOL1xGtbXzta83/vukmyM9PoqTEbaepyY07GDzYfR42bXLv23ffbR7/V17u3q/vvuved717u7EIxx7rPnsDB7rnr1zpPhNLl8IvfpHEn/6UzrJlB38mxo51n5np0904F2th8mT3/j/qqObPxNSpbt1TprhxFDt3cuAzPXGiG7/yrW/B/PnplJS4MYt79rixPeG6Ro1yr3n9+ubP9IgRbrzMmjXNn+lBg1z9q1a5du+95z732dnwla8k8eyz7nWkproxFx984GbIX7nShcaJE91rPvpot979+90ZiuHP9Pr1bnzM0Ue7+o84wo2r2b7drae4GO66K4lf/SqdTZuav5Nyc914nMpKV+uiRe4PR2qq+x4oKHDbGDzY7YPS0ub6BwxwbVescNtYtgxuvTWJxx9PP1D/ihUuCIfrnzbN7deGhub6p051+3/PnubvqUmT3PdUba3bJ4sWuXEndXUwa1YSzz2XTnGxGzO0f797fx13nGuXleXeB2vXNv/+hwxx9bb8ns3IcGNjVq50z33/fbffhgxxJy/MnZvO2rXpJCe7bb/3nqu/tNS9p8Pfs0cd5bZVX+9ey+LF7v1WU+P+E9HW34ljjnG13nFHEnffnU5NTfN7qq2/E8OHuzFLLesfONDV2t7fiZ//PIk//zmd99932w1/JiZNcuvpzGdi61b3dyL8nho71r3vrr02ieefT++Rz0RWFnz4oVu2bJkb1zRmDHz5y0k880z6Ie+pjn4mJk1y4xdbfibGj3f7aNMmuPBCt32I/F3f1Z6llitIlFu7KioqLGArKio60rzL7rvvvqiuP5G2EYvt1NbWWsDW1tZGdTvaJ/G1jUTa77HaTiJsI1b73drE+H0l0jY6sO+7lTMS4pBcSLsvZP369QcOzYwKx2hJeIk8nYREpv0eTNrvwdWBfa9Lo3RUWlraQfcSDGlpadx8883a7wGj/R5M2u/BFe19H6geJv3PQ0REJLA0D1NHtZh+IMN2d/CXiIiIBEbQApMBBtCN0wpFREQkeAIVmERERES6IlCDvkVERES6QoFJEoIx5hZjjG11q2nxuAm1qTLG7DHGvGqMOdLLmqXzjDGnGmP+EdqP1hjzuVaPt7ufjTFpxpjfGmM2G2N2GWP+bozRPCNxrgP7/pE2vgPebtVG+95HjDGzjTHFxpgdxpiNxpi/GmOOaNUmZp95BSZJJB8CmS1uU1s8dgMwC/gebo7+GuAlY8yAWBcp3dIPeB+3H9vSkf18L/B54GLgE0B/4J/GmCQknrW37wH+zcHfAee1evxetO/95DSgCDgROAt3dZL5xph+LdrE7jPf3ZkvddMtHm7ALcB7ER4zQDXwoxbL0oDtwLe9rl23Lu9zC3yuM/sZyAD2ARe1aJMFNALneP2adOvavg8tewT462Geo33v8xswLLTvTw39HNPPvHqYJJFMCHXLlhljnjLGjA0tzwNGAvPDDa219cBrwEke1CnR0ZH9fCyQ0qpNFbAcvRcSwczQoZtVxpjfG2OGt3hM+97/MkL3W0P3Mf3MKzBJongH+BpwDvBN3IfoLWPMkNC/ATa0es6GFo+J/3VkP48E9llrtx2mjfjTXOBS4AzgB7jDM/8xxoSnfda+97HQtEB3A29Ya5eHFsf0M5/cmcYi8cpaO7fFj8uMMQuBj4HLgPDAz9ZzaJg2lon/dWU/673gc9bap1v8uNwYswhYB5wPPH+Yp2rf+8N9wFG4MUitxeQzrx4mSUjW2l3AMmACbhAgHPq/ieEc+j8T8a+O7OcaINUYM+gwbSQBWGurcYFpQmiR9r1PGWN+C3wGON1au77FQzH9zCswSUIKdcNPxg0ILMN9aM5q8Xgq7gyMtzwpUKKhI/t5MbC/VZtMYAp6LySU0OH4HNx3AGjf+05oyoD7gAuAM6y1Za2axPQzr0NykhCMMXcC/wDKcf9z+AmQDjxqrbXGmHuBm4wxpUApcBOwG3jCm4qlK4wx/YHxLRblGWOmAVutteXt7Wdrba0x5g/AXcaYLbjBo3fieiMXxOyFSKcdbt+HbrcAz+EC0hjgNmAz8AJo3/tUEfBl4LPADmNMuCep1lq7pyPf7T26370+TVA33XriBjwFVOFOH63EfXHmt3jc4L5Qq4G9uLMopnhdt26d3s8zceMOWt8e6eh+BnoDvwW2hL5Y/wHkeP3adOv6vgf6APOAjaHvgHWh5Tmt1qF976NbhP1tga+3aBOzz7yuJSciIiLSDo1hEhEREWmHApOIiIhIOxSYRERERNqhwCQiIiLSDgUmERERkXYoMImIiIi0Q4FJREREpB0KTCIiIiLtUGASERERaYcCk4iIiEg7FJhERERE2qHAJCIiItKO/w+FfU9rUxW0vQAAAABJRU5ErkJggg==\n", "text/plain": [ "Graphics object consisting of 201 graphics primitives" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Plot series terms versus computed density\n", "K. = QQ['U']\n", "P. = K[['X']]\n", "\n", "# Set the value of n to test\n", "N = 200\n", "mser = (1 - X)/(1 - 2*X - (U-1)*X*(1-X) + O(X^(N+1)))\n", "uvals = mser[N]\n", "\n", "plt = point([])\n", "for k in range(N):\n", " plt += point([k,uvals[k]/2^N])\n", "\n", "print(\"The following plot shows c(s,n)/2^n vs computed density when n =\", N)\n", "plt + plot(nu.subs(n==N),s,0,N/2,color='red')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Question 5\n", "See Section 4.1.4 of [An Invitation to Analytic Combinatorics](https://melczer.ca/files/Melczer-SubmittedManuscript.pdf) for the solutions to (a) - (d) and [the corresponding Sage notebook](https://melczer.ca/files/TextbookCode/Chapter4/Example4-3-5-6-GeneratingWalkGroupsAndGFs.html)." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The critical points in the main diagonal direction are\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left[x = 1, y = 1, t = \\left(\\frac{1}{4}\\right)\\right], \\left[x = \\left(-1\\right), y = \\left(-1\\right), t = \\left(-\\frac{1}{4}\\right)\\right]\\right]\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[\\left[x = 1, y = 1, t = \\left(\\frac{1}{4}\\right)\\right], \\left[x = \\left(-1\\right), y = \\left(-1\\right), t = \\left(-\\frac{1}{4}\\right)\\right]\\right]$$" ], "text/plain": [ "[[x == 1, y == 1, t == (1/4)], [x == -1, y == -1, t == (-1/4)]]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Introduce variables and function\n", "var('x,y,t')\n", "vari = [x,y,t]\n", "R = [1,1,1]\n", "\n", "ST = [[-1,0],[1,0],[0,-1],[0,1]]\n", "S = add([x^i*y^j for (i,j) in ST])\n", "G = (1+x)*(1+y)\n", "H = 1-t*x*y*S\n", "\n", "# There are two critical points\n", "d = len(vari)\n", "criteqs = [H] + [R[j]*vari[0]*diff(H,vari[0]) - R[0]*vari[j]*diff(H,vari[j]) for j in range(1,d)]\n", "print(\"The critical points in the main diagonal direction are\")\n", "show(solve(criteqs,vari))" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dominant asymptotic behaviour of the main diagonal is\n" ] }, { "data": { "text/html": [ "\\[\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{4}{\\pi \\left(\\frac{1}{4}\\right)^{n} n}\\]" ], "text/latex": [ "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{4}{\\pi \\left(\\frac{1}{4}\\right)^{n} n}$$" ], "text/plain": [ "4/(pi*(1/4)^n*n)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "CP = [x==1, y==1, t == 1/4]\n", "print(\"Dominant asymptotic behaviour of the main diagonal is\")\n", "detM = getHes(H,[1,1,1],[x,y,t],CP).determinant().factor()\n", "ASM = leadingASM(G,H,detM,[1,1,1],[x,y,t],CP)\n", "show(ASM)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 9.4", "language": "sage", "name": "sagemath-9.4" }, "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.9.5" } }, "nbformat": 4, "nbformat_minor": 5 }