back to table · edit · history · where entries came from · files · download
4483 bytes, as of the version from 2026-11-10 23:40 (current). Recorded here, not run.
"""Independent rational-breakpoint check for the Baez-Duarte R_N errors.
sage -python independent_check.py stored.json [result.json]
Input is a NumberDB document or an offline mapping N -> encoded interval.
This does not import generate.py. It builds Mobius values using the divisor
identity sum_{d|n} mu(d) = [n=1], builds the periodic step function from exact
jumps, and integrates its square and logarithmic cross term by endpoint
values of digamma and the deflated Hurwitz-zeta expansion at s=1.
See derivation.md. Shared with the generator: Arb, exact rationals, and the
mathematical definition. Different arithmetic construction, kernel and cross
term implementation. This is builder evidence, not an independent actor verdict.
"""
import json,sys
from functools import lru_cache
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ
from sage.rings.real_arb import RealBallField
from sage.rings.real_mpfi import RealIntervalField
from sage.rings.complex_arb import ComplexBallField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from decimal import Decimal
BITS=500
R=RealBallField(BITS);C=ComplexBallField(BITS)
P=PolynomialRing(C,'s');s=P.gen()
@lru_cache(None)
def endpoint(a):
z=(1+s)._zeta_series(2,a=C(a),deflate=True)
assert z[0].imag().contains_zero() and z[1].imag().contains_zero()
# z[0] = -psi(a); z[1] = -gamma_1(a)
return -z[0].real(),z[1].real()
def weights(n):
mu=[0]*(n+1)
for j in range(1,n+1):
mu[j]=(1 if j==1 else 0)-sum(mu[d] for d in range(1,j) if j%d==0)
M=[sum(mu[:j+1]) for j in range(n+1)]
w={k:QQ(M[n//k])/k for k in range(1,n) if M[n//k]}
assert sum((k*v for k,v in w.items()),QQ(0))==0
return w
def direct(n):
n=int(n);w=weights(n);events={QQ(0):QQ(0),QQ(1):QQ(0)}
for k,v in w.items():
for j in range(1,k):
a=QQ(j)/k;events[a]=events.get(a,QQ(0))-v
points=sorted(events);c=QQ(0);norm=R(0);minus_cross=R(0)
for a,b in zip(points,points[1:]):
c+=events[a]
if not c:continue
assert a>=QQ(1)/n
pa,za=endpoint(a);pb,zb=endpoint(b);delta=pb-pa
norm+=R(c*c)*delta/n
minus_cross+=R(c)*((R(n).log()+1)*delta+zb-za)/n
return 2+norm+2*minus_cross
def control3():
# F_3=-1/2 on [1/2,1); the digamma duplication identity gives
# norm=log(2)/6 and the elementary cancelled Mellin cross term below.
l2=R(2).log()
return 2+l2/6-2*l2*(1-R.euler_constant()+R(3).log()-l2/2)/3
def compare(encoded,value):
if not value.is_finite() or value.accuracy()<380:return False
dec=Decimal(str(encoded).rstrip('?'))
num,den=dec.as_integer_ratio(); center=QQ(num)/den
ulp=QQ(10)**int(dec.as_tuple().exponent)
return bool(value.lower()>=center-ulp and value.upper()<=center+ulp)
def run(input_path,result_path):
doc=json.loads(Path(input_path).read_text())
if 'Numbers' in doc:
rows={str(k):v['number'] for k,v in doc['Numbers'].items()}
else:rows=doc
result={'checked_entries':0,'total_entries':len(rows),'failures':[],
'known_answer_control':False,'negative_control_rejected':False,
'method':'rational-breakpoint integration with deflated Hurwitz zeta, 500-bit Arb',
'comparisons':[]}
v3=direct(3)
result['known_answer_control']=bool(direct(1)==2 and direct(2)==2 and (v3-control3()).contains_zero() and v3.accuracy()>380)
result['negative_control_rejected']=not compare('1.573641014547909618621335725814',v3)
result['invalid_ball_controls']=not compare('2.00',R('nan')) and not compare('2.00',R(RealIntervalField(BITS)(-100,100)))
for n,encoded in rows.items():
value=direct(int(n));ok=compare(encoded,value)
result['checked_entries']+=1
result['comparisons'].append({'N':int(n),'ok':ok,'ball':str(value),'accuracy_bits':int(value.accuracy())})
if not ok:result['failures'].append(int(n))
Path(result_path).write_text(json.dumps(result,indent=2)+'\n')
if int(n)%10==0:print('checked through',n,flush=True)
assert result['known_answer_control'] and result['negative_control_rejected'] and result['invalid_ball_controls'] and not result['failures'],result['failures']
print('Independent check:',result['checked_entries'],'passed; controls passed')
if __name__=='__main__':run(sys.argv[1],sys.argv[2] if len(sys.argv)>2 else 'check-results.json')