restart:
F := 1/(1-x-y-z-w+a*x*y*z*w);# This is the annihilating differential equation with a as a parameter
# (proven using creative telescoping)
DEQ := -w^2*(a^4*w^4+4*a^3*w^3+6*a^2*w^2+4*a*w-256*w+1)*(3*a*w-1)^2*diff(f(w),w,w,w)-3*w*(3*a*w-1)*(6*a^5*w^5+15*a^4*w^4+8*a^3*w^3-6*a^2*w^2-384*a*w^2-6*a*w+384*w-1)*diff(f(w),w,w)-(a*w+1)*(63*a^5*w^5-3*a^4*w^4-66*a^3*w^3+18*a^2*w^2+720*a*w^2+19*a*w-816*w+1)*diff(f(w),w)+(-9*a^6*w^5+3*a^5*w^4+6*a^4*w^3-18*a^3*w^2+360*a^2*w^2-13*a^2*w+384*a*w-a+24)*f(w);
DOP := -w^2*(a^4*w^4+4*a^3*w^3+6*a^2*w^2+4*a*w-256*w+1)*(3*a*w-1)^2*Dw^3-3*w*(3*a*w-1)*(6*a^5*w^5+15*a^4*w^4+8*a^3*w^3-6*a^2*w^2-384*a*w^2-6*a*w+384*w-1)*Dw^2-(a*w+1)*(63*a^5*w^5-3*a^4*w^4-66*a^3*w^3+18*a^2*w^2+720*a*w^2+19*a*w-816*w+1)*Dw-9*a^6*w^5+3*a^5*w^4+6*a^4*w^3-18*a^3*w^2+360*a^2*w^2-13*a^2*w+384*a*w-a+24:# Get initial values
start := add(w^k*coeftayl(F, [x, y, z, w] = [0, 0, 0, 0], [k,k,k,k]),k=0..6);# Check the GF is annihilated (up to order of power series truncation)
simplify(series(eval(subs(f(w)=start,DEQ)),w,7));# Convert to recurrence (with initial conditions)
with(gfun):
REC := {diffeqtorec(DEQ,f(w),u(n)), seq(u(k)=coeff(start,w,k),k=0..5)}:
map(factor,REC[1]);# Verify recurrence on later term in sequence
PR := rectoproc(REC,u(n)):
simplify(PR(10)) - coeftayl(F, [x,y,z,w]=[0,0,0,0], [10,10,10,10]);# So we have the right differential equation
# The exponential growth must be the reciprocal of a root of the following
degree(DOP,Dw);
lpol := coeff(DOP,Dw,3);# Note there is degeneration when a=27!
# The singularity at (3aw-1) = (81w-1), which would lead to an exponential
# growth of 81^n, becomes removable (it is a factor in every term of the differential operator)
map(collect,factor(subs(a=27,DOP)),Dw);# The exponential growth is the expected pair of conjugate algebraic number of modulus 9
map(factor,(81*w^4+14*w^3+w^2)*Dw^3+(486*w^3+63*w^2+3*w)*Dw^2+(567*w^2+48*w+1)*Dw+81*w+3);# We now find the exponential growth for generic a>30
# (and then check a=30)a := 'a';
H := a*w*x*y*z-w-x-y-z+1;
factor(Groebner[Basis]([H,seq(diff(H,x)*x-diff(H,k)*k,k=[y,z,w])],plex(x,y,z,w))):
r1,r2,r3,r4 := solve(%[1],w,explicit):
Groebner[Basis]([H,seq(x*diff(H,x)-k*diff(H,k),k=[y,z,w])],plex(x,y,z,w));# crit point defined by
pa := a*w^4-4*w+1;# No repeated roots unless a=27
discrim(pa,w);# At minus infinity and infinity the Sturm sequences are [sgn(a),-sgn(a),-1,sgn(27-a)] and [sgn(a),sgn(a),1,sgn(27-a)]
# At zero the Sturm sequence is [1,-1,-1,1]
# They have the same number of sign alternations when a>27, so there are never real roots
# When 0<a<27, there are 2 - 0 = 2 real roots
# When a<0, there are 3 - 1 = 2 real roots
p0 := pa;
p1 := diff(pa,w);
p2 := -rem(p0,p1,w);
p3 := -rem(p1,p2,w);# Calculate the first term of asymptotic contribution at pt (will be 0 if G vanishes there)
# Assume smooth minimal critical and diff(H,vars[-1])(pt) <> 0
smoothASM := proc(G,H,vars,pt)
local N,i,j,M,HES,CC,U,lambda,sbs:
# Get the Hessian determinant of the phase implicitly
N := nops(vars):
for i from 1 to N do for j from 1 to N do U[i,j] := vars[i]*vars[j]*diff(H,vars[i],vars[j]) od:od:
lambda := x*diff(H,x):
for i from 1 to N-1 do for j from 1 to N-1 do
if i <> j then M[i,j] := 1 + 1/lambda*(U[i,j]-U[i,N]-U[j,N]+U[N,N]):
else M[i,j] := 2 + 1/lambda*(U[i,i]-2*U[i,N]+U[N,N]):
fi:
od:od:
HES := LinearAlgebra[Determinant](Matrix([seq([seq(M[i,j],i=1..N-1)],j=1..N-1)]));
# Get leading constant
CC := simplify(-G/vars[-1]/diff(H,vars[-1]) * HES^(-1/2) * (2*Pi)^((1-N)/2));
sbs := seq(vars[j]=pt[j],j=1..N):
return evala(1/mul(j,j=pt))^n * n^((1-N)/2) * evala(subs(sbs, CC)):
end:alias(rr = RootOf(pa,w,index=1));ASM := simplify(smoothASM(1,H,[x,y,z,w],[rr,rr,rr,rr])) assuming a::real;###########################
# Specialize to the a=27 case
###########################H := 27*w*x*y*z-w-x-y-z+1;# Use recurrence of diagonal to get coefficients efficiently
with(gfun):
rec := (2+n)^3*u(n+2) + (3+2*n)*(17+21*n+7*n^2)*u(n+1) + 81*(1+n)^3*u(n);
PR := rectoproc({rec,u(0)=1,u(1)=-3},u(n)):rectodiffeq({rec,u(0)=1,u(1)=-3},u(n),y(w));# Check initial terms match (takes 10 seconds)
seq(PR(k),k=0..7);
seq(coeff(coeff(coeff(coeff(mtaylor(1/H,[x,y,z,w],30),x,k),y,k),z,k),w,k),k=0..7);# Critical points
factor(Groebner[Basis]([H,seq(diff(H,x)*x-diff(H,k)*k,k=[y,z,w])],plex(x,y,z,w)));# We treat the smooth points as if they were minimal
pw := 3*w^2+2*w+1;r1,r2 := solve(pw);ASM := smoothASM(1,H,[x,y,z,w],[r1,r1,r1,r1]) + smoothASM(1,H,[x,y,z,w],[r2,r2,r2,r2]);# Simplify algebraic constants
ASM2 := (11621/20420-I*sqrt(30803599)*(1/20420))*(-7-(4*I)*sqrt(2))^n/n^(3/2)/Pi^(3/2) + (11621/20420+I*sqrt(30803599)*(1/20420))*(-7+(4*I)*sqrt(2))^n/n^(3/2)/Pi^(3/2);evalf((11621/20420-I*sqrt(30803599)*(1/20420))/Pi^(3/2));# Numerically we have asymptotic constant (at n=1000)
evalf(PR(10000)/9^10000*(10000)^(3/2));# Our calculations give the following constant
evalf(subs(n=10000,ASM)/9^10000*10000^(3/2));# Alex Raichev's code gives the following contributions
AA1 := -3*sqrt(2)*81^r*(1/(I*sqrt(2)-1)^4)^r/(8*Pi^(3/2)*sqrt(-1/(I*sqrt(2)+4)^3)*r^(3/2)*(I*sqrt(2)+4)*(I*sqrt(2)-1)):
AA2 := -3*sqrt(2)*81^r*(1/(I*sqrt(2)+1)^4)^r/(8*Pi^(3/2)*sqrt(1/(I*sqrt(2)-4)^3)*r^(3/2)*(I*sqrt(2)-4)*(I*sqrt(2)+1)):
# The aymptotics at r=1000 matches our calculation
evalf(subs(r=10000,(AA1+AA2))/9^10000*10000^(3/2));# We seem to be off by a factor of 3
.3511685182 / .1170627761;# Sage package of Mezzarobba gives
ASMm := evalf((0.543449606382202 + 0.259547320313100*I)*(-7.000000000000000 +
5.656854249492380*I)^n/(sqrt(Pi)*n^(3/2)) + (0.543449606382202 -
0.259547320313100*I)*(-7.000000000000000 -
5.656854249492380*I)^n/(sqrt(Pi)*n^(3/2)));evalf(subs(n=10000,ASMm*n^(3/2)/9^n));# Find when coeffs are positive/negative
Pp := Z^4+16*Z^2-98;
cc := -3*sqrt(2)*rho*(rho^2+16)*(1/14) + 3*I*rho;fsolve(Pp,Z);
evalf(subs(rho=2.174378546,cc));
evalf(sqrt((2*I)*sqrt(2)-8)*(sqrt(2)+5*I));-Re(exp(I*theta)*cc) assuming rho>0,theta::real;3*cos(theta)*sqrt(2)*rho^3*(1/14)+24*cos(theta)*sqrt(2)*rho*(1/7)+3*sin(theta)*rho;SGN := cos(n*theta)*rho^2+16*cos(n*theta)+7*sqrt(2)*sin(n*theta);SGN2 := sqrt(16*rho^2+452) * sin(n*theta + arctan((rho^2+16)/(7*sqrt(2))));subs(theta=4*arctan(sqrt(2)),sin(4*n*arctan(sqrt(2))+arctan((1/14)*(rho^2+16)*sqrt(2))));