checks.py

back to table · edit · history · where entries came from · files · download

10012 bytes, as of the version from 2026-11-02 07:16 (current). Recorded here, not run.

"""Independent checks for numberdb.org/T671.

sage -python checks.py              # offline, using frozen companion documents
sage -python checks.py --stored     # repeat on the live draft (API key needed)

The checker enumerates with Sage BinaryQF, evaluates J_m as a sum of Hecke
translates, and evaluates j through modular lambda, not the generator's Faber
polynomials or modular_j call. It also reconstructs complete primitive Hilbert
class polynomials, checks untwisted traces and rational CM values against
frozen NumberDB companions, and checks printed twisted coefficients.
Only the proposed results being tested are imported from generate.py.
"""
import json
import os
import sys
from pathlib import Path
from functools import lru_cache
from math import gcd
from decimal import Decimal
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.quadratic_forms.binary_qf import BinaryQF_reduced_representatives
import generate as gen

BITS = 1024
C = ComplexBallField(BITS)
R = RealBallField(BITS)
EPS = QQ(1)/ZZ(10)**150
DATA = json.loads(Path(__file__).with_name('companions.json').read_text())
P = PolynomialRing(ZZ,'x')
x = P.gen()


def narrow(z,limit=EPS):
    z=C(z)
    assert z.real().is_finite() and z.imag().is_finite()
    assert z.real().rad()<limit and z.imag().rad()<limit
    return z


def same(z,w):
    z=narrow(z); w=narrow(w)
    assert (z-w).contains_zero(), (z,w)


def j_other(tau):
    lam=tau.modular_lambda()
    return 256*(1-lam+lam**2)**3/(lam**2*(1-lam)**2)


@lru_cache(None)
def forms(N,primitive=False):
    return tuple(tuple(int(v) for v in Q) for Q in
                 BinaryQF_reduced_representatives(-N,primitive_only=primitive))


