back to table · edit · history · where entries came from · files · download
9600 bytes, as of the version from 2026-11-01 18:35 (current). Recorded here, not run.
"""Independent arithmetic and embedding checks for numberdb.org/T657.
sage -python checks.py
sage -python checks.py --stored
For a private draft, --stored requires NUMBERDB_API_KEY. The first command
checks generated values without writing. The second reads all stored values
and applies the same independent tests. Both print a JSON evidence record.
Numerical checking uses real Euler products, not Arb's modular_eta. With
q=exp(-2*pi*sqrt(r)), truncate each Euler product after N=160 factors. Its tail
is in [1-q^(N+1)/(1-q),1]. Widening 1 by that bound is conservative and all
arithmetic is outward-rounded real ball arithmetic at 768 bits.
For B, independently substitute m=C/(x-b) into the canonical polynomial,
clear its denominator and eliminate j using PARI polresultant. For A,
eliminate B from its already checked minimal polynomial and the quadratic
A=B+C/(B-b). Require exact divisibility, irreducibility, vanishing on the
finite-product enclosure, and nonvanishing of the quotient after removing
all copies of the selected factor. Thus no other irreducible factor can
contain the specified value. Controls with known and deliberately wrong
answers precede these tests. Input H polynomials are checked with PARI polclass
and FLINT's certified Hilbert routine; canonical equations with PARI's data.
"""
import json
import os
import sys
from fractions import Fraction
from pathlib import Path
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
from sage.libs.pari import pari
from sage.schemes.elliptic_curves.cm import hilbert_class_polynomial
from sage.env import SAGE_VERSION
import generate as G
BITS = 768
TERMS = 160
RB = RealBallField(BITS)
ZX = PolynomialRing(ZZ, 'x')
x = ZX.gen()
INPUTS = json.loads((Path(__file__).resolve().parent / 'cm-inputs.json').read_text())
def precise(v, minimum=600):
assert v.is_finite(), 'nonfinite independent enclosure'
assert v.accuracy() >= minimum, 'independent enclosure too wide'
return v
def euler(q):
partial = RB(1)
power = q
for n in range(1, TERMS+1):
partial *= 1-power
power *= q
omitted_bound = q**(TERMS+1)/(1-q)
return partial * RB(1).add_error(omitted_bound)
def independent_values(p, r):
t = 2*RB.pi()*RB(r).sqrt()
q = (-t).exp()
ratio = t.exp()*(euler(q)/euler(q**p))**(24//(p-1))
b, c = 24//(p-1), ZZ(p)**(12//(p-1))
return precise(ratio+b+c/ratio), precise(ratio+b)
def certify_factor(annihilator, selected, value):
precise(value)
ann, selected = ZX(annihilator).change_ring(QQ), ZX(selected).change_ring(QQ)
assert selected.is_monic()
assert bool(pari(selected).polisirreducible())
quotient, remainder = ann.quo_rem(selected)
assert remainder == 0, 'candidate does not divide exact annihilator'
multiplicity = 1
while quotient.degree() >= selected.degree():
nxt, rem = quotient.quo_rem(selected)
if rem != 0:
break
quotient = nxt
multiplicity += 1
residual, exclusion = selected(value), quotient(value)
assert residual.is_finite() and exclusion.is_finite()
assert residual.contains_zero(), 'candidate misses the CM embedding'
assert not exclusion.contains_zero(), 'another factor may contain the CM embedding'
return multiplicity
def direct_b_annihilator(p, r):
# PARI's polynomial variables have priorities; explicitly name elimination.
b, c = 24//(p-1), p**(12//(p-1))
h = INPUTS['hilbert'][str(r)]
phi = INPUTS['canonical'][str(p)]
expr = ('polresultant(subst(%s,x,j),'
'subst(subst(%s,y,j),x,%d/(x-%d))*(x-%d)^%d,j)'
% (h, phi, c, b, b, p+1))
answer = pari(expr)
return ZX(str(answer))
def direct_a_annihilator(pb, p):
b, c = 24//(p-1), p**(12//(p-1))
expr = ('polresultant(subst(%s,x,t),(t-%d)^2-(x-%d)*(t-%d)+%d,t)'
% (str(pb).replace('**','^'), b, b, b, c))
return ZX(str(pari(expr)))
def exact_decimal(text):
f = Fraction(str(text))
return QQ(f.numerator)/QQ(f.denominator)
def within_companion(value, text):
text = str(text)
if '.' not in text and 'e' not in text.lower():
assert ZZ(text) in value
return
centre = exact_decimal(text)
mantissa, _, exponent = text.lower().partition('e')
places = len(mantissa.split('.')[1]) if '.' in mantissa else 0
unit = QQ(10)**(int(exponent or '0')-places)
# Positive differences certify containment of the full independent ball.
assert value-(centre-unit) > 0 and (centre+unit)-value > 0
def controls():
assert euler(RB(0)) == 1
a, b = independent_values(2, 1)
assert 544 in a and 536 in b
assert 536 not in a and 544 not in b
certify_factor((x-3)*(x+2), x-3, RB(3))
rejected = 0
for ann, poly, ball in [((x-3)*(x+2),x+2,RB(3)),
((x-3)*(x+2),x-4,RB(3)),
((x-3)*(x+2),x-3,RB(3).add_error(10)),
((x-3)*(x+2),x-3,RB('nan'))]:
try:
certify_factor(ann,poly,ball)
except AssertionError:
rejected += 1
assert rejected == 4
return {'known_square_lattice':True,'Euler_empty_product':True,
'known_linear_factor':True,'deliberately_wrong_or_wide_rejected':rejected}
def read_values(stored):
if not stored:
generator = G.MonstrousCMPolynomials()
return {(d['g'],int(d['r'])):generator.value(d,100)['number'] for d in generator.enumerate()}
doc = numberdb.table('T657')
numbers = doc['Numbers']
if isinstance(numbers,list):
values = {(str(row['params']['g']),int(row['params']['r'])):ZX(row['number']) for row in numbers}
else:
values = {}
for g, rows in numbers.items():
for r, entry in rows.items():
if r in ('comment','param-latex'):
continue
text = entry['number'] if isinstance(entry,dict) else entry
values[(g,int(r))] = ZX(text)
expected = [(g,r) for g in G.CLASSES for r in range(1,11)]
assert list(values) == expected, 'stored enumeration differs'
return values
def main(stored=False):
result = {'table':'T657','stored':stored,'controls':controls(),
'skill_manifest':'agents/runs/20261101T181018Z-build-skills/manifest.json',
'skill_sha256':'a2b28d9953264d7b0ca62663785c2952190c5be6d6d3597fd695d9a54d606278',
'sage':SAGE_VERSION,'pari':str(pari('version()')),
'numberdb':numberdb.__version__,'real_product_bits':BITS,'product_terms':TERMS}
values = read_values(stored)
assert len(values) == 100
# Controls precede every full-range comparison.
assert ZX(str(pari(-4).polclass())) == x-1728
assert ZX(hilbert_class_polynomial(-4,algorithm='arb')) == x-1728
for r in range(1,11):
wanted = ZX(INPUTS['hilbert'][str(r)])
assert ZX(str(pari(-4*r).polclass())) == wanted
assert ZX(hilbert_class_polynomial(-4*r,algorithm='arb')) == wanted
for p in G.PRIMES:
canonical = pari('ellmodulareqn(%d)' % p)
assert canonical[1] == 0
assert G.QMJ(str(canonical[0])) == G.QMJ(INPUTS['canonical'][str(p)])
result['hilbert_PARI_and_arb_matches'] = 10
result['canonical_PARI_matches'] = 5
details = []
for p in G.PRIMES:
for r in range(1,11):
ba, bb = independent_values(p,r)
pa, pb = values[(str(p)+'A',r)], values[(str(p)+'B',r)]
assert pa.base_ring() is ZZ and pb.base_ring() is ZZ
mb = certify_factor(direct_b_annihilator(p,r),pb,bb)
ma = certify_factor(direct_a_annihilator(pb,p),pa,ba)
# Exact A/B relation in the quotient Q[x]/(P_B).
b, c = 24//(p-1), ZZ(p)**(12//(p-1))
d = pa.degree()
composed = sum((pa[i]*((x-b)*x+c)**i*(x-b)**(d-i) for i in range(d+1)),ZX.zero())
assert composed.quo_rem(pb)[1] == 0
for suffix,poly,ball,mult in [('A',pa,ba,ma),('B',pb,bb,mb)]:
g = str(p)+suffix
within_companion(ball,INPUTS['cm_values'][g][str(r)])
production = G.quotient_ball(p,r)
production = production+b+c/production if suffix=='A' else production+b
assert production.overlaps(ball)
details.append({'g':g,'r':r,'degree':int(poly.degree()),
'characters':len(str(poly)),'factor_multiplicity':mult,
'product_accuracy_bits':int(ball.accuracy()),
'generator_accuracy_bits':int(production.accuracy()),
'embedding_ball':str(ball),
'polynomial':str(poly)})
assert values[('2A',1)] == x-544 and values[('2B',1)] == x-536
result.update(entries=100,exact_AB_relations=50,certified_factor_checks=100,
companion_enclosures=100,minimum_product_accuracy_bits=min(d['product_accuracy_bits'] for d in details),
minimum_generator_accuracy_bits=min(d['generator_accuracy_bits'] for d in details),
longest=max(details,key=lambda d:d['characters']),details=details)
print('EVIDENCE '+json.dumps(result,ensure_ascii=False),flush=True)
return result
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN'):
numberdb.configure(api_key=sys.stdin.read().strip())
main(stored='--stored' in sys.argv)