back to table · edit · history · where entries came from · files · download
9877 bytes, as of the version from 2026-10-31 22:58 (current). Recorded here, not run.
"""Independent integral checks for Selberg errors -- numberdb.org/T642.
sage -python checks.py > numerical-evidence.json
sage -python checks.py --stored > stored-evidence.json
The second command reads T642 and the companion tables; it never writes.
Download generate.py beside this file. NUMBERDB_API_KEY must see the draft.
The campaign runs all Sage commands via agents/sage.sh.
No cotangent evaluation or generator kernel is used in the independent path.
For gcd(j,k)=d, j=dp, k=dq, integrate {pu}{qu} psi_1(u)/(dpq)
over [0,1]. Each piece is quadratic. Terms n=0..15 of psi_1 are
integrated exactly in elementary functions; the tail is expanded through
u^111 and bounded explicitly. See derivation.md for the remainder proof.
"""
import sys,os,json,time,functools,hashlib
from pathlib import Path
from fractions import Fraction
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ
from sage.rings.real_arb import RealBallField
from sage.version import version as sage_version
from numberdb._write import to_text
from generate import SelbergErrors, kernel, weights, mobius_values
BITS=640
K=16
M=112
R=RealBallField(BITS)
# Each coefficient encloses (-1)^m(m+1) sum_{n=K}^infinity n^(-m-2).
@functools.lru_cache(None)
def tail_coefficients():
return tuple((-1)**m*(m+1)*(R(m+2).zeta()-
sum((R(n)**(-m-2) for n in range(1,K)),R(0)))
for m in range(M))
@functools.lru_cache(None)
def integral_product(p,q):
cuts=sorted({QQ(a)/p for a in range(p+1)}|{QQ(b)/q for b in range(q+1)})
total=R(0)
moments=[R(0) for _ in range(M)]
for left,right in zip(cuts,cuts[1:]):
middle=(left+right)/2
ip,iq=(p*middle).floor(),(q*middle).floor()
a=ZZ(p*q);b=-ZZ(p*iq+q*ip);c=ip*iq
lo,hi=R(left),R(right)
for n in range(K):
if n==0 and left==0:
assert b==0 and c==0
total+=a*(hi-lo)
else:
low,high=lo+n,hi+n
total+=a*(hi-lo)+(b-2*a*n)*(high/low).log()+\
(a*n*n-b*n+c)*(1/low-1/high)
lp,hp=lo,hi
for m in range(M):
moments[m]+=a*(hp*hi*hi-lp*lo*lo)/(m+3)+\
b*(hp*hi-lp*lo)/(m+2)+c*(hp-lp)/(m+1)
lp*=lo;hp*=hi
total+=sum((coef*moment for coef,moment in zip(tail_coefficients(),moments)),R(0))
# |P(u)|<=1, hence the integrated alternating-series remainder is
# <= sum n^(-M-2) <= K^(-M-2)+K^(-M-1)/(M+1).
bound=R(QQ(1)/ZZ(K)**(M+2)+QQ(1)/(ZZ(M+1)*ZZ(K)**(M+1)))
return total.add_error(bound)
def gram(j,k):
d=int(ZZ(j).gcd(ZZ(k)));p,q=j//d,k//d
p,q=sorted((p,q))
return integral_product(p,q)/(d*p*q)
def mu_by_divisors(n):
# Independent exact trial-division count, including square detection.
if n==1:return 1
count=0
for p in range(2,n+1):
if n%p or any(p%d==0 for d in range(2,p)):continue
if n%(p*p)==0:return 0
count+=1
return (-1)**count
def solve(G,b):
"""Independent Gaussian elimination in balls, with positive-pivot checks."""
N=len(b);H=[list(G[i][:N])+[b[i]] for i in range(N)]
pivots=[]
for k in range(N):
assert H[k][k]>0
pivots.append(H[k][k])
for i in range(k+1,N):
f=H[i][k]/H[k][k]
for j in range(k+1,N+1):H[i][j]-=f*H[k][j]
x=[R(0) for _ in range(N)]
for i in range(N-1,-1,-1):
x[i]=(H[i][N]-sum((H[i][j]*x[j] for j in range(i+1,N)),R(0)))/H[i][i]
return x,pivots
def decimal_interval(text):
text=str(text)
f=Fraction(text)
mant,_,ex=text.lower().partition('e')
scale=(int(ex) if ex else 0)-(len(mant.split('.')[1]) if '.' in mant else 0)
unit=QQ(10)**scale
center=QQ(f.numerator)/QQ(f.denominator)
return center-unit,center+unit
def strict_inside(x,text):
lo,hi=decimal_interval(text)
return x.is_finite() and x.lower()>lo and x.upper()<hi
def healthy(x,threshold=QQ(1)/10**110):
assert x.is_finite() and x.rad()<threshold
def live_numbers(tid):
doc=numberdb.table(tid)
block=doc['Numbers']
return {int(k):(v['number'] if isinstance(v,dict) else v) for k,v in block.items()}
def main():
start=time.monotonic()
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
key=sys.stdin.read().strip()
if '=' in key and key.split('=',1)[0].isupper():key=key.split('=',1)[1].strip().strip('\'"')
numberdb.configure(api_key=key)
generator=SelbergErrors()
# Known nonzero answers before comparisons.
C=(2*R.pi()).log()-R.euler_constant()
unit=gram(1,1);healthy(unit)
assert unit.overlaps(C) and not (unit+QQ(1)/100).overlaps(C)
exact2=3+(2*R.pi()).log()-3*R.euler_constant()
first=generator.value({'N':'2'},100);healthy(first)
assert first.overlaps(exact2)
wrong_sign=1-2*(1-R.euler_constant())+C
assert not wrong_sign.overlaps(first)
for n in range(1,31):assert mu_by_divisors(n)==mobius_values(30)[n]
# The target integral is independent of cotangents. Its elementary
# derivation is attached, and finite harmonic partial sums supply a control.
H=sum((QQ(1)/n for n in range(1,1001)),QQ(0))
partial=R(1001).log()-R(H)+QQ(1000)/1001
target=1-R.euler_constant()
assert partial < target and target < partial+QQ(1)/1001
G=[[gram(j,k) for k in range(1,31)] for j in range(1,31)]
prodR,prodG,prodb=kernel(30,numberdb.bits(100,losing=128))
for i in range(30):
for j in range(30):
healthy(G[i][j]);assert G[i][j].overlaps(prodG[i][j])
assert prodG[i][j].overlaps(prodG[j][i])
b=[(R(n).log()+target)/n for n in range(1,31)]
records=[];max_radius=0;min_accuracy=10000
stored_mode='--stored' in sys.argv or os.environ.get('SELBERG_STORED')=='1'
actual=live_numbers('T642') if stored_mode else None
companion=live_numbers('T641') if stored_mode else None
if actual is not None:assert list(actual)==list(range(2,31))
for N in range(2,31):
# Independent coefficients use log(N/n)/log(N), not 1-log(n)/log(N).
w=[R(1) if n==1 else R(mu_by_divisors(n))*R(QQ(N)/n).log()/R(N).log() for n in range(1,N+1)]
assert w[0]==1 and w[-1].is_zero()
known_w=weights(N,prodR)
for i in range(N):assert w[i].overlaps(known_w[i])
e=1+2*sum((w[i]*b[i] for i in range(N)),R(0))+\
sum((w[i]*G[i][j]*w[j] for i in range(N) for j in range(N)),R(0))
result=generator.value({'N':str(N)},100)
healthy(e);healthy(result)
assert e.overlaps(result)
text=to_text(result,100)
assert strict_inside(e,text) and strict_inside(result,text)
c,pivots=solve(G,b[:N])
d2=1-sum((b[i]*c[i] for i in range(N)),R(0))
healthy(d2)
assert 0<d2<e
q=[w[i]+c[i] for i in range(N)]
excess=sum((q[i]*G[i][j]*q[j] for i in range(N) for j in range(N)),R(0))
healthy(excess);assert excess>0 and excess.overlaps(e-d2)
if actual is not None:
assert strict_inside(e,actual[N]) and strict_inside(result,actual[N])
assert strict_inside(d2,companion[N])
# A bad decimal at ten last places must be rejected.
lo,hi=decimal_interval(actual[N]);ulp=(hi-lo)/2
assert not strict_inside(e+10*ulp,actual[N])
storedlo,storedhi=decimal_interval(actual[N])
dlo,dhi=decimal_interval(companion[N])
assert storedlo>dhi
rec={'N':N,'stored_decimal' if actual else 'decimal':actual[N] if actual else text,
'generator_ball':str(result),'independent_ball':str(e),
'generator_accuracy_bits':int(result.accuracy()),
'independent_radius':str(e.rad()),'minimum_ball':str(d2),
'excess_ball':str(excess),'positive_pivots':len(pivots)}
records.append(rec)
max_radius=max(max_radius,e.rad());min_accuracy=min(min_accuracy,result.accuracy())
evidence={'table':'T642','mode':'stored' if stored_mode else 'pre-fill',
'skill_manifest':'agents/runs/20261031T225043Z-build-skills/manifest.json',
'skill_sha256':{'local':'a2b28d9953264d7b0ca62663785c2952190c5be6d6d3597fd695d9a54d606278',
'live':'a2b28d9953264d7b0ca62663785c2952190c5be6d6d3597fd695d9a54d606278'},
'versions':{'sage':sage_version,'numberdb':numberdb.__version__,'python':sys.version.split()[0]},
'parameters':{'N':{'category':'discrete_index','extent':'2 <= N <= 30'}},
'generator_working_bits':int(prodR.precision()),'independent_working_bits':BITS,
'direct_trigamma_terms':K,'tail_expansion_terms':M,
'entry_count':len(records),'gram_entries_checked':900,'symmetric_gram_entries':465,
'min_generator_accuracy_bits':int(min_accuracy),'max_independent_radius':str(max_radius),
'known_answer_controls':['G(1,1)=log(2*pi)-gamma','E(2)=3+log(2*pi)-3*gamma',
'target integral between harmonic partial sum and explicit 1/1001 tail'],
'negative_controls':['G(1,1)+1/100 disjoint','wrong residual sign rejected']+
(['ten-last-place stored perturbation rejected for all 29 entries'] if actual else []),
'checks':['finite narrow balls','independent integral agreement',
'all weights including exact endpoints','Gram symmetry',
'strict containment in 100-digit last-place intervals',
'E_N > independently solved d_N^2','E_N-d_N^2=(w+c)^T G (w+c)'],
'elapsed_seconds':round(time.monotonic()-start,3),'records':records,
'source_sha256':{name:hashlib.sha256(Path(__file__).with_name(name).read_bytes()).hexdigest()
for name in ['generate.py','checks.py']}}
print(json.dumps(evidence,indent=2))
if __name__=='__main__':main()