back to table · edit · history · where entries came from · files · download
9178 bytes, as of the version from 2026-11-01 08:20 (current). Recorded here, not run.
"""Independent numerical checks for Values of the Riesz function, T645.
Run: sage -python checks.py [--stored]
Keep generate.py, check.py, skill-manifest.json and T644.json beside this file.
The optional --stored uses the authenticated NumberDB client to read T645;
otherwise all results are checked before any write. Output is JSON evidence.
No website mutations are performed by this script.
The independent method splits the Mobius representation at N=16. It evaluates
n<=N as exponentials and expands the rest through k=99. The residual coefficient
is 1/zeta(2k+2)-sum_{n<=N}mu(n)/n^(2k+2). Taylor's integral remainder for e^-t
on t>=0 is bounded by t^K/K!, hence the omitted tail of R(x) is at most
x^(K+1)/(K!*(2K+1)*N^(2K+1)), K=100. Reciprocal zeta is recomputed through
exact Bernoulli numbers and Machin's pi identity, never Arb zeta. All independent
arithmetic uses 900-bit balls. This method shares no evaluation code with riesz.
The coefficient EGF is separately checked using finite differences at 1600 bits.
The bound |c_k|<=zeta(2)<2 follows from its Mobius representation. Thus 2 times
the absolute exponential-series majorant, multiplied by x*e^-x, bounds its tail.
"""
import json, sys, pathlib, hashlib, math
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 numberdb._write import to_text
from generate import RieszValues, riesz, truncation
from check import exactness, measure, names_its_rings, agrees_with
import sage.version
HERE = pathlib.Path(__file__).resolve().parent
F = RealBallField(900)
N, K = 16, 100
def decimal_rational(text):
text = str(text)
mant, *ex = text.lower().split('e')
exponent = int(ex[0]) if ex else 0
if '.' in mant:
whole, frac = mant.split('.')
return QQ(int(whole+frac)) * QQ(10)**(exponent-len(frac))
return QQ(int(mant)) * QQ(10)**exponent
def written_bounds(text):
text = str(text)
if '.' not in text and 'e' not in text.lower():
q = QQ(text); return q,q,QQ(0)
mant, *ex = text.lower().split('e')
exponent = int(ex[0]) if ex else 0
ulp = QQ(10)**(exponent - (len(mant.split('.')[1]) if '.' in mant else 0))
centre = decimal_rational(text)
return centre-ulp, centre+ulp, ulp
def enclosed(text, value):
lo, hi, ulp = written_bounds(text)
if ulp == 0:
return value == lo
return bool(value.is_finite() and F(lo) < value and value < F(hi))
def mobius(n):
sign, p = 1, 2
while p*p <= n:
if n % p == 0:
n //= p; sign = -sign
if n % p == 0:
return 0
p += 1
return -sign if n > 1 else sign
def atan_inverse(q):
# Alternating rational terms with an explicit next-term error.
total=F(0); j=0
target=QQ(1)/ZZ(10)**285
while True:
total += F(QQ((-1)**j) / ((2*j+1)*ZZ(q)**(2*j+1)))
j+=1
tail=QQ(1)/((2*j+1)*ZZ(q)**(2*j+1))
if tail < target:
return total + F(0,tail)
def bernoullis(last):
B=[QQ(1)]
for n in range(1,last+1):
B.append(-sum((QQ(math.comb(n+1,j))*B[j] for j in range(n)),QQ(0))/QQ(n+1))
return B
PI=16*atan_inverse(5)-4*atan_inverse(239)
B=bernoullis(2*K)
MU=[0]+[mobius(n) for n in range(1,N+1)]
assert MU[1:7] == [1,-1,-1,0,-1,1]
assert B[:5] == [QQ(1),-QQ(1)/2,QQ(1)/6,QQ(0),-QQ(1)/30]
assert PI.is_finite() and PI.accuracy()>880 and PI>3 and PI<4
RECIP=[2*F(math.factorial(2*k+2))/((2*PI)**(2*k+2)*F(abs(B[2*k+2]))) for k in range(K)]
RESIDUAL=[RECIP[k]-F(sum((QQ(MU[n])/ZZ(n)**(2*k+2) for n in range(1,N+1)),QQ(0))) for k in range(K)]
def independent(x):
x=QQ(x)
if x==0: return ZZ(0)
small=sum((F(QQ(MU[n])/n**2)*(-F(x)/n**2).exp() for n in range(1,N+1)),F(0))
tailpoly=F(0)
for k in reversed(range(K)):
tailpoly=tailpoly*(-F(x))+RESIDUAL[k]/math.factorial(k)
remainder=x**(K+1)/QQ(math.factorial(K)*(2*K+1)*N**(2*K+1))
return F(x)*(small+tailpoly)+F(0,remainder)
def egf_values(xs):
E=RealBallField(1600)
M,_=truncation(QQ(max(xs)),135)
diffs=[1/E(2*j+2).zeta() for j in range(M)]
coeff=[]
for j in range(M):
coeff.append(diffs[0]); diffs=[diffs[k]-diffs[k+1] for k in range(len(diffs)-1)]
results={}
for x in xs:
if x==0: results[x]=ZZ(0);continue
t,s=E(1),E(0)
for j,c in enumerate(coeff):
s+=t*c; t*=E(x)/(j+1)
bound=2*E(x)*t/(1-E(x)/(M+1))*(-E(x)).exp()
value=E(x)*(-E(x)).exp()*s+E(0,bound.upper())
assert value.is_finite() and value.accuracy()>450
results[x]=value
return results,coeff,M
def main():
# Controls precede the whole-range comparison.
assert enclosed('0.500000',F(QQ(1)/2))
assert not enclosed('0.500000',F(QQ(500010)/10**6))
assert (1/F(2).zeta()).overlaps(6/PI**2)
assert not (1/F(2).zeta()).overlaps(6/PI**2+F(QQ(1)/1000))
assert riesz(0)==0 and independent(0)==0
assert independent(1)>0 and independent(2)<0
generator=RieszValues(); params=list(generator.enumerate())
assert params==[{'x':str(x)} for x in range(101)]
values={int(p['x']):generator.value(p,100) for p in params}
assert not exactness(values)
assert not names_its_rings(str(HERE/'generate.py'))
egf, coefficients, M=egf_values(list(values))
# Independent Bernoulli inputs against the ingredients, not stored decimals.
assert all(RECIP[k].overlaps(1/F(2*k+2).zeta()) for k in range(K))
# Closed forms fix the signs and normalization of the EGF coefficients.
c0=6/PI**2; c1=c0-90/PI**4; c2=c0-180/PI**4+945/PI**6
assert all(coefficients[k].overlaps(v) for k,v in enumerate([c0,c1,c2]))
companion=json.loads((HERE/'T644.json').read_text())['Numbers']
for k in range(101):
assert enclosed(companion[str(k)],F(coefficients[k]))
if '--stored' in sys.argv:
document=numberdb.table('T645'); numbers=document['Numbers']
assert list(numbers)==[str(x) for x in range(101)]
else:
numbers={str(x):to_text(v,100) for x,v in values.items()}
rows=[]; independent_values={}
for x,v in values.items():
text=numbers[str(x)]
if isinstance(text,dict):text=text['number']
ind=independent(x); independent_values[x]=ind
assert enclosed(text,v) and enclosed(text,ind) and enclosed(text,egf[x]),x
if x:
assert v.is_finite() and v.accuracy()>350
assert ind.is_finite() and ind.accuracy()>600
assert v.overlaps(ind) and v.overlaps(egf[x])
# Cislo-Wolf (33), including the full lower-order terms.
k=F(x)
bound=3*PI.sqrt()/16/k**(QQ(3)/2)+F(27)/2*(-F(3)).exp()/k**2+15*PI.sqrt()/32/k**(QQ(5)/2)+128*(-F(4)).exp()/k**3
lo,hi,ulp=written_bounds(text)
stored_ball=F((lo+hi)/2,ulp)
clo,chi,cu=written_bounds(companion[str(x)])
ck=F((clo+chi)/2,cu)
assert abs(stored_ball/k-ck)<bound
# Every decimal is sensitive to a deliberate ten-ulp shift.
shifted=ind+10*F(ulp)
assert not enclosed(text,shifted)
rows.append({'x':str(x),'number':text,'primary_accuracy_bits':int(v.accuracy()) if x else None,'independent_accuracy_bits':int(ind.accuracy()) if x else None,'primary_radius':str(v.rad()) if x else '0','independent_radius':str(ind.rad()) if x else '0','terms':truncation(QQ(x),100)[0] if x else 0,'three_enclosures_inside_written_interval':True})
# agrees_with is intended for exact equality; apply it to independently
# rounded rational decimal centers (all 100 significant digits agree).
assert not agrees_with({x:decimal_rational(to_text(v,100)) for x,v in values.items()},lambda x:decimal_rational(to_text(independent_values[x],100)))
evidence={'table':'T645','stage':'build','stored_values_checked':'--stored' in sys.argv,'skill_manifest':'agents/runs/20261101T075931Z-build-skills/manifest.json','skills':json.loads((HERE/'skill-manifest.json').read_text()),'versions':{'sage':sage.version.version,'numberdb':numberdb.__version__},'entries':len(rows),'measurement_raw_balls':measure(values),'measurement_written':measure(numbers),'primary_bits_range':[numberdb.bits(100,losing=82),numberdb.bits(100,losing=280)],'minimum_primary_accuracy_bits':min(r['primary_accuracy_bits'] for r in rows if r['x']!='0'),'minimum_independent_accuracy_bits':min(r['independent_accuracy_bits'] for r in rows if r['x']!='0'),'independent_N':N,'independent_K':K,'independent_working_bits':900,'egf_terms':M,'egf_working_bits':1600,'checks':['exact origin','known rational nonzero control','deliberate-error controls','all 100 reciprocal zeta factors from exact Bernoulli numbers and Machin pi','all 101 accelerated Mobius values','all 101 EGF values with tail bound','first three coefficient closed forms','T644 coefficient intervals k=0..100','Cislo-Wolf equation (33) for k=1..100','every primary and independent enclosure contained in written last-place interval','all 100 ten-ulp perturbations rejected'], 'rows':rows}
print(json.dumps(evidence,indent=2))
if __name__=='__main__':main()