back to table · edit · history · where entries came from · files · download
10151 bytes, as of the version from 2026-11-01 08:46 (current). Recorded here, not run.
"""Independent checks for T646; run: sage -python checks.py.
No website writes. With --stored and NUMBERDB_API_KEY (or a stdin key via
NUMBERDB_KEY_FROM_STDIN=1), also check the saved decimal interval and all
257 stored coefficients of T644. Print the evidence JSON to stdout.
The main independent calculation uses the arithmetic kernel mu(n)/(2*n*n-1).
For n>A, expand its geometric series through m=L and evaluate reciprocal
zeta values from exact Bernoulli numbers and a Machin-series enclosure of pi.
No Arb zeta or pi function is used in that independent calculation.
"""
import json
import os
import sys
from math import comb, factorial
from decimal import Decimal
from pathlib import Path
import numberdb.sage as numberdb
from numberdb._write import to_text
from numberdb._generate import _stored_entries
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.version import version as sage_version
from generate import AlternatingSum
PRECISION = 700
SPLIT = 16
TERMS = 60
POLICY_MANIFEST = 'agents/runs/20261101T082247Z-build-skills/manifest.json'
POLICY_SHA256 = 'a2b28d9953264d7b0ca62663785c2952190c5be6d6d3597fd695d9a54d606278'
def bernoulli_list(n):
values = [QQ(1)]
for k in range(1, n+1):
values.append(-sum((QQ(comb(k+1, j))*values[j] for j in range(k)), QQ(0))/QQ(k+1))
return values
def machin_pi(field):
# The alternating arctangent series has the sign and magnitude of its
# first omitted term. Use a symmetric bound, with exact rationals.
def arctan_inverse(q):
total = QQ(0)
j = 0
while True:
term = QQ(1) / ((2*j+1)*q**(2*j+1))
if term < QQ(1)/2**760:
return field(total, term), j
total += (-1)**j * term
j += 1
a, na = arctan_inverse(5)
b, nb = arctan_inverse(239)
return 16*a - 4*b, [na, nb]
def mobius_sieve(limit):
mu = [1]*(limit+1)
isprime = [True]*(limit+1)
mu[0] = 0
for p in range(2,limit+1):
if not isprime[p]:
continue
for n in range(p,limit+1,p):
mu[n] *= -1
isprime[n] = False
for n in range(p*p,limit+1,p*p):
mu[n] = 0
return mu
def reciprocal_even_zeta(m, field, bernoulli, pi):
return field(QQ(2*factorial(2*m))/abs(bernoulli[2*m]))/(2*pi)**(2*m)
def independent_constant(field, bernoulli, pi, mu):
answer = sum((field(QQ(mu[n])/(2*n*n-1)) for n in range(1,SPLIT+1)), field(0))
for m in range(1,TERMS+1):
prefix = sum((QQ(mu[n])/n**(2*m) for n in range(1,SPLIT+1)), QQ(0))
answer += (reciprocal_even_zeta(m,field,bernoulli,pi)-field(prefix))/2**m
tail = (QQ(1)/(2**(TERMS+1)*(2*TERMS+1)*SPLIT**(2*TERMS+1))
/(1-QQ(1)/(2*(SPLIT+1)**2)))
return answer + field(0,tail), tail
def decimal_rational(text):
value = Decimal(str(text))
numerator, denominator = value.as_integer_ratio()
return QQ(numerator)/denominator
def written_interval(text, field):
centre = decimal_rational(text)
exponent = Decimal(text).as_tuple().exponent
unit = QQ(10)**exponent
return field(centre,unit), centre, unit
def inside_written(ball, text, field):
_, centre, unit = written_interval(text,field)
return bool(ball.is_finite() and ball > field(centre-unit) and ball < field(centre+unit))
def records(document):
# Use the package reader for both nested and flat API documents.
names=list(document.get('Parameters') or {})
result=[]
for identity, text in _stored_entries(document).items():
parts=identity.split(',') if names else []
assert len(parts)==len(names)
result.append({'params':dict(zip(names,parts)),'number':text})
return result
def run(stored=False):
field=RealBallField(PRECISION)
bernoulli=bernoulli_list(2*TERMS)
assert bernoulli[:5] == [QQ(1),-QQ(1)/2,QQ(1)/6,QQ(0),-QQ(1)/30]
assert bernoulli[12] == -QQ(691)/2730
pi, atan_terms=machin_pi(field)
assert pi.is_finite() and pi.accuracy()>680
assert pi.overlaps(field.pi())
# Exact Gaussian arithmetic verifies Machin's angle identity:
# (5+i)^4/(239+i) has equal positive real and imaginary parts.
z=(QQ(1),QQ(0))
for _ in range(4): z=(5*z[0]-z[1],z[0]+5*z[1])
z=(239*z[0]+z[1],239*z[1]-z[0])
assert z[0]==z[1] and z[0]>0
mu=mobius_sieve(10000)
assert mu[1:11]==[1,-1,-1,0,-1,1,-1,0,0,1]
divisor_sums=[0]*10001
for d in range(1,10001):
for multiple in range(d,10001,d): divisor_sums[multiple]+=mu[d]
assert divisor_sums[1]==1 and all(v==0 for v in divisor_sums[2:])
for m,expected in [(1,field(6)/pi**2),(2,field(90)/pi**4),(3,field(945)/pi**6)]:
assert reciprocal_even_zeta(m,field,bernoulli,pi).overlaps(expected)
generator=AlternatingSum()
assert list(generator.enumerate())==[{}]
primary=generator.value({},100)
extra=generator.value({},130)
independent,tail=independent_constant(field,bernoulli,pi,mu)
for ball in (primary,extra,independent):
assert ball.is_finite() and ball.accuracy()>380
assert primary.overlaps(extra) and primary.overlaps(independent)
text=to_text(primary,100)
assert inside_written(primary,text,field)
assert inside_written(independent,text,field)
assert inside_written(extra,text,field)
interval,centre,unit=written_interval(text,field)
assert not inside_written(independent+field(10*unit),text,field)
assert not inside_written(field(0,QQ(1)),text,field)
# A slow, unaccelerated arithmetic sum, with an entirely elementary
# remainder bound: sum(n>N) 1/(2n^2-1) <= sum(n>N) 1/n^2 < 1/N.
direct=sum((field(QQ(mu[n])/(2*n*n-1)) for n in range(1,10001)),field(0))
assert abs(direct-independent)<field(QQ(1)/10000)
# Partial alternating sums of fresh c_k and the finite arithmetic kernel
# must agree within the tail. No infinite rearrangement is used here.
max_k=256
reciprocal=[field(1)/field(2*j+2).zeta() for j in range(max_k+1)]
coeff=[]
for k in range(max_k+1):
coeff.append(sum(((-1)**j*comb(k,j)*reciprocal[j] for j in range(k+1)),field(0)))
partial_controls={}
for K in (0,1,8,32,256):
partial=sum(((-1)**k*coeff[k] for k in range(K+1)),field(0))
arithmetic=sum((field(QQ(mu[n])/(2*n*n-1))
*(1-(-field(1-QQ(1)/(n*n)))**(K+1))
for n in range(1,10001)),field(0))
bound=QQ(2)/10000
assert abs(partial-arithmetic)<field(bound)
partial_controls[str(K)]={'absolute_error_bound':str(bound),'passed':True}
stored_checks={}
if stored:
doc=numberdb.table('T646')
rows=records(doc)
assert len(rows)==1 and rows[0]['params']=={}
stored_text=str(rows[0]['number'])
assert inside_written(primary,stored_text,field)
assert inside_written(independent,stored_text,field)
assert inside_written(extra,stored_text,field)
report=generator.verify(sample=None)
assert report.ok
companion=records(numberdb.table('T644'))
assert len(companion)==257
stored_coeff={int(row['params']['k']):written_interval(str(row['number']),field)[0] for row in companion}
assert set(stored_coeff)==set(range(257))
companion_text={int(row['params']['k']):str(row['number']) for row in companion}
for k in range(257):
assert inside_written(coeff[k],companion_text[k],field)
for K in (0,1,8,32,256):
stored_partial=sum(((-1)**k*stored_coeff[k] for k in range(K+1)),field(0))
arithmetic=sum((field(QQ(mu[n])/(2*n*n-1))
*(1-(-field(1-QQ(1)/(n*n)))**(K+1))
for n in range(1,10001)),field(0))
assert abs(stored_partial-arithmetic)<field(QQ(2)/10000)
stored_checks={'tid':'T646','entries':1,'stored_text':stored_text,
'independent_enclosure_contained':True,'verify':str(report),
'T644_coefficients_checked':257,'stored_partial_sum_checks':partial_controls}
# A source digit prefix is an external check, not a certification method.
oeis_prefix='0.782527985325384234576688474283784076805800979652599159992647208782607376719900352292716'
assert text.startswith(oeis_prefix)
result={
'table':'T646','skill_manifest':POLICY_MANIFEST,'skill_local_sha256':POLICY_SHA256,
'skill_live_sha256':POLICY_SHA256,'entries':1,'parameters':[],
'value_100_digits':text,'primary':{'terms':400,'working_bits':numberdb.bits(100,losing=80),
'relative_accuracy_bits':int(primary.accuracy()),'ball':str(primary),
'tail_bound_exact':str(QQ(1)/2**400)},
'independent':{'method':'Möbius kernel; exact Bernoulli recurrence; Machin pi enclosure',
'working_bits':PRECISION,'split':SPLIT,'geometric_terms':TERMS,'atan_terms':atan_terms,
'relative_accuracy_bits':int(independent.accuracy()),'ball':str(independent),
'tail_bound_exact':str(tail)},
'higher_precision_relative_accuracy_bits':int(extra.accuracy()),
'primary_independent_and_130_digit_balls_inside_written_interval':True,
'controls':['Bernoulli B0..B4 and B12','Machin exact angle identity and pi enclosure',
'Mobius known values and divisor convolution through 10000','zeta(2), zeta(4), zeta(6)',
'ten-last-place perturbation rejected','radius-one ball rejected'],
'direct_mobius_terms':10000,'direct_mobius_tail_bound':'1/10000',
'partial_sum_checks':partial_controls,'oeis_prefix_digits':len(oeis_prefix)-2,
'stored_checks':stored_checks,'versions':{'sage':sage_version,'numberdb':numberdb.__version__}}
print(json.dumps(result,indent=2,ensure_ascii=False))
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
key=sys.stdin.read().strip()
if key.startswith('NUMBERDB_'): key=key.split('=',1)[1].strip().strip('\"\'')
numberdb.configure(api_key=key)
run(stored='--stored' in sys.argv or os.environ.get('NUMBERDB_CHECK_STORED')=='1')