restart:
F := 1/(1-x-y-z-w+a*(x*y*z+w*x*y+w*x*z+w*y*z));normal(subs(w=w/x/y/z,F)/x/y/z);DEQ := 16*a^2*(552960*a^19*w^10+1399680*a^18*w^9+4952064*a^17*w^9+1636992*a^16*w^8-25116672*a^15*w^8+26784*a^15*w^7-5417472*a^14*w^7-139968*a^14*w^6+34619392*a^13*w^7+1082592*a^13*w^6-4632064*a^12*w^6-197496*a^12*w^5-3923968*a^11*w^6-474336*a^11*w^5-9180*a^11*w^4+17260800*a^10*w^5+489048*a^10*w^4-32710656*a^9*w^5+478560*a^9*w^4+26730*a^9*w^3-17890560*a^8*w^4-66736*a^8*w^3+33013760*a^7*w^4+522*a^8*w^2-1902624*a^7*w^3+10614*a^7*w^2+10432000*a^6*w^3-226448*a^6*w^2-13942784*a^5*w^3+555*a^6*w+1338016*a^5*w^2-12010*a^5*w-3245568*a^4*w^2+6*a^5+91048*a^4*w+2813952*a^3*w^2-220*a^4-310176*a^3*w+2342*a^3+475008*a^2*w-10888*a^2-258048*a*w+23328*a-18816)*f(w)
+(115015680*a^21*w^11+229201920*a^20*w^10+660996096*a^19*w^10+6770688*a^18*w^9-4600430592*a^17*w^9+5465088*a^17*w^8-1089472512*a^16*w^8-13450752*a^16*w^7+8701607936*a^15*w^8+74772480*a^15*w^7+466993152*a^14*w^7-7384320*a^14*w^6-5421793280*a^13*w^7-182454272*a^13*w^6+177408*a^13*w^5+2401054720*a^12*w^6+46640640*a^12*w^5-3036413952*a^11*w^6+1728*a^12*w^4+323813376*a^11*w^5+5059776*a^11*w^4-3998883840*a^10*w^5-26426624*a^10*w^4+6678249472*a^9*w^5+144144*a^10*w^3-447541248*a^9*w^4-1996352*a^9*w^3+2923573248*a^8*w^4+72*a^9*w^2-25014784*a^8*w^3-4170711040*a^7*w^4+13320*a^8*w^2+333819904*a^7*w^3-2053056*a^7*w^2-1154105344*a^6*w^3+204*a^7*w+25030912*a^6*w^2+1240399872*a^5*w^3-75376*a^6*w-116103168*a^5*w^2+12*a^6+1118784*a^5*w+235763712*a^4*w^2-1244*a^5-6454784*a^4*w-175767552*a^3*w^2+20248*a^4+17522688*a^3*w-141824*a^3-21540864*a^2*w+506112*a^2+8847360*a*w-906240*a+645120)*diff(f(w),w)
+(175767552*a^21*w^12+302469120*a^20*w^11+648806400*a^19*w^11-325312512*a^18*w^10-6402473984*a^17*w^10-24514560*a^17*w^9-1332391936*a^16*w^9-12572928*a^16*w^8+15289286656*a^15*w^9+30898176*a^15*w^8+2209120256*a^14*w^8-131328*a^14*w^7-15039463424*a^13*w^8-150333440*a^13*w^7+1476288*a^13*w^6+1238999040*a^12*w^7+30110976*a^12*w^6+2328363008*a^11*w^7+8640*a^12*w^5+596514816*a^11*w^6+5944128*a^11*w^5-5632180224*a^10*w^6-19426048*a^10*w^5+8337489920*a^9*w^6+184464*a^10*w^4-970385408*a^9*w^5-6126464*a^9*w^4+5805187072*a^8*w^5+1080*a^9*w^3-43254528*a^8*w^4-8022654976*a^7*w^5-274656*a^8*w^3+791603200*a^7*w^4-2024064*a^7*w^3-2951577600*a^6*w^4-5916*a^7*w^2+57041152*a^6*w^3+3279028224*a^5*w^4-122416*a^6*w^2-332052480*a^5*w^3-12*a^6*w+2842304*a^5*w^2+776015872*a^4*w^3-4476*a^5*w-19492096*a^4*w^2-647233536*a^3*w^3+79680*a^4*w+61704192*a^3*w^2-20*a^4-559296*a^3*w-92725248*a^2*w^2+288*a^3+1960704*a^2*w+53084160*a*w^2-1536*a^2-3434496*a*w+3584*a+2396160*w-3072)*diff(f(w),w$2)
+(w*(66060288*a^21*w^12+101744640*a^20*w^11+136314880*a^19*w^11-222756864*a^18*w^10-2186280960*a^17*w^10-25509888*a^17*w^9-315326464*a^16*w^9-3345408*a^16*w^8+6366953472*a^15*w^9+18456576*a^15*w^8+1243873280*a^14*w^8+3269376*a^14*w^7-8206155776*a^13*w^8+12132352*a^13*w^7+833472*a^13*w^6-673644544*a^12*w^7+4225792*a^12*w^6+4139778048*a^11*w^7+6048*a^12*w^5+155656192*a^11*w^6+1444800*a^11*w^5-1397882880*a^10*w^6-1073408*a^10*w^5+1541406720*a^9*w^6+15120*a^10*w^4-421650432*a^9*w^5-2665728*a^9*w^4+2416246784*a^8*w^5-24277248*a^8*w^4-3229614080*a^7*w^5-246608*a^8*w^3+424599552*a^7*w^4-248448*a^7*w^3-1588330496*a^6*w^4-6748*a^7*w^2+30968064*a^6*w^3+1764753408*a^5*w^4-42000*a^6*w^2-201891840*a^5*w^3-48*a^6*w+1632192*a^5*w^2+499548160*a^4*w^3-2764*a^5*w-12103424*a^4*w^2-434110464*a^3*w^3+54064*a^4*w+40394752*a^3*w^2-25*a^4-378176*a^3*w-64946176*a^2*w^2+360*a^3+1291520*a^2*w+40894464*a*w^2-1920*a^2-2187264*a*w+4480*a+1474560*w-3840))*diff(f(w),w$3)
+(w^2*(4*a^4*w^2-8*a^2*w-a+4)*(256*a^4*w^3-512*a^2*w^2-64*a*w+256*w-1)*(6144*a^13*w^7+8640*a^12*w^6+28672*a^11*w^6+8832*a^10*w^5-104448*a^9*w^5+240*a^9*w^4-16064*a^8*w^4-216*a^8*w^3+73728*a^7*w^4+1888*a^7*w^3-17664*a^6*w^3-444*a^6*w^2+26624*a^5*w^3-1344*a^5*w^2-12*a^5*w+21568*a^4*w^2-64*a^4*w-36864*a^3*w^2+1440*a^3*w-5*a^3-5504*a^2*w+52*a^2+6144*a*w-176*a+192))*diff(f(w),w$4):degree(w^2*(4*a^4*w^2-8*a^2*w-a+4)*(256*a^4*w^3-512*a^2*w^2-64*a*w+256*w-1)*(6144*a^13*w^7+8640*a^12*w^6+28672*a^11*w^6+8832*a^10*w^5-104448*a^9*w^5+240*a^9*w^4-16064*a^8*w^4-216*a^8*w^3+73728*a^7*w^4+1888*a^7*w^3-17664*a^6*w^3-444*a^6*w^2+26624*a^5*w^3-1344*a^5*w^2-12*a^5*w+21568*a^4*w^2-64*a^4*w-36864*a^3*w^2+1440*a^3*w-5*a^3-5504*a^2*w+52*a^2+6144*a*w-176*a+192),w);# Get initial values
start := add(w^k*coeftayl(F, [x, y, z, w] = [0, 0, 0, 0], [k,k,k,k]),k=0..7):# Check the GF is annihilated (up to order of power series truncation)
simplify(series(eval(subs(f(w)=start,DEQ)),w,6));# 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]):# So we have the right differential equation
# The exponential growth must be the reciprocal of a root of the following
lpol := coeff(DEQ,diff(f(w),w$4),1);# Again at a=64/27 there is degeneration, and leading poly is
w^2*(256*w-81)*(16777216*w^2-5971968*w+216513)*(8388608*w^2+12773376*w-2381643);
fsolve(%);# Only singular when a = 64/27
a := 'a';
H := 1-x-y-z-w+a*(x*y*z+w*x*y+w*x*z+w*y*z);
Groebner[Basis]([seq(diff(H,k),k=[x,y,z,w]),H],plex(a,x,y,z,w));# By GWS, there is a minimal critical point (smooth unless a=64/27) given by the following
factor(Groebner[Basis]([z-w,z-x,z-y,H,seq(diff(H,x)*x-diff(H,k)*k,k=[y,z,w])],plex(x,y,z,w,a)));
pa := %[1];
r1,r2,r3 := solve(pa,w,explicit):# No repeated roots unless a=0 or 64/27
discrim(pa,w);# At minus infinity and infinity the Sturm sequences are [-sgn(a),sgn(a),-1,sgn(64-27a)] and [sgn(a),sgn(a),1,sgn(64-27a)]
# At 0 the Sturm sequence is [1,-1,-1,sgn(64-27a)]
# When a>64/27, there is 2-1 = 1 real root, which is negative
# When 0<a<64/27, there are 3 - 0 = 3 real roots, of which 2 are positive
# When a<0, there is 2 - 1 = 1 real root, which is positive
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;seq(coeftayl(1/subs(a=1,H),[x,y,z,w]=[0,0,0,0],[k,k,k,k]),k=0..6);###########################
# Specialize to the a=64/27 case
###########################H := subs(a=64/27,H);DOP := w^2*(256*w-81)*(16777216*w^2-5971968*w+216513)*(8388608*w^2+12773376*w-2381643)*Dw^4+3*w*(126100789566373888*w^5+145702882866364416*w^4-144989442627600384*w^3+37699919662546944*w^2-3508882176945024*w+69613650565965)*Dw^3+(12*(83879543059775488*w^5+127612480615612416*w^4-91463559904493568*w^3+18580583017193472*w^2-1333811087342316*w+13922730113193))*Dw^2+(72*(9147936743096320*w^4+18313981068312576*w^3-9089964406996992*w^2+1351389647038464*w-58491439678683))*Dw+50665495807918080*w^3+141873283866820608*w^2-42412899154526208*w+4304856990449664:
D64 := w^2*(256*w-81)*(16777216*w^2-5971968*w+216513)*(8388608*w^2+12773376*w-2381643)*(diff(f(w), w, w, w, w))+3*w*(126100789566373888*w^5+145702882866364416*w^4-144989442627600384*w^3+37699919662546944*w^2-3508882176945024*w+69613650565965)*(diff(f(w), w, w, w))+(12*(83879543059775488*w^5+127612480615612416*w^4-91463559904493568*w^3+18580583017193472*w^2-1333811087342316*w+13922730113193))*(diff(f(w), w, w))+(72*(9147936743096320*w^4+18313981068312576*w^3-9089964406996992*w^2+1351389647038464*w-58491439678683))*(diff(f(w), w))+(50665495807918080*w^3+141873283866820608*w^2-42412899154526208*w+4304856990449664)*f(w);with(gfun):
rec := diffeqtorec(D64,f(w),u(n));
PR := rectoproc({rec,seq(u(k)=coeftayl(1/H,[x,y,z,w]=[0,0,0,0],[k,k,k,k]),k=0..5)},u(n)):# Check initial terms match (takes 10 seconds)
seq(PR(k),k=0..7);
seq(coeftayl(1/H,[x,y,z,w]=[0,0,0,0],[k,k,k,k]),k=0..7);# Leading coefficient of DEQ
lc := coeff(DOP,Dw,4);# Leading critical points
factor(Groebner[Basis]([z-w,x-w,y-w,H,seq(diff(H,x)*x-diff(H,k)*k,k=[y,z,w])],plex(x,y,z,w)));GB := factor(Groebner[Basis]([1-ht*x*y*z*w,H,seq(diff(H,x)*x-diff(H,k)*k,k=[y,z,w])],plex(x,y,z,w,ht))):
ph := GB[1];solve(ph);
evalf([%]);# Highest hight is non-smooth lacuna point
solve(subs(ht=4096/81,GB));# Get points of second highest height
rho := 4096/297+16384*sqrt(3)*(1/2673);
convert(map(evala,[solve(subs(ht=rho,GB),explicit)]),set):
sbs := select(a->evala(subs(a,x*y*z*w*rho))=1,%);ASM := simplify(add(smoothASM(1,H,[w,x,y,z],subs(j,[w,x,y,z])),j=sbs));# Numerically we have asymptotic constant (at n=1000)
# This agrees with the differential equation analysis
evalf(PR(10000)/rho^10000*(10000)^(3/2));# Our calculations give the following constant
evalf(-3*sqrt(88+66*sqrt(3))*(sqrt(3)-4/3)/(11*Pi^(3/2)));