back to table · edit · history · where entries came from · files · download
7447 bytes, as of the version from 2026-11-02 16:48. Recorded here, not run.
"""Independent checks for T681. Run: sage -python checks.py.
Keep generate.py and check.py beside this file. Optional --stored reads the
private draft using NUMBERDB_API_KEY, or NUMBERDB_KEY_FROM_STDIN=1. No writes.
The all-entry Hurwitz-zeta recomputation is independent of the generator's
polylogarithm algorithm, but shares LS Theorem 6; direct quadrature in
quadrature.py checks that theorem against the defining integral separately.
"""
import sys
import os
import json
from functools import lru_cache
from math import gcd
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.real_mpfi import RealIntervalField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.arith.misc import moebius
from numberdb._write import to_text
from check import exactness, agrees_with, names_its_rings, measure
from generate import HigherCyclotomicMahler, moment, second_coefficient
PRECISION = 512
R = RealBallField(PRECISION)
def primitive_angles(n):
return [QQ(a) / n for a in range(n) if gcd(a, n) == 1]
def divisor_coefficient(n):
ds = ZZ(n).divisors()
return sum((QQ(moebius(n // a) * moebius(n // b) * gcd(a, b) ** 2)
/ (12 * a * b) for a in ds for b in ds), QQ(0))
def root_second_coefficient(n):
angles = primitive_angles(n)
return sum(((a - b)**2 - abs(a - b) + QQ(1)/6
for a in angles for b in angles), QQ(0)) / 2
@lru_cache(maxsize=None)
def hurwitz_values(b, s):
return tuple(R(s).zeta(R(QQ(r)/b)) for r in range(1, b + 1))
@lru_cache(maxsize=None)
def fourier_hurwitz(delta, s):
"""Group the absolutely convergent Fourier series by residue class."""
a, b = int(delta.numerator()), int(delta.denominator())
total = R(0)
for r, zeta in enumerate(hurwitz_values(b, s), start=1):
theta = 2 * R.pi() * (QQ((a * r) % b) / b)
trig = theta.cos() if s == 3 else theta.sin()
total += trig * zeta
return total / ZZ(b)**s
def third_from_hurwitz(n):
angles = primitive_angles(n)
d = len(angles)
csum, ssum = R(0), R(0)
# Deliberately uncompressed and unreflected; indexes are root ranks.
for j in range(d):
for i in range(j):
delta = angles[j] - angles[i]
csum += fourier_hurwitz(delta, 3)
ssum += (d * delta - (j - i)) * fourier_hurwitz(delta, 2)
return -(R(3)/2) * d*d * R(3).zeta() - 3*d*csum - 3*R.pi()*ssum
def meaningful_overlap(a, b, bits=350):
assert a.is_finite() and b.is_finite()
assert a.accuracy() >= bits and b.accuracy() >= bits
assert a.overlaps(b), (a, b)
def rational_text(text):
"""Exact rational parsing of the decimal spelling, without float coercion."""
from decimal import Decimal
d = Decimal(text)
sign, digits, exponent = d.as_tuple()
mantissa = int(''.join(map(str, digits))) * (-1 if sign else 1)
return QQ(mantissa) * QQ(10)**int(exponent)
def decimal_contains(text, ball):
from decimal import Decimal
centre = rational_text(text)
unit = QQ(10)**int(Decimal(text).as_tuple().exponent)
interval = RealIntervalField(ball.parent().precision())(ball)
return (centre-unit <= interval.lower().exact_rational()
and interval.upper().exact_rational() <= centre+unit)
def main():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
token = sys.stdin.read().strip()
if '=' in token and token.split('=',1)[0].isidentifier():
token = token.split('=',1)[1].strip().strip('\"\'')
os.environ['NUMBERDB_API_KEY'] = token
generator = HigherCyclotomicMahler()
selection = list(generator.enumerate())
independent_ns = [n for n in range(1, 106, 2) if moebius(n) != 0]
assert selection == [{'n':n,'k':k} for n in independent_ns for k in (2,3)]
assert not names_its_rings('generate.py')
# Positive and signed controls before believing any agreement measurement.
meaningful_overlap(moment(1,2), R.pi()**2/12)
meaningful_overlap(moment(1,3), -R(3)*R(3).zeta()/2)
meaningful_overlap(fourier_hurwitz(QQ(0),3), R(3).zeta())
meaningful_overlap(fourier_hurwitz(QQ(1)/2,3), -R(3)*R(3).zeta()/4)
assert abs(fourier_hurwitz(QQ(1)/2,2)) < R(2)**(-450)
print('CONTROL: gamma-ratio second and third moments; Fourier endpoints pass',flush=True)
coefficients = {n: second_coefficient(n) for n in independent_ns}
assert not agrees_with(coefficients, divisor_coefficient)
assert not agrees_with(coefficients, root_second_coefficient)
rows = []
values = {}
minimum_bits = 10000
for params in selection:
n,k = params['n'],params['k']
value = moment(n,k)
reference = R.pi()**2*divisor_coefficient(n) if k==2 else third_from_hurwitz(n)
meaningful_overlap(value,reference)
minimum_bits = min(minimum_bits,int(value.accuracy()))
text = to_text(value,digits=100)
# The whole certified computation ball must fit in the written interval.
assert decimal_contains(text,value)
values[(n,k)] = value
rows.append({'params':params,'number':text,'accuracy_bits':int(value.accuracy()),
'check_accuracy_bits':int(reference.accuracy()),
'independent_number':to_text(reference,digits=110)})
if k==3: print('CHECKED',n,'both moments, third sign',('positive' if value > 0 else 'negative' if value < 0 else 'undetermined'),flush=True)
assert not exactness(values)
# Exact polynomial identities underlying the reduction of every order.
P = PolynomialRing(ZZ,'x'); x = P.gen()
for n in range(1,421):
fac = ZZ(n).factor()
rad = ZZ(1)
for p,e in fac: rad *= p
assert P.cyclotomic_polynomial(n) == P.cyclotomic_polynomial(rad)(x**(n//rad))
for n in independent_ns:
sign = -1 if n==1 else 1
assert P.cyclotomic_polynomial(2*n) == sign*P.cyclotomic_polynomial(n)(-x)
if ZZ(n).is_prime():
assert P.cyclotomic_polynomial(n) == sum((x**j for j in range(n)),P(0))
print('Exact reduction identities through n=420, all stored odd representatives, and prime geometric sums pass',flush=True)
evidence={'table':'T681','entries':len(rows),'orders':independent_ns,
'minimum_accuracy_bits':minimum_bits,'measure':measure(values),'values':rows,
'skill_manifest':'agents/runs/20261102T162707Z-build-skills/manifest.json',
'skill_sha256':'a2b28d9953264d7b0ca62663785c2952190c5be6d6d3597fd695d9a54d606278'}
print('EVIDENCE_JSON '+json.dumps(evidence),flush=True)
if '--stored' in sys.argv or os.environ.get('CHECK_STORED')=='1':
doc=numberdb.table('T681')
records=doc['Numbers']
assert isinstance(records,list)
assert [(int(r['params']['n']),int(r['params']['k'])) for r in records] == list(values)
for r in records:
n,k=int(r['params']['n']),int(r['params']['k'])
reference=R.pi()**2*divisor_coefficient(n) if k==2 else third_from_hurwitz(n)
assert decimal_contains(r['number'],reference), (n,k)
report=generator.verify(sample=None)
print('STORED_VERIFY',report,flush=True)
assert report.ok
print('STORED: all written intervals contain the independent enclosures',flush=True)
if __name__=='__main__':
main()