checks.py

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

8061 bytes, as of the version from 2026-11-09 12:11 (current). Recorded here, not run.

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

Run: sage -python checks.py
Reads generate.py only to compare its outputs. The alternate implementation
uses Python Fraction, nilpotent-parameter period products and Lagrange inversion;
it calls none of the generator's arithmetic helpers. The line counts are
independent Grassmannian intersection numbers. Localization controls are
transcribed from the Macaulay2 documentation, not recomputed by that package.
With NUMBERDB_CHECK_STORED=1, reapply the comparisons to the authenticated draft.
"""
import json
import os
import platform
import shutil
import sys
from fractions import Fraction as R
from math import prod
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.env import SAGE_VERSION
from check import exactness, measure, agrees_with, names_its_rings
from numberdb._generate import _stored_entries
from generate import ProjectiveCalabiYauGW, TARGETS, MAX_DEGREE, yukawa_series


def as_qq(text):
    value=R(str(text))
    return QQ(value.numerator)/QQ(value.denominator)


def conv(a, b, n):
    out = [R(0)] * (n+1)
    for i, ai in enumerate(a):
        for j, bj in enumerate(b):
            if i+j <= n:
                out[i+j] += ai*bj
    return out


def quotient(a, b, n):
    rem = list(a) + [R(0)] * (n+1-len(a))
    answer = [R(0)]*(n+1)
    for i in range(n+1):
        answer[i] = rem[i]/b[0]
        for j in range(1, min(len(b), n+1-i)):
            rem[i+j] -= answer[i]*b[j]
    return answer


def exp_by_powers(a, n):
    assert a[0] == 0
    out = [R(1)] + [R(0)]*n
    term = out[:]
    for k in range(1, n+1):
        term = [v/R(k) for v in conv(term,a,n)]
        out = [u+v for u,v in zip(out,term)]
    return out


def second_period(ell, n):
    # a_m(rho) = prod_l prod_{j=1}^{l*m}(j+l*rho) /
    #             prod_{j=1}^{m}(j+rho)^(r+4), modulo rho^3.
    jet = [R(1), R(0), R(0)]
    columns = [[R(1)], [R(0)], [R(0)]]
    for m in range(1,n+1):
        for l in ell:
            for j in range(l*(m-1)+1,l*m+1):
                jet = conv(jet,[R(j),R(l)],2)
        inv = [R(1,m),-R(1,m*m),R(1,m*m*m)]
        for _ in range(len(ell)+4):
            jet = conv(jet,inv,2)
        for i in range(3):
            columns[i].append(jet[i])
    f, g, b = columns
    h = quotient(g,f,n)
    j = quotient(b,f,n)
    hh = conv(h,h,n)
    a = [u-v/R(2) for u,v in zip(j,hh)]
    deriv = [R(k)*a[k] for k in range(1,n+1)]
    values = {}
    for d in range(1,n+1):
        exponential = exp_by_powers([-R(d)*v for v in h[:d]],d-1)
        coefficient = sum((deriv[i]*exponential[d-1-i] for i in range(d)),R(0))
        values[d] = R(prod(ell),d*d)*coefficient
    return values, columns


def mobius(n):
    sign=1
    p=2
    while p*p<=n:
        if n%p==0:
            n//=p
            sign=-sign
            if n%p==0:
                return 0
        p+=1
    return -sign if n>1 else sign


def line_count(ell):
    # Integration on Gr(2,n+1): [u^n v^(n-1)] (u-v)*c_top.
    # Coefficient list records powers of u; total degree is implicit.
    coefficients=[1]
    for l in ell:
        for i in range(l+1):
            nxt=[0]*(len(coefficients)+1)
            for k,v in enumerate(coefficients):
                nxt[k] += (l-i)*v
                nxt[k+1] += i*v
            coefficients=nxt
    n=len(ell)+3
    return coefficients[n-1]-coefficients[n]


def partitions(n, minimum=1):
    if n==0:
        yield ()
    for k in range(minimum,n+1):
        for rest in partitions(n-k,k):
            yield (k,)+rest


def check_all(stored_document=None):
    assert set(TARGETS)=={tuple(x+1 for x in p) for p in partitions(4)}
    # Known nonzero controls before the comparisons.
    assert conv([R(1),R(2)],[R(3),R(4)],2)==[3,10,8]
    assert quotient([R(1)],[R(1),R(-1)],4)==[1]*5
    assert exp_by_powers([R(0),R(1),R(0),R(0)],3)==[1,1,R(1,2),R(1,6)]
    assert line_count((5,))==2875
    gen=ProjectiveCalabiYauGW()
    params=list(gen.enumerate())
    generated={(p['multidegree'],int(p['d'])):gen.value(p,100) for p in params}
    values=generated
    if stored_document is not None:
        # The client reader handles both flat storage and nested API exports.
        values={}
        for identity,written in _stored_entries(stored_document).items():
            multidegree,degree=identity.rsplit(',',1)
            values[(multidegree,int(degree))]=as_qq(written)
        assert list(values)==list(generated), 'Stored enumeration or order differs'
    assert not exactness(values)
    assert not names_its_rings(Path(__file__).with_name('generate.py'))
    independent={}
    controls=json.loads(Path(__file__).with_name('localization-controls.json').read_text())['values']
    checked_controls=0
    gv_integrality=0
    template_parameters=((R(1,5),R(2,5),R(3,5),R(4,5)),
                         (R(1,2),R(1,2),R(1,4),R(3,4)),
                         (R(1,3),R(2,3),R(1,3),R(2,3)),
                         (R(1,2),R(1,2),R(1,3),R(2,3)),
                         (R(1,2),)*4)
    for ell,alphas in zip(TARGETS,template_parameters):
        key=','.join(map(str,ell));kappa=prod(ell);scale=prod(l**l for l in ell)
        alt,columns=second_period(ell,MAX_DEGREE)
        independent.update({(key,d):QQ(v.numerator)/QQ(v.denominator) for d,v in alt.items()})
        assert values[(key,1)]==line_count(ell)
        for d,v in enumerate(controls[key],1):
            assert values[(key,d)]==as_qq(v),(key,d,'localization')
            checked_controls+=1
        low=yukawa_series(ell,MAX_DEGREE+4)
        high=yukawa_series(ell,MAX_DEGREE+8)
        assert low[:MAX_DEGREE+1]==high[:MAX_DEGREE+1]
        assert low[0]==kappa
        # Compare factorial periods to the exact T679/T680 coefficient convention.
        coeff=R(1);harm=R(0)
        for m in range(MAX_DEGREE+1):
            assert columns[0][m]==coeff*scale**m
            assert columns[1][m]==coeff*scale**m*harm
            if m<MAX_DEGREE:
                coeff*=prod(R(m)+a for a in alphas)/R(m+1)**4
                harm+=sum((1/(R(m)+a) for a in alphas),R(0))-R(4,m+1)
        # Lambert extraction, integrality, and reconstruct all stored coefficients.
        instantons={}
        for d in range(1,MAX_DEGREE+1):
            coefficient=QQ(d)**3*values[(key,d)]
            residual=coefficient-sum((QQ(k)**3*instantons[k] for k in range(1,d) if d%k==0),QQ(0))
            instantons[d]=residual/QQ(d)**3
            assert instantons[d].denominator()==1
            inverted=sum((QQ(mobius(k))*values[(key,d//k)]/QQ(k)**3 for k in range(1,d+1) if d%k==0),QQ(0))
            assert instantons[d]==inverted
            assert sum((QQ(k)**3*instantons[k] for k in range(1,d+1) if d%k==0),QQ(0))==low[d]
            gv_integrality+=1
    assert not agrees_with(values,lambda key:independent[key])
    result={'entries':len(values),'series_comparisons':len(independent),
            'localization_documentation_comparisons':checked_controls,
            'grassmannian_line_comparisons':len(TARGETS),
            'truncations':[MAX_DEGREE+4,MAX_DEGREE+8],
            'integral_lambert_coefficients':gv_integrality,
            'mobius_comparisons':gv_integrality,
            'hypergeometric_coefficients_checked':len(TARGETS)*(MAX_DEGREE+1)*2,
            'size':measure(values),'stored':stored_document is not None,
            'python':platform.python_version(),'sage':SAGE_VERSION,
            'numberdb':numberdb.__version__ if hasattr(numberdb,'__version__') else 'see client version',
            'Macaulay2_available':shutil.which('M2') is not None}
    return result,values


if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
    doc=numberdb.table('T688') if os.environ.get('NUMBERDB_CHECK_STORED')=='1' else None
    result,values=check_all(doc)
    print(json.dumps(result,indent=2,default=str))
    print('VALUES_JSON '+json.dumps([{'params':{'multidegree':k[0],'d':str(k[1])},'number':str(v)} for k,v in values.items()]))