{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 7.6 (Distribution of Leaves in Planar Trees)\n", "We use Sturm's theorem to verify minimality of a critical point in a parametrized direction. \n", "*Requirements: None*" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{{\\left(u y z - y z - 2 \\, y + 1\\right)} y}{u y z - u z - y z - y + 1}$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{{\\left(u y z - y z - 2 \\, y + 1\\right)} y}{u y z - u z - y z - y + 1}\n", "\\end{math}" ], "text/plain": [ "(u*y*z - y*z - 2*y + 1)*y/(u*y*z - u*z - y*z - y + 1)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Recall from Example 3.11 in Chapter 3 that the number of binary trees on n nodes with\n", "# k leaves is the coefficient of u^k*y^n*z^n in the power series expansion of the following\n", "var('u,y,z,s,λ,t')\n", "F = (u*y*z - y*z - 2*y + 1)*y/(u*y*z - u*z - y*z - y + 1)\n", "G,H = F.numerator_denominator()\n", "r = [s,1,1]\n", "show(F)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The ideal generated by [a1*a3 - s, a1*s^2 - 2*a1*s + a1 - s^2, a2 - s, a3*s - s^2 + 2*s - 1, b1, b2, b3, lambdaR - s + 1, lambdaI] has (for fixed s) a single solution:\n" ] }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[a_{1} = \\frac{s^{2}}{s^{2} - 2 \\, s + 1}, a_{2} = s, a_{3} = \\frac{s^{2} - 2 \\, s + 1}{s}, b_{1} = 0, b_{2} = 0, b_{3} = 0, \\mathit{lambdaR} = s - 1, \\mathit{lambdaI} = 0\\right]$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left[a_{1} = \\frac{s^{2}}{s^{2} - 2 \\, s + 1}, a_{2} = s, a_{3} = \\frac{s^{2} - 2 \\, s + 1}{s}, b_{1} = 0, b_{2} = 0, b_{3} = 0, \\mathit{lambdaR} = s - 1, \\mathit{lambdaI} = 0\\right]\n", "\\end{math}" ], "text/plain": [ "[a1 == s^2/(s^2 - 2*s + 1),\n", " a2 == s,\n", " a3 == (s^2 - 2*s + 1)/s,\n", " b1 == 0,\n", " b2 == 0,\n", " b3 == 0,\n", " lambdaR == s - 1,\n", " lambdaI == 0]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# We split the variables into real/imaginary components and define the extended critical point system\n", "var('a1,a2,a3,b1,b2,b3,lambdaR,lambdaI')\n", "assume(a1,'real'); assume(a2,'real'); assume(a3,'real')\n", "assume(b1,'real'); assume(b2,'real'); assume(b3,'real')\n", "Hc = H.subs(u=a1+I*b1, y=a2+I*b2, z=a3+I*b3)\n", "Hr = Hc.real()\n", "Hi = Hc.imag()\n", "CPab1 = [diff(Hr,v)*v + diff(Hr,w)*w - R*lambdaR for (v,w,R) in zip([a1,a2,a3],[b1,b2,b3],r)]\n", "CPab2 = [diff(Hi,v)*v + diff(Hi,w)*w - R*lambdaI for (v,w,R) in zip([a1,a2,a3],[b1,b2,b3],r)]\n", "\n", "# Compute a Gröbner basis for the critical point equations in the a and b variables\n", "R = PolynomialRing(QQ,9,'a1,a2,a3,b1,b2,b3,lambdaR,lambdaI,s',order='lex')\n", "J = R.ideal([R(k) for k in [Hr,Hi] + CPab1 + CPab2])\n", "GB = R.ideal([R(k) for k in [Hr,Hi] + CPab1 + CPab2]).groebner_basis()\n", "\n", "# Note there are a finite number of (complex) solutions! \n", "[cp] = solve([SR(k) for k in GB],[a1,a2,a3,b1,b2,b3,lambdaR,lambdaI])\n", "print(\"The ideal generated by {} has (for fixed s) a single solution:\".format(GB))\n", "show(cp)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# To determine minimality, we add the extra equations using x and y variables\n", "# First, define xj + I*yj as a root of H\n", "# Then encode that the coordinate-wise moduli of xj + I*yj is a multiple of crit pt aj + I*bj\n", "# Finally, add the extra equations derived from critical point method\n", "var('x1,x2,x3,y1,y2,y3,ν')\n", "Hrxy = Hr.subs(a1==x1,a2==x2,a3==x3,b1==y1,b2==y2,b3==y3)\n", "Hixy = Hi.subs(a1==x1,a2==x2,a3==x3,b1==y1,b2==y2,b3==y3)\n", "mod = [X^2+Y^2 - t*(A^2+B^2) for (X,Y,A,B) in zip([x1,x2,x3],[y1,y2,y3],[a1,a2,a3],[b1,b2,b3])]\n", "nuxy = [(Y - ν*X)*diff(Hrxy,X) - (X + ν*Y)*diff(Hrxy,Y) for (X,Y) in zip([x1,x2,x3],[y1,y2,y3])]\n", "\n", "# Define the ideal with all these equations\n", "S = PolynomialRing(QQ,17,'a1,a2,a3,b1,b2,b3,x1,x2,x3,y1,y2,y3,lambdaR,lambdaI,ν,s,t',order='lex')\n", "Jxy = S.ideal([S(k) for k in list(GB) + [Hrxy,Hixy] + mod + nuxy])" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}(t - 1) \\cdot (s^{4} t^{2} - 3 s^{2} t + 2 s t - t + 1) \\cdot (s^{4} t^{3} - s^{4} t^{2} + 4 s^{3} t^{2} - 11 s^{2} t^{2} - s^{2} t + 6 s t^{2} + 6 s t - t^{2} - 2 t - 1) \\cdot (s^{4} t^{3} - s^{4} t^{2} + 4 s^{3} t^{2} - 3 s^{2} t^{2} - s^{2} t + 2 s t^{2} + 2 s t - t^{2} - 2 t - 1) \\cdot (s^{4} t^{3} - s^{4} t^{2} + 8 s^{3} t^{2} - 11 s^{2} t^{2} + 3 s^{2} t + 6 s t^{2} - 6 s t - t^{2} + 2 t - 1)$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}(t - 1) \\cdot (s^{4} t^{2} - 3 s^{2} t + 2 s t - t + 1) \\cdot (s^{4} t^{3} - s^{4} t^{2} + 4 s^{3} t^{2} - 11 s^{2} t^{2} - s^{2} t + 6 s t^{2} + 6 s t - t^{2} - 2 t - 1) \\cdot (s^{4} t^{3} - s^{4} t^{2} + 4 s^{3} t^{2} - 3 s^{2} t^{2} - s^{2} t + 2 s t^{2} + 2 s t - t^{2} - 2 t - 1) \\cdot (s^{4} t^{3} - s^{4} t^{2} + 8 s^{3} t^{2} - 11 s^{2} t^{2} + 3 s^{2} t + 6 s t^{2} - 6 s t - t^{2} + 2 t - 1)\n", "\\end{math}" ], "text/plain": [ "(t - 1) * (s^4*t^2 - 3*s^2*t + 2*s*t - t + 1) * (s^4*t^3 - s^4*t^2 + 4*s^3*t^2 - 11*s^2*t^2 - s^2*t + 6*s*t^2 + 6*s*t - t^2 - 2*t - 1) * (s^4*t^3 - s^4*t^2 + 4*s^3*t^2 - 3*s^2*t^2 - s^2*t + 2*s*t^2 + 2*s*t - t^2 - 2*t - 1) * (s^4*t^3 - s^4*t^2 + 8*s^3*t^2 - 11*s^2*t^2 + 3*s^2*t + 6*s*t^2 - 6*s*t - t^2 + 2*t - 1)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Compute the corresponding Gröbner basis. This can take a minute or so.\n", "# (Note: the same Gröbner computation is faster in Maple)\n", "GB2 = Jxy.groebner_basis()\n", "\n", "# The values of t at the solutions of the system satisfy the following\n", "tPol = GB2[-1].factor()\n", "show(tPol)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# To prove minimality for any s in (0,1) we must show that this polynomial\n", "# has no root with t in (0,1). We do this using Sturm sequences, examining\n", "# each factor -- other than t-1, which trivially has no root in (0,1) -- \n", "# seperately. Details of Sturm sequences are given in the textbook.\n", "\n", "# Function to compute a Sturm sequence\n", "def Sturm(poly):\n", " ST = []\n", " q = [poly, diff(poly,t)]\n", " while q[-1] != 0:\n", " _,r = q[-2].quo_rem(q[-1])\n", " q += [-r]\n", " return q[:-1]\n", "\n", "# Function to count sign alternations of sequence at a point\n", "def alternations(L,pt):\n", " alt = 0\n", " L2 = [p.subs(s=pt) for p in L if p.subs(s=pt) != 0]\n", " for k in range(len(L2)-1):\n", " if sign(L2[k]*L2[k+1])<0: alt += 1\n", " return alt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### First Factor" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This factor is s^4*t^2 + (-3*s^2 + 2*s - 1)*t + 1\n" ] }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|Sturm|\\phantom{\\verb!x!}\\verb|sequence|\\phantom{\\verb!x!}\\verb|at|\\phantom{\\verb!x!}\\verb|t=0|\\phantom{\\verb!x!}\\verb|is| \\left[1, -3 s^{2} + 2 s - 1, \\frac{\\frac{5}{4} s^{4} - 3 s^{3} + \\frac{5}{2} s^{2} - s + \\frac{1}{4}}{s^{4}}\\right]$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|Sturm|\\phantom{\\verb!x!}\\verb|sequence|\\phantom{\\verb!x!}\\verb|at|\\phantom{\\verb!x!}\\verb|t=0|\\phantom{\\verb!x!}\\verb|is| \\left[1, -3 s^{2} + 2 s - 1, \\frac{\\frac{5}{4} s^{4} - 3 s^{3} + \\frac{5}{2} s^{2} - s + \\frac{1}{4}}{s^{4}}\\right]\n", "\\end{math}" ], "text/plain": [ "'Sturm sequence at t=0 is' [1, -3*s^2 + 2*s - 1, (5/4*s^4 - 3*s^3 + 5/2*s^2 - s + 1/4)/s^4]" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|Sturm|\\phantom{\\verb!x!}\\verb|sequence|\\phantom{\\verb!x!}\\verb|at|\\phantom{\\verb!x!}\\verb|t=1|\\phantom{\\verb!x!}\\verb|is| \\left[s^{4} - 3 s^{2} + 2 s, 2 s^{4} - 3 s^{2} + 2 s - 1, \\frac{\\frac{5}{4} s^{4} - 3 s^{3} + \\frac{5}{2} s^{2} - s + \\frac{1}{4}}{s^{4}}\\right]$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\verb|Sturm|\\phantom{\\verb!x!}\\verb|sequence|\\phantom{\\verb!x!}\\verb|at|\\phantom{\\verb!x!}\\verb|t=1|\\phantom{\\verb!x!}\\verb|is| \\left[s^{4} - 3 s^{2} + 2 s, 2 s^{4} - 3 s^{2} + 2 s - 1, \\frac{\\frac{5}{4} s^{4} - 3 s^{3} + \\frac{5}{2} s^{2} - s + \\frac{1}{4}}{s^{4}}\\right]\n", "\\end{math}" ], "text/plain": [ "'Sturm sequence at t=1 is' [s^4 - 3*s^2 + 2*s,\n", " 2*s^4 - 3*s^2 + 2*s - 1,\n", " (5/4*s^4 - 3*s^3 + 5/2*s^2 - s + 1/4)/s^4]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# First factor \n", "p1 = QQ[s].fraction_field()[t](tPol[1][0])\n", "q = Sturm(p1)\n", "C0 = [k.subs(t=0) for k in q]\n", "C1 = [k.subs(t=1) for k in q]\n", "print(\"This factor is\",p1)\n", "show(\"Sturm sequence at t=0 is\", C0)\n", "show(\"Sturm sequence at t=1 is\", C1)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# No roots of a numerator or denominator for either Sturm sequence lie in (0,1)\n", "rts = flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C0])\n", "rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C0])\n", "rts += flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C1])\n", "rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C1])\n", "rts01 = list(filter(lambda x: x>0 and x<1, rts))\n", "rts01" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Thus, the signs are the same for all s in (0,1), including both endpoints \n", "# This means the number of roots of this factor with t in (0,1] for s in (0,1) is\n", "alternations(C1,1/2) - alternations(C0,1/2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Second Factor" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This factor is s^4*t^3 + (-s^4 + 4*s^3 - 11*s^2 + 6*s - 1)*t^2 + (-s^2 + 6*s - 2)*t - 1\n" ] } ], "source": [ "p2 = QQ[s].fraction_field()[t](tPol[2][0])\n", "q = Sturm(p2)\n", "C0 = [k.subs(t=0) for k in q]\n", "C1 = [k.subs(t=1) for k in q]\n", "print(\"This factor is\",p2)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0.3542486889354094?, 0.3851648071345041?, 0.3941724074329110?]" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The sign conditions of the Sturm sequence can only change at roots of the numerator\n", "# or denominator of its entries -- i.e., for s in (0,1) it can only change at these values\n", "rts = flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C0])\n", "rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C0])\n", "rts += flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C1])\n", "rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C1])\n", "rts01 = list(filter(lambda x: x>0 and x<1, rts))\n", "rts01 = sorted(set(rts01))\n", "rts01" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0, 0, 0, 0]" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# For s in (0,1) between these points the sign conditions are constant,\n", "# so we can test the sign conditions by picking values of s between each point\n", "# This shows there are no roots for t in (0,1) and s in (0,1) when s not in rts01\n", "pts = [0.1, 0.36, 0.39, 0.5]\n", "[alternations(C1,pt) - alternations(C0,pt) for pt in pts]" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[[7.140984766613703?], [5.388399446915442?], [5.018578328500075?]]" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Manually show there are no roots for t in (0,1) when s in rts01 by solving the poly\n", "[AA[t](SR(p2).subs(s=rt)).roots(multiplicities=False) for rt in rts01]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Third Factor" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This factor is s^4*t^3 + (-s^4 + 4*s^3 - 3*s^2 + 2*s - 1)*t^2 + (-s^2 + 2*s - 2)*t - 1\n" ] } ], "source": [ "p3 = QQ[s].fraction_field()[t](tPol[3][0])\n", "q = Sturm(p3)\n", "C0 = [k.subs(t=0) for k in q]\n", "C1 = [k.subs(t=1) for k in q]\n", "print(\"This factor is\",p3)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0.7254470427491280?]" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Find roots in (0,1)\n", "rts = flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C0])\n", "rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C0])\n", "rts += flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C1])\n", "rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C1])\n", "rts01 = list(filter(lambda x: x>0 and x<1, rts))\n", "rts01 = sorted(set(rts01))\n", "rts01" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0, 0]" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# For s in (0,1) on each side of this point the sign conditions are constant,\n", "# so we can test the sign conditions by picking two values of s in (0,1)\n", "# This shows there are no roots for t in (0,1) and s in (0,1) when s not the element of rts01\n", "pts = [0.1,0.9]\n", "[alternations(C1,pt) - alternations(C0,pt) for pt in pts]" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[[2.148303747649492?]]" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Manually show there are no roots for t in (0,1) when s in rts01 by solving the poly\n", "[AA[t](SR(p3).subs(s=rt)).roots(multiplicities=False) for rt in rts01]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Fourth (Final) Factor" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This factor is s^4*t^3 + (-s^4 + 8*s^3 - 11*s^2 + 6*s - 1)*t^2 + (3*s^2 - 6*s + 2)*t - 1\n" ] } ], "source": [ "# Final factor \n", "p4 = QQ[s].fraction_field()[t](tPol[4][0])\n", "q = Sturm(p4)\n", "C0 = [k.subs(t=0) for k in q]\n", "C1 = [k.subs(t=1) for k in q]\n", "print(\"This factor is\",p4)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0.2370545200314687?,\n", " 0.2403702750125248?,\n", " 0.2548844570965663?,\n", " 0.2616705580491203?,\n", " 0.3720842206893895?,\n", " 0.4226497308103743?]" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Find roots in (0,1)\n", "rts = flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C0])\n", "rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C0])\n", "rts += flatten([QQ[s](SR(k).numerator()).roots(AA, multiplicities=False) for k in C1])\n", "rts += flatten([QQ[s](SR(k).denominator()).roots(AA, multiplicities=False) for k in C1])\n", "rts01 = list(filter(lambda x: x>0 and x<1, rts))\n", "rts01 = sorted(set(rts01))\n", "rts01" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[0, 0, 0, 0, 0, 0, 0]" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# For s in (0,1) on each side of this point the sign conditions are constant,\n", "# so we can test the sign conditions by picking two values of s in (0,1)\n", "# This shows there are no roots for t in (0,1) and s in (0,1) when s not the element of rts01\n", "pts = [0.1, 0.24, 0.25, 0.26, 0.3, 0.4, 0.9]\n", "[alternations(C1,pt) - alternations(C0,pt) for pt in pts]" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[[1.6630008475?, 13.799307214583437?],\n", " [1.6745190182?],\n", " [1.725018869517811?],\n", " [1.748599118643231?],\n", " [2.084201751623?],\n", " [2.170435189413383?]]" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Manually show there are no roots for t in (0,1) when s in rts01 by solving the poly\n", "[AA[t](SR(p4).subs(s=rt)).roots(multiplicities=False) for rt in rts01]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Computing Asymptotics\n", "\n", "Combining everything above, we have shown that the single critical point is minimal for all $s \\in (0,1)$. We end by computing asymptotics after importing our code for smooth ACSV from Chapter 5. If you just want to compute asymptotics *assuming* minimality, you only need to run the code in this section." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "# Import code for smooth asymptotics\n", "# 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": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{{\\left(u y z - y z - 2 \\, y + 1\\right)} y}{u y z - u z - y z - y + 1}$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{{\\left(u y z - y z - 2 \\, y + 1\\right)} y}{u y z - u z - y z - y + 1}\n", "\\end{math}" ], "text/plain": [ "(u*y*z - y*z - 2*y + 1)*y/(u*y*z - u*z - y*z - y + 1)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Restate basic parameters to make this section independent of above\n", "var('u,y,z,s')\n", "F = (u*y*z - y*z - 2*y + 1)*y/(u*y*z - u*z - y*z - y + 1)\n", "show(F)\n", "G,H = F.numerator_denominator()\n", "r = [s,1,1]\n", "vars = [u,y,z]" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The asymptotic expansion for the [s,1,1] diagonal begins\n" ] }, { "data": { "text/html": [ " $\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\left(\\frac{{\\left(s - 1\\right)}^{2 \\, s}}{{\\left(s^{2} - 2 \\, s + 1\\right)} s^{2 \\, s}}\\right)^{n}}{n} \\Bigg( \\frac{1}{2 \\, \\pi n} + O\\Bigg( \\frac{1}{n^{2}} \\Bigg)\\Bigg)$ " ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\frac{\\left(\\frac{{\\left(s - 1\\right)}^{2 \\, s}}{{\\left(s^{2} - 2 \\, s + 1\\right)} s^{2 \\, s}}\\right)^{n}}{n} \\Bigg( \\frac{1}{2 \\, \\pi n} + O\\Bigg( \\frac{1}{n^{2}} \\Bigg)\\Bigg)\n", "\\end{math}" ], "text/plain": [ "((s - 1)^(2*s)/((s^2 - 2*s + 1)*s^(2*s)))^n/n \\Bigg( 1/2/(pi*n) + O\\Bigg( n^(-2) \\Bigg)\\Bigg)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Get critical point (shown minimal above)\n", "[CP] = critpt(H,r,vars)\n", "\n", "# Get and print asymptotics\n", "M = 2\n", "g = vars[-1].subs(solve(H,vars[-1]))\n", "ex,pw,se = smoothContrib(G,H,r,vars,CP,M,g)\n", "print(\"The asymptotic expansion for the [s,1,1] diagonal begins\")\n", "disp_asm(ex,pw,se,M)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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 }