independent_check.py

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

6625 bytes, as of the version from 2026-11-10 22:53 (current). Recorded here, not run.

"""Independent formula checks for T710, including stored-polynomial mode.

sage -python independent_check.py --polynomials polynomials.json --output checks.json
Input is a JSON object mapping denominator strings to polynomial strings.
No import of generate.py is made. The checker uses rational polynomials modulo
Phi_L for the finite Estermann double sum, not cotangent rational functions
or NumberField.minpoly. Products over all automorphisms give a characteristic
polynomial; its squarefree part gives the minimal polynomial. Sage polynomial
arithmetic and factorization are shared infrastructure, not a separate CAS.
Direct Arb cotangent balls independently check the positive-angle real branch.
Known cases k=2,3,4,6, wrong constants, wrong scale/sign at k=4, and an
integral Python float are controls. The Galois sign action and Vasyunin inverse
conversion are checked for every reduced numerator in the stored range.
"""
import argparse
import json
import sys
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.real_arb import RealBallField

QT = PolynomialRing(QQ, 't')
QX = PolynomialRing(QQ, 'x')
ZX = PolynomialRing(ZZ, 'x')
R = RealBallField(512)


def estermann(k, h=1):
    L = int(ZZ(k).lcm(4))
    phi = QT.cyclotomic_polynomial(L)
    # Accumulate by exponent before polynomial reduction: exact finite D(0,h/k).
    c = [QQ(0)]*L
    for r in range(1, k+1):
        for s in range(1, k+1):
            c[(h*r*s*(L//k)) % L] += (QQ(1)/2-QQ(r)/k)*(QQ(1)/2-QQ(s)/k)
    D = QT(c)
    beta = (-2*QT.gen()**(L//4)*(D-QQ(1)/4)) % phi
    return L, phi, beta


def eval_mod(p, beta, phi):
    v = QT.zero()
    for c in reversed(p.list()):
        v = (v*beta+c) % phi
    return v


def primitive(q):
    p = ZX(q*q.denominator())
    p = ZX(p/p.content())
    return p if p.leading_coefficient()>0 else -p


def orbit_polynomial(L, phi, beta):
    # Multiply (x-beta(t^a)) in (Q[t]/Phi_L)[x], including repeated conjugates.
    coeff = [QT.one()]
    for a in range(1,L):
        if ZZ(a).gcd(L) != 1:
            continue
        b = beta(QT.gen()**a) % phi
        new = [QT.zero()]*(len(coeff)+1)
        for j,c in enumerate(coeff):
            new[j] = (new[j]-b*c) % phi
            new[j+1] = (new[j+1]+c) % phi
        coeff = new
    assert all(c.degree() <= 0 for c in coeff)
    characteristic = QX([c[0] for c in coeff])
    minimal = characteristic // characteristic.gcd(characteristic.derivative())
    return primitive(minimal), characteristic


def direct(k,h=1):
    return -sum(R(QQ(m)/k)*(R.pi()*QQ(m*h)/k).cot() for m in range(1,k))


def nonvacuous(v):
    return v.is_finite() and v.rad_as_ball() < R(10)**(-140)


def row(k,p):
    assert p.content()==1 and p.leading_coefficient()>0 and p.is_irreducible()
    assert all(type(c).__name__=='Integer' for c in p.list())
    L,phi,beta=estermann(k)
    assert eval_mod(p,beta,phi)==0
    independent,characteristic=orbit_polynomial(L,phi,beta)
    assert independent==p
    multiplicity=phi.degree()//p.degree()
    assert characteristic==QX(p).monic()**multiplicity
    v=direct(k)
    assert nonvacuous(v) and v>0
    residual=p(v)
    assert residual.is_finite() and residual.contains_zero()
    assert residual.rad_as_ball()<R(10)**(-90)
    # Enclosing endpoints have opposite signs; derivative interval excludes zero.
    eps=R(2)**(-400)
    left,right=v-eps,v+eps
    assert nonvacuous(left) and nonvacuous(right)
    assert p(left)*p(right)<0
    bracket=v.add_error(eps)
    assert not p.derivative()(bracket).contains_zero()
    assert not (p+1)(v).contains_zero()
    assert eval_mod(p+1,beta,phi)!=0
    count=0
    for h in range(1,k):
        if ZZ(h).gcd(k)!=1:
            continue
        a=next(a for a in range(h,L,k) if ZZ(a).gcd(L)==1)
        signed=(-1)**((a-1)//2)
        _,_,bh=estermann(k,h)
        assert beta(QT.gen()**a)%phi == signed*bh
        transformed=ZX([c*(-1)**(j+p.degree()) for j,c in enumerate(p.list())]) if signed<0 else p
        assert eval_mod(transformed,bh,phi)==0
        # Direct modular-inverse Vasyunin conversion, independently reindexed.
        inv=int(ZZ(h).inverse_mod(k))
        V=sum(R(QQ((m*h)%k)/k)*(R.pi()*QQ(m)/k).cot() for m in range(1,k))
        cv=direct(k,inv)
        assert nonvacuous(V) and nonvacuous(cv) and (V+cv).contains_zero()
        count+=1
    return {'k':k,'degree':int(p.degree()),'L':L,'ambient_degree':int(phi.degree()),
            'characteristic_multiplicity':int(multiplicity),'chars':len(str(p)),
            'value_ball':str(v),'radius':str(v.rad()),'accuracy_bits':int(v.accuracy()),
            'polynomial_residual':str(residual),'checked_residues':count,
            'root_bracket_halfwidth':'2^-400 plus direct ball radius',
            'exact_annihilation':True,'orbit_polynomial_equal':True,'irreducible':True}


def run(values):
    known={2:'x',3:'27*x^2-1',4:'2*x-1',6:'27*x^2-49'}
    for k,s in known.items():
        L,phi,beta=estermann(k)
        expected=ZX(s)
        assert eval_mod(expected,beta,phi)==0
        assert orbit_polynomial(L,phi,beta)[0]==expected
    _,phi,beta=estermann(4)
    for wrong in ['2*x+1','x-1','4*x-1']:
        assert eval_mod(ZX(wrong),beta,phi)!=0
    assert not nonvacuous(R('nan'))
    assert not nonvacuous(R(0).add_error(R(1)))
    assert set(values)=={str(k) for k in range(3,25)}
    rows=[]
    for key,text in values.items():
        result=row(int(key),ZX(text)); rows.append(result)
        print('CHECKED',key,'degree',result['degree'],'chars',result['chars'],flush=True)
    return {'table':'T710','checked_entries':len(rows),'total_entries':len(values),
            'failures':[],'known_answer_control':True,'negative_control_rejected':True,
            'working_bits':512,'rows':rows,'known_polynomials':known,
            'independence':'Estermann double sum and full orbit product; no generator import or minpoly call. Shared Sage polynomial arithmetic and irreducibility backend.',
            'skill_manifest':'../20261110T221640Z-build-skills/manifest.json',
            'skill_bundle_sha256':'d6661c2f21613524889fa3a57d3c95f8330344ca62e5a856211c08a282b7dbab'}

if __name__=='__main__':
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--polynomials',required=True)
    parser.add_argument('--output',required=True)
    args=parser.parse_args()
    with open(args.polynomials) as f: values=json.load(f)
    result=run(values)
    with open(args.output,'w') as f:json.dump(result,f,indent=2)