@lru_cache(None)
def at_form(N,Q,m):
    a,b,c=Q
    tau=C(QQ(-b)/(2*a),R(N).sqrt()/(2*a))
    # Double cosets of integral determinant-m matrices, with no division by m.
    total=C(0)
    for e in range(1,m+1):
        if m%e==0:
            for r in range(e):
                total += j_other(((m//e)*tau+r)/e)-744
    return narrow(total)


def chi_other(D,Q):
    a,b,c=Q
    if gcd(D,gcd(a,gcd(b,c)))!=1:
        return ZZ(0)
    signs=set()
    for u in range(-4,5):
        for v in range(-4,5):
            n=a*u*u+b*u*v+c*v*v
            if n and gcd(n,D)==1:
                signs.add(ZZ(D).kronecker(n))
    assert len(signs)==1,(D,Q,signs)
    return signs.pop()


def weight_other(Q):
    a,b,c=Q
    if b==0 and a==c: return ZZ(2)
    if a==b and b==c: return ZZ(3)
    return ZZ(1)


@lru_cache(None)
def trace_other(m,D,d):
    total=C(0)
    for Q in forms(D*d):
        total += chi_other(D,Q)*at_form(D*d,Q,m)/weight_other(Q)
    return narrow(total)


def val(v):
    return v['number'] if isinstance(v,dict) else v


def controls():
    same(j_other(C(0,1)),1728)
    same(j_other(C(QQ(-1)/2,R(3).sqrt()/2)),0)
    assert not (narrow(j_other(C(0,1)))-1729).contains_zero()
    for z in (C(R('0 +/- 1')),C(R('nan'))):
        try: narrow(z)
        except AssertionError: pass
        else: raise AssertionError('Bad enclosure passed the control')
    assert forms(15)==((1,1,4),(2,1,2))
    assert [chi_other(5,Q) for Q in forms(15)]==[1,-1]
    assert chi_other(5,(5,5,5))==0
    # A leading coefficient divisible by D must not be mistaken for zero.
    assert chi_other(5,(5,5,1))==1
    print('Known-answer and deliberately wrong/wide/nonfinite controls passed',flush=True)


def check_companions(pair_list):
    Ns=sorted({D*d for D,d in pair_list})
    primitive_Ns=set(Ns)
    for N in Ns:
        for a,b,c in forms(N):
            f=gcd(a,gcd(b,c));primitive_Ns.add(N//(f*f))
    poly_coeffs=0
    for N in sorted(primitive_Ns):
        coeffs=[C(1)]
        for Q in forms(N,True):
            j=at_form(N,Q,1)+744
            out=[C(0)]*(len(coeffs)+1)
            for k,c in enumerate(coeffs):
                out[k]-=j*c;out[k+1]+=c
            coeffs=out
        expected=P(str(val(DATA['T129']['Numbers'][str(-N)])))
        assert len(coeffs)==expected.degree()+1
        for k,c in enumerate(coeffs):
            same(c,expected[k]);poly_coeffs+=1
    rational=0
    for disc,v in DATA['T83']['Numbers'].items():
        N=-int(disc);qs=forms(N,True);assert len(qs)==1
        same(at_form(N,qs[0],1)+744,ZZ(val(v)));rational+=1
    ordinary=0
    for N in Ns:
        if N>200:continue
        for m in range(1,5):
            s=sum((at_form(N,Q,m)/weight_other(Q) for Q in forms(N)),C(0))
            same(s,ZZ(val(DATA['T643']['Numbers'][str(m)][str(N)])));ordinary+=1
    for m in range(2,5):
        assert gen.faber(m,x)==P(DATA['T635']['Numbers']['1A'][str(m)].replace('X','x'))
    # Exact finite Laurent arithmetic verifies the principal part from stored j.
    from sage.rings.laurent_series_ring import LaurentSeriesRing
    L=LaurentSeriesRing(ZZ,'q',default_prec=12);q=L.gen()
    J=sum((ZZ(v)*q**int(k) for k,v in DATA['T79']['Numbers'].items() if -1<=int(k)<=10 and int(k)!=0),L(0)).add_bigoh(11)
    for m in range(1,5):
        S=gen.faber(m,J)
        for k in range(-m,1):assert S[k]==(1 if k==-m else 0)
    return dict(hilbert_polynomials=len(primitive_Ns),hilbert_coefficients=poly_coeffs,rational_CM=rational,ordinary_traces=ordinary,faber_polynomials=3,principal_parts=4)


def ordered_records(doc):
    block=doc['Numbers']
    if isinstance(block,list):return block
    out=[]
    for m,ds in block.items():
        for D,dds in ds.items():
            for d,norms in dds.items():
                for norm,v in norms.items():out.append({'params':dict(m=m,D=D,d=d,normalisation=norm),'number':val(v)})
    return out


def decimal_bounds(text):
    dec=Decimal(str(text))
    t=dec.as_tuple();n=ZZ(''.join(map(str,t.digits)))
    if t.sign:n=-n
    unit=QQ(10)**t.exponent
    center=n*unit
    return center-unit,center+unit


def main(stored=False):
    controls()
    # Independent enumeration uses squarefree criteria instead of the generator.
    Ds=[]
    for D in range(2,101):
        if D%4==1 and ZZ(D).is_squarefree():Ds.append(D)
        elif D%4==0 and (D//4)%4 in (2,3) and ZZ(D//4).is_squarefree():Ds.append(D)
    pair_list=[(D,d) for D in Ds for d in range(3,301//D+1) if d%4 in (0,3)]
    assert pair_list==list(gen.pairs())
    groups={'fundamental_coprime':0,'fundamental_shared':0,'order_coprime':0,'order_shared':0}
    for D,d in pair_list:
        group=('fundamental' if ZZ(-d).is_fundamental_discriminant() else 'order')+('_coprime' if gcd(D,d)==1 else '_shared')
        groups[group]+=1
        assert set(forms(D*d))==set(gen.reduced_forms(D*d))
        assert sum((QQ(chi_other(D,Q))/weight_other(Q) for Q in forms(D*d)),QQ(0))==0
        for Q in forms(D*d):
            chi=chi_other(D,Q);assert chi==gen.character(D,Q)
            a,b,c=Q
            assert chi_other(D,(a,b+2*a,a+b+c))==chi
            assert chi_other(D,(c,-b,a))==chi
    assert groups=={'fundamental_coprime':84,'fundamental_shared':32,'order_coprime':13,'order_shared':10}
    # Three explicit coefficients, including the common-prime cases of Zagier ยง7.
    for D,d,A in [(5,3,-85995),(8,4,18473000),(5,15,-292658282496)]:
        same(trace_other(1,D,d)/R(D).sqrt(),A)
    print('Enumeration, character invariance, constant cancellation, and printed values passed',flush=True)
    companion_checks=check_companions(pair_list)
    print('Companion checks',json.dumps(companion_checks),flush=True)
    triples=[];min_accuracy=10**9;min_trace_accuracy=10**9;max_rad=R(0)
    for m in range(1,5):
        for D,d in pair_list:
            independent=trace_other(m,D,d)/R(D).sqrt()
            A=gen.coefficient(m,D,d)
            same(independent,A)
            b=gen.coefficient_ball(m,D,d)
            assert A in b
            min_accuracy=min(min_accuracy,int(b.accuracy()))
            max_rad=max(max_rad,R(b.rad()))
            t=gen.TwistedTraces().value(dict(m=m,D=D,d=d,normalisation='trace'),100)['number']
            if A:min_trace_accuracy=min(min_trace_accuracy,int(t.accuracy()))
            triples.append({'m':m,'D':D,'d':d,'A':str(A)})
        print('Independent Hecke sums passed for m='+str(m),flush=True)
    expected_order=[(m,D,d,norm) for m in range(1,5) for D,d in pair_list for norm in ('trace','coefficient')]
    assert [(p['m'],p['D'],p['d'],p['normalisation']) for p in gen.TwistedTraces().enumerate()]==expected_order
    readback=0
    if stored:
        doc=numberdb.table('T671');records=ordered_records(doc)
        got=[]
        for r in records:
            p=r['params'];m,D,d=(int(p[k]) for k in ('m','D','d'));norm=p['normalisation'];got.append((m,D,d,norm))
            A=gen.coefficient(m,D,d);text=str(r['number'])
            if norm=='coefficient' or not A:assert text==str(A)
            else:
                lo,hi=decimal_bounds(text);t=trace_other(m,D,d).real()
                assert t>lo and t<hi
                converted=R(D).sqrt()*A
                assert converted>lo and converted<hi
            readback+=1
        assert got==expected_order
    report={'table':'T671','pairs':len(pair_list),'D_count':len(Ds),'groups':groups,'entries':len(expected_order),'independent_Hecke_sums':len(triples),'minimum_coefficient_accuracy_bits':min_accuracy,'maximum_coefficient_radius':str(max_rad),'minimum_output_trace_accuracy_bits':min_trace_accuracy,'stored_entries_checked':readback,**companion_checks}
    print('REPORT_JSON '+json.dumps(report),flush=True)
    print('COEFFICIENTS_JSON '+json.dumps(triples),flush=True)


if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        key=sys.stdin.read().strip()
        if '=' in key and key.split('=',1)[0].isidentifier():key=key.split('=',1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY']=key
    main('--stored' in sys.argv or os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1')