"""Independent period-integral checks for T640. Run: sage -python checks.py For a read-back check also set NUMBERDB_READBACK=1 and NUMBERDB_API_KEY, or run checks.py --stored. Reads only; never writes to NumberDB. The campaign wrapper uses NUMBERDB_KEY_FROM_STDIN=1 for a key piped on stdin. Proof of the independent remainder bound: derivation.md. """ import os import sys import json import time from functools import lru_cache from fractions import Fraction from math import gcd import numberdb.sage as numberdb from sage.rings.rational_field import QQ from sage.rings.real_arb import RealBallField from sage.env import SAGE_VERSION from numberdb._write import to_text from numberdb._compare import _claim, digits_of from generate import FractionalPartAutocorrelation, autocorrelation, vasyunin BITS = 512 M = 16 L = 100 R = RealBallField(BITS) def qfraction(x): """Convert a Python exact Fraction to a Sage rational without division floats.""" x = Fraction(x) return QQ(x.numerator)/QQ(x.denominator) def healthy(x, decimal_radius=110): assert x.is_finite(), 'A nonfinite ball cannot be evidence' assert x.rad_as_ball() < R(10)**(-decimal_radius), 'Comparison enclosure too wide' def primitive_integral(a, b, c, left, right, n): """Integral of (a*u²+b*u+c)/(n+u)² on [left,right].""" B, D = b-2*a*n, a*n*n-b*n+c if n == 0 and left == 0: assert B == 0 and D == 0 return R(a*(right-left)) lo, hi = left+n, right+n return (R(a*(right-left)) + R(B)*R(hi/lo).log() + R(D*(1/lo-1/hi))) @lru_cache(maxsize=None) def tail_coefficients(): return tuple((-1)**r*(r+1)*R(r+2).zeta(R(M)) for r in range(L+1)) def period_integral(p, q): """No cotangents, Euler constant or rational-A formula are used here.""" breaks = sorted({QQ(j)/QQ(q) for j in range(q+1)} | {QQ(k)/QQ(p) for k in range(p+1)}) finite = R(0) moments = [QQ(0)]*(L+1) for left, right in zip(breaks,breaks[1:]): midpoint = (left+right)/2 floorq, floorp = (q*midpoint).floor(), (p*midpoint).floor() a, b, c = QQ(p*q), QQ(-q*floorp-p*floorq), QQ(floorq*floorp) for n in range(M): finite += primitive_integral(a,b,c,left,right,n) # Reuse exact rational powers; moments are integrated before rounding. lp, rp = left, right for r in range(L+1): moments[r] += (a*(rp*right*right-lp*left*left)/(r+3) +b*(rp*right-lp*left)/(r+2)+c*(rp-lp)/(r+1)) lp *= left rp *= right total = (finite + sum((R(mu)*co for mu,co in zip(moments,tail_coefficients())),R(0)))/q bound = (QQ(L+2)/QQ(M)**(L+3) + QQ(1)/QQ(M)**(L+2))/q total = total.add_error(R(bound)) healthy(total) return total def interval_from_text(text): center, radius = _claim(str(text)) center, radius = qfraction(center), qfraction(radius) return R(center).add_error(R(radius)) def contains_written(text, computed): center, radius = _claim(str(text)) center, radius = qfraction(center), qfraction(radius) # Strict comparison certifies containment, not merely overlapping midpoints. assert abs(computed-R(center)) < R(radius), 'Written interval does not contain enclosure' def independent_c0(h,k): if k == 1: return R(0) return -sum((R(QQ(m)/QQ(k)) * (R.pi()*R(QQ((m*h) % k)/QQ(k))).cot() for m in range(1,k)),R(0)) def run(stored=False): started=time.monotonic() # Exact known-answer controls come before the measurements. assert QQ(1)/2 in primitive_integral(QQ(0),QQ(0),QQ(1),QQ(0),QQ(1),1) assert QQ(1) in primitive_integral(QQ(1),QQ(0),QQ(0),QQ(0),QQ(1),0) unit=period_integral(1,1) C=(2*R.pi()).log()-R.euler_constant() healthy(C) assert unit.overlaps(C) assert not (unit+R(QQ(1)/10**90)).overlaps(C) print('Known-answer primitive and unit controls; shifted-value rejection: passed',flush=True) gen=FractionalPartAutocorrelation() params=list(gen.enumerate()) expected=sorted({Fraction(p,q) for p in range(1,31) for q in range(1,31)}, key=lambda x:(max(x.numerator,x.denominator),x.numerator,x.denominator)) assert [x['lambda'] for x in params]==[str(x) for x in expected] assert len(params)==555 document = numberdb.table('T640') if stored else None if stored: from numberdb._generate import _stored_entries rows=_stored_entries(document) assert list(rows)==[x['lambda'] for x in params], 'Stored identity order differs' assert list(document['Parameters'])==['lambda'] assert document['Display properties']['number-header']==r'$A(\lambda)$' else: rows={} out=[]; values={}; independent={}; written={} minimum=10000;maxrad=R(0);maxcheckrad=R(0) for ix,pa in enumerate(params): lam=QQ(pa['lambda']);p,q=int(lam.numerator()),int(lam.denominator()) raw=gen.value(pa,100);value=raw['number'] if isinstance(raw,dict) else raw healthy(value,125) other=period_integral(p,q) assert value.overlaps(other), 'Independent integral disagreement at '+str(lam) text=to_text(value,100) assert digits_of(text)>=100 contains_written(text,value) contains_written(text,other) if stored: storedtext=rows[str(lam)]['number'] assert digits_of(storedtext)>=100 contains_written(storedtext,other) contains_written(storedtext,value) text=storedtext values[str(lam)]=value;independent[str(lam)]=other;written[str(lam)]=text minimum=min(minimum,int(value.accuracy())) maxrad=max(maxrad,value.rad_as_ball());maxcheckrad=max(maxcheckrad,other.rad_as_ball()) out.append({'lambda':str(lam),'number':text,'computed_ball':str(value), 'independent_ball':str(other),'accuracy_bits':int(value.accuracy())}) if (ix+1)%50==0: print('Independent integrals:',ix+1,'of 555',flush=True) for lam,value in values.items(): inverse=str(1/QQ(lam)) assert value.overlaps(R(QQ(lam))*values[inverse]) assert independent[lam].overlaps(R(QQ(lam))*independent[inverse]) if stored: assert interval_from_text(written[lam]).overlaps(R(QQ(lam))*interval_from_text(written[inverse])) # Each conversion used in the finite formula, with every reduced residue. count=0 for k in range(2,31): for h in range(1,k): if gcd(h,k)!=1:continue inverse=pow(h,-1,k) v=vasyunin(h,k,BITS);other=-independent_c0(inverse,k) healthy(v);healthy(other);assert v.overlaps(other);count+=1 for j in range(1,31): for k in range(1,31): gjk=interval_from_text(written[str(QQ(j)/QQ(k))])/j gkj=interval_from_text(written[str(QQ(k)/QQ(j))])/k assert gjk.overlaps(gkj) assert (interval_from_text(written['1'])/j).overlaps(C/j) assert interval_from_text(written['1']).overlaps(C) assert not interval_from_text(written['1']).overlaps(C+R(QQ(1)/10**90)) evidence={'table':'T640','stage':'build-readback' if stored else 'build-offline', 'skill_manifest':'agents/runs/20261031T203213Z-build-skills/manifest.json', 'skill_sha256':'a2b28d9953264d7b0ca62663785c2952190c5be6d6d3597fd695d9a54d606278', 'entries':len(out),'main_bits':numberdb.bits(100,losing=128),'check_bits':BITS, 'minimum_accuracy_bits':minimum,'maximum_generator_radius':str(maxrad), 'maximum_independent_radius':str(maxcheckrad),'tail_M':M,'tail_L':L, 'all_independent_integrals_agree':True,'written_intervals_contain_both_enclosures':True, 'reciprocal_checks':555,'modular_inverse_checks':count,'gram_symmetry_checks':900, 'controls':['primitive constant 1/(1+u)^2 -> 1/2','removable endpoint u^2/u^2 -> 1', 'period integral at 1 -> log(2pi)-gamma','10^-90 perturbation rejected'], 'seconds':time.monotonic()-started,'versions':{'sage':SAGE_VERSION,'numberdb':numberdb.__version__,'python':sys.version}, 'values':out} print('EVIDENCE_JSON '+json.dumps(evidence),flush=True) return evidence if __name__=='__main__': if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1': token=sys.stdin.read().strip() if '=' in token and token.split('=',1)[0].isupper(): token=token.split('=',1)[1].strip().strip('\'"') os.environ['NUMBERDB_API_KEY']=token run('--stored' in sys.argv or os.environ.get('NUMBERDB_READBACK')=='1')