back to table · edit · history · where entries came from · files · download
4526 bytes, as of the version from 2026-11-02 01:06 (current). Recorded here, not run.
"""Independent Domb-moment continuation checks for numberdb.org/T664.
Run: sage -python independent_checks.py > independent-evidence.log
Uses no code from generate.py. The moment integral's Taylor continuation is
started in its absolutely convergent exterior, routed through the upper half
plane, and landed on each ordinary real logarithmic-potential boundary value.
Rogers (0704.2438), (2.11), and BSWZ (1103.2995), (37), give the starting
series and differential equation. No eta products or nome solving are used.
At the regular-singular arguments k=0,12 the independent hypergeometric
endpoint integral is used; k=-3 additionally uses the Dirichlet evaluation.
k=-4 uses the exact polynomial factorization and Smyth's theorem.
"""
import numberdb.sage as numberdb
import mpmath as m
import json,math,time
m.mp.dps=145
ORDER=400
# Domb coefficients, computed as exact integers by their recurrence.
D=[1,4]
for n in range(2,650):
numerator=2*(2*n-1)*(5*n*n-5*n+2)*D[-1]-64*(n-1)**3*D[-2]
assert numerator%n**3==0
D.append(numerator//n**3)
# Independent combinatorial check of recurrence coefficients.
for n in range(35):
direct=sum(math.comb(n,j)**2*math.comb(2*j,j)*math.comb(2*n-2*j,n-j) for j in range(n+1))
assert D[n]==direct
assert D[:5]==[1,4,28,256,2716]
def exterior(u):
u=m.mpc(u)
f=m.log(u)-m.fsum(D[n]/(n*u**n) for n in range(1,len(D)))
return [f,m.fsum(D[n]/u**(n+1) for n in range(len(D))),
-m.fsum((n+1)*D[n]/u**(n+2) for n in range(len(D)))/2,
m.fsum((n+1)*(n+2)*D[n]/u**(n+3) for n in range(len(D)))/6]
# sum_{d=1}^4 A[d](u) d^d g(u)/du^d = 0.
A={1:[-4,1],2:[64,-68,7],3:[0,192,-90,6],4:[0,0,64,-20,1]}
def step(a,c,h):
pol={d:[sum(p[t]*math.comb(t,j)*a**(t-j) for t in range(j,len(p))) for j in range(len(p))] for d,p in A.items()}
cc=list(c)
for n in range(ORDER-3):
val=m.mpc(0)
for d in range(1,5):
for j,p in enumerate(pol[d]):
if j>n or (d==4 and j==0):continue
t=n-j+d
val+=p*math.prod(range(t-d+1,t+1))*cc[t]
cc.append(-val/(pol[4][0]*math.prod(range(n+1,n+5))))
result=[]
for d in range(4):
total=m.mpc(0)
for j in range(len(cc)-1,d-1,-1):total=total*h+cc[j]*math.comb(j,d)
result.append(total)
return result
def travel(a,c,b):
b=m.mpc(b)
while abs(b-a)>m.mpf('1e-140'):
distance=min(abs(a),abs(a-4),abs(a-16))
h=b-a
if abs(h)>m.mpf('0.35')*distance:h*=m.mpf('0.35')*distance/abs(h)
c=step(a,c,h);a+=h
return a,c
start=time.time()
a=m.mpc(32);initial=exterior(a)
b,c=travel(a,initial,64)
err=abs(c[0]-exterior(64)[0]);assert err<m.mpf('1e-130')
print('CONTROL '+json.dumps({'exterior_continuation_32_to_64_error':m.nstr(err,12),'moments_checked':35}),flush=True)
# Cache a point safely above all three regular singularities.
a,base=travel(m.mpc(32),initial,32+16j)
a,base=travel(a,base,16+16j)
values={}
for k in sorted(set(range(-8,17))|{-36},key=lambda x:(abs(x),x<0)):
if k in [-4,0,12]:continue
u=k+4
b,c=travel(a,base,u+16j)
b,c=travel(b,c,u)
values[str(k)]=m.nstr(c[0].real,140)
print('INDEPENDENT '+json.dumps({'k':k,'value':values[str(k)],'method':'Domb ODE continuation'}),flush=True)
# A positive real 2F1 integral avoids analytic continuation and nome choices
# for the endpoint. Clausen: F=3F2(1/3,1/2,2/3;1,1;z)=2F1(1/6,1/3;1;z)^2.
# Integral (F-1)/z = 5F4(4/3,3/2,5/3,1,1;2,2,2,2;1)/9.
a1=m.mpf(1)/6;b1=m.mpf(1)/3
def integrand(z):
if not z:return m.mpf(1)/9
return m.expm1(2*m.log(m.hyp2f1(a1,b1,1,z)))/z
integral=m.quad(integrand,[0,m.mpf(1)/2,1])
fcc=16*m.log(2)/15+8*m.log(3)/5-8*integral/15
for k,value in [(12,fcc),(0,fcc/4),(-4,7*m.zeta(3)/m.pi**2)]:
values[str(k)]=m.nstr(value,140)
print('INDEPENDENT '+json.dumps({'k':k,'value':values[str(k)],'method':'hypergeometric endpoint integral' if k!=-4 else 'factorization and Smyth theorem'}),flush=True)
L=(m.zeta(2,m.mpf(1)/3)-m.zeta(2,m.mpf(2)/3))/9
special=6*m.sqrt(3)*L/(5*m.pi)
err=abs(special-m.mpf(values['-3']));assert err<m.mpf('1e-125')
# Endpoint derivative: from the exterior moment period and Gauss's 2F1 value.
F1=(m.gamma(m.mpf(1)/2)/(m.gamma(m.mpf(5)/6)*m.gamma(m.mpf(2)/3)))**2
watson=9*m.gamma(m.mpf(1)/3)**6/(m.power(2,m.mpf(14)/3)*m.pi**4)
assert abs(F1-watson)<m.mpf('1e-135')
print('SPECIAL '+json.dumps({'dirichlet_minus3_error':m.nstr(err,12),'watson_error':m.nstr(abs(F1-watson),12),'fcc':m.nstr(fcc,140)}),flush=True)
# A control that a deliberately wrong value fails the desired comparison.
assert abs(special-(m.mpf(values['-3'])+m.mpf('1e-90')))>m.mpf('1e-110')
print('SECONDS',time.time()-start,flush=True)