{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 7.3 (Apéry Revisited)\n", "We run our algorithm and compute the Kronecker representation \"by hand\" using Sage's capability to compute Gröbner bases. Note that Example 7.1 uses the Maple package of Melczer and Salvy to automatically compute this.\n", "\n", "*Requirements: None*" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# We will need our procedure to get the Hessian matrix appearing in asymptotics\n", "# This code copied from Chapter 5\n", "def getHes(H,vars):\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] = 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" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[-(x*y*z + y*z + y + z + 1)*w*(x + 1)*(y + 1)*(z + 1) + 1,\n", " -(x*y*z + y*z + y + z + 1)*w*(x + 1)*(y + 1)*(z + 1) - λ,\n", " -(w*(x + 1)*(y + 1)*y*(z + 1)*z + (x*y*z + y*z + y + z + 1)*w*(y + 1)*(z + 1))*x - λ,\n", " -((x*z + z + 1)*w*(x + 1)*(y + 1)*(z + 1) + (x*y*z + y*z + y + z + 1)*w*(x + 1)*(z + 1))*y - λ,\n", " -((x*y + y + 1)*w*(x + 1)*(y + 1)*(z + 1) + (x*y*z + y*z + y + z + 1)*w*(x + 1)*(y + 1))*z - λ,\n", " -(t^3*x*y*z + t^2*y*z + t*y + t*z + 1)*(t*x + 1)*(t*y + 1)*(t*z + 1)*t*w + 1]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Set basic info for the Apéry example\n", "var('w,x,y,z,λ,t,u')\n", "vars = [w,x,y,z,λ,t]\n", "H = 1 - w*(1 + x)*(1 + y)*(1 + z)*(x*y*z + y*z + y + z + 1)\n", "\n", "# Form the extended critical point system\n", "C = [H] + [v * diff(H,v) - λ for v in vars[0:-2]]\n", "S = C + [H.subs([v==t*v for v in vars])]\n", "S" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "u^14 + 1144*u^13 + 561458*u^12 + 153334416*u^11 + 25196311153*u^10 + 2498103299744*u^9 + 139498034620672*u^8 + 3517543585569916*u^7 + 11318250382599528*u^6 + 16555760259324420*u^5 + 19429225239935936*u^4 + 19058949700404944*u^3 + 9535423903452608*u^2 - 2094865315255056*u - 85832780580609892" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Define the multivariate polynomial ring and compute a Gröbner basis\n", "R = PolynomialRing(QQbar,7,'w,x,y,z,λ,t,u',order='lex')\n", "S = S + [u - (w+t)]\n", "GB = R.ideal([R(k) for k in S]).groebner_basis()\n", "\n", "# The Kronecker is representation parametrized by the roots of the following polynomial\n", "P = QQ[u](GB[-1])\n", "P" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The polynomial P is square-free: its gcd with its derivative is zero\n", "# (Note: we need to cast between polynomial and symbolic rings to compute desired quantities)\n", "QQ[u](P).gcd(QQ[u](SR(P).diff(u)))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}{\\mathbf{#1}}\\verb|In|\\phantom{\\verb!x!}\\verb|the|\\phantom{\\verb!x!}\\verb|Kronecker|\\phantom{\\verb!x!}\\verb|rep.|\\phantom{\\verb!x!}\\verb|the|\\phantom{\\verb!x!}\\verb|variable|\\phantom{\\verb!x!}\\verb|\t| w \\phantom{\\verb!x!}\\verb|is|\\phantom{\\verb!x!}\\verb|represented|\\phantom{\\verb!x!}\\verb|by\t| -1148 u^{13} - 1126240 u^{12} - 461022128 u^{11} - 100896914288 u^{10} - 12481652051316 u^{9} - 833437695269152 u^{8} - 24278547801823360 u^{7} - 78457363285316984 u^{6} - 114190314674867744 u^{5} - 134092823841986328 u^{4} - 132015613891611904 u^{3} - 65877266035071264 u^{2} + 11516697886625280 u + 604881444905678176 \\phantom{\\verb!x!}\\verb|divided|\\phantom{\\verb!x!}\\verb|by\t| 14 u^{13} + 14872 u^{12} + 6737496 u^{11} + 1686678576 u^{10} + 251963111530 u^{9} + 22482929697696 u^{8} + 1115984276965376 u^{7} + 24622805098989412 u^{6} + 67909502295597168 u^{5} + 82778801296622100 u^{4} + 77716900959743744 u^{3} + 57176849101214832 u^{2} + 19070847806905216 u - 2094865315255056\n", "\\end{math}" ], "text/plain": [ "'In the Kronecker rep. the variable \\t' w ' is represented by\\t' -1148*u^13 - 1126240*u^12 - 461022128*u^11 - 100896914288*u^10 - 12481652051316*u^9 - 833437695269152*u^8 - 24278547801823360*u^7 - 78457363285316984*u^6 - 114190314674867744*u^5 - 134092823841986328*u^4 - 132015613891611904*u^3 - 65877266035071264*u^2 + 11516697886625280*u + 604881444905678176 ' divided by\\t' 14*u^13 + 14872*u^12 + 6737496*u^11 + 1686678576*u^10 + 251963111530*u^9 + 22482929697696*u^8 + 1115984276965376*u^7 + 24622805098989412*u^6 + 67909502295597168*u^5 + 82778801296622100*u^4 + 77716900959743744*u^3 + 57176849101214832*u^2 + 19070847806905216*u - 2094865315255056" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Compute the Q polynomials for the Kronecker representation from the Gröbner basis\n", "Pd = diff(P,u)\n", "Q = {}\n", "for v in vars:\n", " [Qv] = list(filter(lambda p: SR(p).has(v), GB))\n", " _,Q[v] = QQ[u](Pd * v.subs(solve(SR(Qv),v))).quo_rem(P)\n", "sys = [SR(Q[v]/Pd) for v in vars]\n", "\n", "# For instance, here is the first variable in the Kronecker representation\n", "show(\"In the Kronecker rep. the variable \\t\",vars, \" is represented by\\t\",Q[vars],\" divided by\\t\",Pd)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[13, 13, 12, 12, 13, 13]" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# All polynomials in the Kronecker representation have degree close to degree(P) = 14\n", "[Q[v].degree() for v in vars]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[18, 17, 17, 17, 17, 18]" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# They have coefficients that are at most 18 digits long\n", "# As expected, this is an *order of magnitude* smaller than the constants\n", "# appearing in the Gröbner basis, which have constants with ~100 digits\n", "[Q[v].norm(infinity).log(10).round() for v in vars]" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "99" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "max(map(lambda p: abs(QQ(p).numerator()),GB.coefficients())).log(10).round()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(u^2 + 162*u - 167) * (u^12 + 982*u^11 + 402541*u^10 + 88286768*u^9 + 10961079084*u^8 + 737152378392*u^7 + 21909849528196*u^6 + 91252409193628*u^5 + 194304964440524*u^4 + 317508355295408*u^3 + 441800743647348*u^2 + 511124563867704*u + 513968745991676)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Note that P factors\n", "P.factor()" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "u^2 + 162*u - 167" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The critical points, where t=1, are given by the values of u which are roots of \n", "Pt = gcd(P,Pd-Q[t])\n", "Pt" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "There is a critical point where the coordinates [w, x, y, z] equal\n", "[-164.0243866176395?, -0.4142135623730951?, -0.7071067811865475?, -0.7071067811865475?]\n", "There is a critical point where the coordinates [w, x, y, z] equal\n", "[0.02438661763951283?, 2.414213562373095?, 0.7071067811865475?, 0.7071067811865475?]\n" ] } ], "source": [ "# There is a critical point with positive coordinates (and one with negative coordinates)\n", "for uval in Pt.roots(QQbar, multiplicities=False):\n", " print(\"There is a critical point where the coordinates {} equal\".format(vars[:-2]))\n", " print([p.subs(u == uval) for p in sys[:-2]])" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "At the solutions of the system, t takes the values [-0.005979808060975540?, 1.000000000000000?, 1.139687851637440?, 1.580884722469457?, 2.409071027725915?, 1.000000000000000?, 0.7665952278602726? - 2.211004338122884?*I, 0.7665952278602726? + 2.211004338122884?*I, -2.040073239787960? - 0.6346682615488171?*I, -2.040073239787960? + 0.6346682615488171?*I, -0.7916404858532640? - 1.499007434480210?*I, -0.7916404858532640? + 1.499007434480210?*I, 0.5032866008950336? - 1.375523098438217?*I, 0.5032866008950336? + 1.375523098438217?*I]\n", "Here are the values lying in (0,1): []\n" ] } ], "source": [ "# We can compute the algebraic values of t at solutions of the system\n", "tvals = [sys[-1].subs(u == uval) for uval in P.roots(QQbar, multiplicities=False)]\n", "print(\"At the solutions of the system, t takes the values {}\".format(tvals))\n", "\n", "# And we can verify none lie in (0,1)\n", "# This means the point with positive coordinates is *minimal*\n", "num_in_zo = list(filter(lambda c: c.is_real() and c>0 and c<1, tvals))\n", "print(\"Here are the values lying in (0,1): {}\".format(num_in_zo))" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The Kronecker representation for the critical points (where t=1 in the extended system) is given by\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}{\\mathbf{#1}}\\left[w = -\\frac{2 \\, {\\left(41 \\, u - 43\\right)}}{u + 81}, x = \\frac{u + 197}{u + 81}, y = \\frac{58}{u + 81}, z = \\frac{58}{u + 81}\\right]\n", "\\end{math}" ], "text/plain": [ "[w == -2*(41*u - 43)/(u + 81),\n", " x == (u + 197)/(u + 81),\n", " y == 58/(u + 81),\n", " z == 58/(u + 81)]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "where u satisfies u^2 + 162*u - 167=0.\n" ] } ], "source": [ "# We can simplify the system to the points where t = 1\n", "# First, compute reciprocal of the new derivative modulo the new P\n", "newP = Pt\n", "newPd = diff(newP,u)\n", "_,invPd,_ = Pd.xgcd(newP)\n", "\n", "# Then update the Q polynomials using polynomial division\n", "newQ = {}\n", "for v in vars[0:-1]:\n", " _,newQ[v] = QQ[u](Q[v] * newPd * invPd).quo_rem(newP)\n", "\n", "# Print the new representation\n", "print(\"The Kronecker representation for the critical points (where t=1 in the extended system) is given by\")\n", "show([v == newQ[v]/newPd for v in vars[0:-2]])\n", "print(\"where u satisfies {}=0.\".format(newP))" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The Kronecker representation for the critical points (where t=1 in the extended system) is given by\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}{\\mathbf{#1}}\\left[T = \\frac{17 \\, u - 15}{u + 81}, h = \\frac{8 \\, {\\left(6 \\, u - 7\\right)}}{u + 81}, g = \\left(-1\\right)\\right]\n", "\\end{math}" ], "text/plain": [ "[T == (17*u - 15)/(u + 81), h == 8*(6*u - 7)/(u + 81), g == -1]" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "where u satisfies u^2 + 162*u - 167=0.\n" ] } ], "source": [ "# Compute the Hessian determinant Htilde for final asymptotic formula\n", "hes = vars[-3]*diff(H,vars[-3])*getHes(H,vars[0:-2])\n", "hesdet = hes.determinant().factor()\n", "\n", "# Add the Hessian, the product of the variables, and the negation of the numerator\n", "# to the Kronecker rep. We use a straightforward but inefficient approach by \n", "# computing another Gröbner basis!\n", "var('T,h,g')\n", "PolySys = [newP, T - w*x*y*z, h - hesdet, g+1] + S\n", "R2 = PolynomialRing(QQ,10,'g,T,h,w,x,y,z,λ,t,u',order='lex')\n", "GB2 = R2.ideal([R2(str(k)) for k in PolySys]).groebner_basis()\n", "for v in [T,h,g]:\n", " [Qv] = list(filter(lambda p: SR(p).has(v), GB2))\n", " _,newQ[v] = QQ[u](newPd * v.subs(solve(SR(Qv),v))).quo_rem(newP)\n", " \n", "# Print the new elements of the representation\n", "print(\"The Kronecker representation for the critical points (where t=1 in the extended system) is given by\")\n", "show([v == newQ[v]/newPd for v in [T,h,g]])\n", "print(\"where u satisfies {}=0.\".format(newP))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Asymptotics are given by evaluating\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}{\\mathbf{#1}}\\frac{\\left(\\frac{u + 81}{17 \\, u - 15}\\right)^{n} \\sqrt{-\\frac{u + 81}{6 \\, u - 7}}}{8 \\, \\left(\\pi n\\right)^{\\frac{3}{2}}}\n", "\\end{math}" ], "text/plain": [ "1/8*((u + 81)/(17*u - 15))^n*sqrt(-(u + 81)/(6*u - 7))/(pi*n)^(3/2)" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "At the root 1.024386617639513? of the polynomial u^2 + 162*u - 167.\n", "This gives dominant asymptotic behaviour\n" ] }, { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}{\\mathbf{#1}}\\frac{\\sqrt{96.08326112068523?} 33.97056274847714?^{n}}{8 \\, \\left(\\pi n\\right)^{\\frac{3}{2}}}\n", "\\end{math}" ], "text/plain": [ "1/8*sqrt(96.08326112068523?)*33.97056274847714?^n/(pi*n)^(3/2)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Get final asymptotic constants\n", "d = 4\n", "A = newQ[g]/newQ[λ]\n", "B = newQ[λ]^(d-1)/newQ[h]*newPd^(2-d)\n", "C = newPd/newQ[T]\n", "\n", "# Find root of P that gives critical point with positive coordinates\n", "[posu] = filter(lambda uval: sys.subs(u == uval)>0, newP.roots(QQbar, multiplicities=False))\n", "\n", "# Print result of algorithm\n", "var('n')\n", "ASM = ((n*2*pi)^(-3/2) * A * sqrt(B) * SR(C)^n).simplify()\n", "print(\"Asymptotics are given by evaluating\")\n", "show(ASM)\n", "print(\"At the root {} of the polynomial {}.\".format(posu,newP))\n", "print(\"This gives dominant asymptotic behaviour\")\n", "show(ASM.subs(u=posu))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "" ], "text/latex": [ "\\begin{math}\n", "\\newcommand{\\Bold}{\\mathbf{#1}}\\frac{\\sqrt{2} \\sqrt{\\pi n} {\\left(12 \\, \\sqrt{2} + 17\\right)}^{n} \\sqrt{17 \\, \\sqrt{2} + 24}}{8 \\, \\pi^{2} n^{2}}\n", "\\end{math}" ], "text/plain": [ "1/8*sqrt(2)*sqrt(pi*n)*(12*sqrt(2) + 17)^n*sqrt(17*sqrt(2) + 24)/(pi^2*n^2)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Because u is the root of a quadratic polynomial, we can express asymptotics in radicals\n", "maxima_calculus('algebraic: true;') # This helps simplify roots\n", "show(ASM.subs(u=posu.radical_expression()).simplify_full())" ] }, { "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 }