back to table · edit · history · where entries came from · files · download
6849 bytes, as of the version from 2026-11-02 11:21 (current). Recorded here, not run.
"""Independent checks for numberdb.org/T675.
sage -python checks.py
sage -python checks.py --stored # authenticated for the private draft
Checks all exact source factors independently with PARI at 1/z, then checks
the Hasse truncation, exact Hensel certificates and the independent power-sum
ratio at full precision. No numerical computation or root-lifting code is
shared by the power-sum method. Its known-root control runs first.
Requires generate.py, frobenius-input.json and check.py beside this file.
"""
import json
import os
import re
import sys
from collections import Counter
from math import gcd
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.libs.pari import pari
from sage.env import SAGE_VERSION
from generate import (CalabiYauUnitRoots, factors, exponent, hensel, evaluate,
R, TEMPLATES, FIBRES)
from check import exactness, measure, names_its_rings, agrees_with
# c holds coefficients of P(T)=1+c1*T+...+cd*T^d in ascending order.
def power_sum_ratio(c,p,n):
degree, modulus = len(c)-1,p**n
s = [degree]
for k in range(1,n+2):
if k <= degree:
v = -k*c[k]-sum(c[j]*s[k-j] for j in range(1,k))
else:
v = -sum(c[j]*s[k-j] for j in range(1,degree+1))
s.append(v%modulus)
assert s[n]%p
return (s[n+1]*pow(s[n],-1,modulus))%modulus
def hasse(template,z,p):
# The p-1 truncation uses rational parameters reduced modulo p;
# its factorial denominators are units throughout this range.
alpha = [QQ(a)/d for d in map(int,template.split('.'))
for a in range(1,d) if gcd(a,d)==1]
residue = lambda x: int(x.numerator())*pow(int(x.denominator()),-1,p)%p
term,total = 1,1
zmod = residue(QQ(z))
for k in range(1,p):
term = term*zmod*pow(pow(k,4,p),-1,p)%p
for a in alpha: term = term*residue(a+k-1)%p
total=(total+term)%p
return total
def parse_stored(document):
out={}
numbers=document['Numbers']
if isinstance(numbers,list):
rows=[(r['params'],r['number']) for r in numbers]
else:
rows=[]
for template,zs in numbers.items():
for z,ps in zs.items():
for p,entry in ps.items():
rows.append(({'template':template,'z':z,'p':p},
entry['number'] if isinstance(entry,dict) else entry))
for params,text in rows:
key=(str(params['template']),str(params['z']),int(params['p']))
m=re.fullmatch(r'\s*(\d+)\s*\+\s*O\((\d+)\^(\d+)\)\s*',text)
assert m, (key,text)
residue,p,n=map(int,m.groups()); assert p==key[2]
out[key]=(residue,n)
return out
def main():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
print('VERSIONS',SAGE_VERSION,str(pari('version()')),numberdb.__version__)
# Exact unit root 2 of (X-2)(X-5)(X-25)(X-125).
X=R.gen(); control=(X-2)*(X-5)*(X-25)*(X-125)
ascending=tuple(int(x) for x in control.list())
assert hensel(ascending,5,40,2)==2
assert power_sum_ratio(tuple(reversed(ascending)),5,40)==2
assert hasse('2.2.2.2','-1',5)==1 # independent manual factor has trace -4
print('PASS known-root 2 Hensel/power-sum control and nonzero Hasse control')
data=factors(); polynomials={}
H={key:pari('hgminit([%s])'%key.replace('.',',')) for key in TEMPLATES}
assert tuple(int(x) for x in pari.hgmeulerfactor(H['5'],QQ(1)/2,3).list())==(1,5,45,135,729)
# Source control is classical z=2, which has PARI t=1/2.
per_template=Counter(); excluded=[]
for (key,z,p),c in data.items():
q=pari.hgmeulerfactor(H[key],QQ(1)/QQ(z),p)
independent=tuple(int(x) for x in q.list())
assert independent==c, (key,z,p,c,independent)
a,b=-c[1],c[2]
assert c==(1,-a,b,-p**3*a,p**6) and b%p==0
assert hasse(key,z,p)==a%p, ('Hasse',key,z,p,a,hasse(key,z,p))
if a%p: per_template[key]+=1
else: excluded.append((key,z,p))
polynomials[(key,z,p)]=R(c)
# A control that would fail if the reciprocal coordinate were forgotten.
assert data[('5','1/2',7)]==(1,-35,805,-12005,117649)
assert data[('5','2',7)]==(1,10,420,3430,117649)
print('PASS all 390 PARI factor comparisons, Hasse congruences and divisibility checks')
print('COUNTS',json.dumps(dict(per_template)), 'excluded',len(excluded))
gen=CalabiYauUnitRoots(); values={}; residues={}; all_n=[]
for params in gen.enumerate():
k=(params['template'],params['z'],int(params['p']))
values[k]=gen.value(params,100)
u=values[k]['number']; n=int(u.precision_absolute())
residues[k]=(int(u.lift()),n); all_n.append(n)
assert not exactness(values) and not names_its_rings(__file__.replace('checks.py','generate.py'))
if '--stored' in sys.argv:
report=gen.verify(sample=None); print(report); assert report.ok
residues=parse_stored(numberdb.table('T675'))
assert set(residues)==set(values)
live_parent=numberdb.table('T674')['Numbers']
for (key,z,p),c in data.items():
entry=live_parent[key][z][str(p)]
assert R(entry['number'] if isinstance(entry,dict) else entry)==R(c)
print('PASS live parent factors still agree with the attached snapshot')
q_values={}
for k,(r,n) in residues.items():
p=k[2]; c=data[k]; a=-c[1]; modulus=p**n
assert n==exponent(p,100) and p**(n-1)<10**100<=modulus
assert r%p==a%p and r%p
Q=tuple(reversed(c)); derivative=[j*v for j,v in enumerate(Q)][1:]
assert evaluate(Q,r)%modulus==0
assert evaluate(derivative,r)%p==pow(a,3,p)!=0
inv=pow(r,-1,modulus)
assert evaluate(c,inv)%modulus==0
assert power_sum_ratio(c,p,n)==r%modulus, ('power sum',k)
# Quantitative low-precision check by enumerating all lifts, not Newton.
current=[a%p]; m=p
for _ in range(3):
current=[v+j*m for v in current for j in range(p)
if sum(Q[i]*(v+j*m)**i for i in range(5))%(m*p)==0]
m*=p
assert len(current)==1 and current[0]==r%m
q_values[k]=r%modulus
assert not agrees_with(q_values,lambda k: power_sum_ratio(data[k],k[2],exponent(k[2],100)))
print('PASS',len(residues),'full-precision roots, inverse zeros, derivative certificates, power-sum ratios and exhaustive four-digit lifts')
print('PRECISION',min(all_n),max(all_n))
print('MEASURE',json.dumps(measure(values)))
print('DATA_JSON',json.dumps([{'params':dict(zip(('template','z','p'),[k[0],k[1],str(k[2])])),
'residue':str(r),'precision':n} for k,(r,n) in residues.items()]))
if __name__=='__main__': main()