back to table · edit · history · where entries came from · files · download
8852 bytes, as of the version from 2026-10-31 21:13 (current). Recorded here, not run.
"""Independent numerical checks for T641; run sage -python checks.py.
Requires generate.py alongside this file. Optional stored.json contains the
API document read back from T641, in which case each written interval is checked.
The independent kernel integrates fractional parts over one period; it uses
neither cotangents nor any numerical kernel code from generate.py. Gaussian
elimination independently recomputes every minimum. See derivation.md for the
analytic integral and remainder bound. JSON evidence is written to stdout.
"""
import json
import os
import platform
from pathlib import Path
from math import gcd
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.version import version as sage_version
import generate
PRECISION=720
K=16
M=125
R=RealBallField(PRECISION)
def require_close(a,b,threshold=110):
assert a.is_finite() and b.is_finite()
assert a.rad_as_ball()<R(10)**(-threshold)
assert b.rad_as_ball()<R(10)**(-threshold)
assert a.overlaps(b), (str(a),str(b))
def segments(p,q):
cuts=sorted(set([QQ(i)/q for i in range(q+1)]+
[QQ(i)/p for i in range(p+1)]))
out=[]
for left,right in zip(cuts,cuts[1:]):
mid=(left+right)/2
a=(q*mid).floor();b=(p*mid).floor()
out.append((left,right,QQ(p*q),QQ(-q*b-p*a),QQ(a*b)))
return out
def integral_rational_kernel(parts,n):
total=R(0)
for left,right,alpha,beta,delta in parts:
# Integrate P(u)/(n+u)^2 exactly apart from enclosed logarithms.
if n==0 and left==0:
assert beta==0 and delta==0
total+=R(alpha*right)
continue
v0=R(n+left); v1=R(n+right)
total+=(R(alpha*(right-left)) + R(beta-2*alpha*n)*(v1/v0).log()
+R(alpha*n*n-beta*n+delta)*(1/v0-1/v1))
return total
@lru_cache(None)
def zeta_tail(m):
return R(m+2).zeta(R(K))
@lru_cache(None)
def integral_A(p,q):
assert gcd(p,q)==1
parts=segments(p,q)
total=sum((integral_rational_kernel(parts,n) for n in range(K)),R(0))
# Exact moments of the polynomial on each open interval.
powers=[(left,right,alpha,beta,delta,left,right)
for left,right,alpha,beta,delta in parts]
for m in range(M+1):
moment=QQ(0)
for idx,(left,right,alpha,beta,delta,lp,rp) in enumerate(powers):
moment+=(alpha*(rp*right**2-lp*left**2)/(m+3)
+beta*(rp*right-lp*left)/(m+2)
+delta*(rp-lp)/(m+1))
powers[idx]=(left,right,alpha,beta,delta,lp*left,rp*right)
total+=((-1)**m)*(m+1)*zeta_tail(m)*R(moment)
# First omitted term bounds the alternating remainder pointwise.
# zeta(M+3,K) <= K^(-M-3)+K^(-M-2)/(M+2).
error=QQ(M+2)/QQ(K)**(M+3)+QQ(1)/QQ(K)**(M+2)
total=total.add_error(R(error))
return total/q
def gauss_distance(G,b,n):
A=[G[i][:n]+[b[i]] for i in range(n)]
for k in range(n):
assert A[k][k]>0
for i in range(k+1,n):
factor=A[i][k]/A[k][k]
for j in range(k+1,n+1):
A[i][j]-=factor*A[k][j]
A[i][k]=R(0)
x=[R(0) for i in range(n)]
for i in range(n-1,-1,-1):
x[i]=(A[i][n]-sum((A[i][j]*x[j] for j in range(i+1,n)),R(0)))/A[i][i]
return 1-sum((b[i]*x[i] for i in range(n)),R(0)),x
def decimal_interval(text):
"""Interpret a NumberDB decimal as exact centre +/- one last-place unit."""
from decimal import Decimal
d=Decimal(str(text)); sign,digits,exponent=d.as_tuple()
mantissa=int(''.join(str(v) for v in digits))*(-1 if sign else 1)
unit=QQ(10)**exponent
center=mantissa*unit
return center-unit,center+unit
def main():
# Known-answer nonzero controls before using the integration machinery.
parts=[(QQ(0),QQ(1),QQ(1),QQ(0),QQ(0))]
require_close(integral_rational_kernel(parts,0),R(1))
require_close(integral_rational_kernel(parts,1),R(3)/2-2*R(2).log())
C=(2*R.pi()).log()-R.euler_constant()
control=integral_A(1,1)
require_close(control,C)
assert not (control+R(QQ(1)/10**80)).overlaps(C)
# Linear-solver control: b=(1,2), G=diag(2,8) gives distance zero.
control_distance,control_x=gauss_distance([[R(2),R(0)],[R(0),R(8)]],[R(1),R(2)],2)
assert control_distance.contains_zero()
assert (control_x[0]-QQ(1)/2).contains_zero() and (control_x[1]-QQ(1)/4).contains_zero()
field,G,b,values,coeffs,pivots,L,y=generate.computation(30,100)
GI=[[R(0) for j in range(30)] for i in range(30)]
for j in range(1,31):
for k in range(1,j+1):
common=gcd(j,k)
GI[j-1][k-1]=GI[k-1][j-1]=integral_A(j//common,k//common)/j
require_close(G[j-1][k-1],GI[j-1][k-1])
bi=[(R(n).log()+1-R.euler_constant())/n for n in range(1,31)]
independent=[]; max_residual=R(0); max_quadratic_radius=R(0)
all_rows=[]
for n in range(1,31):
independent_value,x=gauss_distance(GI,bi,n)
require_close(values[n-1],independent_value)
independent.append(independent_value)
c=coeffs[n-1]
for i in range(n):
residual=sum((G[i][j]*c[j] for j in range(n)),field(0))-b[i]
assert residual.contains_zero() and residual.is_finite()
assert residual.abs()<field(10)**(-110)
if residual.abs().upper()>max_residual.upper(): max_residual=R(residual.abs().upper())
require_close(c[i],x[i])
quadratic=1-2*sum((b[i]*c[i] for i in range(n)),field(0))
quadratic+=sum((c[i]*G[i][j]*c[j] for i in range(n) for j in range(n)),field(0))
require_close(values[n-1],quadratic)
if quadratic.rad_as_ball()>max_quadratic_radius: max_quadratic_radius=R(quadratic.rad_as_ball())
if n>1:
assert values[n-2]>values[n-1]
require_close(values[n-2]-values[n-1],y[n-1]**2/pivots[n-1])
all_rows.append({'N':n,'distance_ball':str(values[n-1]),
'independent_distance_ball':str(independent_value),
'pivot_ball':str(pivots[n-1]),
'coefficients':[str(v) for v in c]})
known_one=1-(1-R.euler_constant())**2/C
require_close(values[0],known_one)
g12=3*C/4-R(2).log()/4
known_two=1-(C*bi[1]**2+(C/2)*bi[0]**2-2*g12*bi[0]*bi[1])/(C*C/2-g12*g12)
require_close(values[1],known_two)
stored_checks=0
stored_path=Path(__file__).with_name('stored.json')
if stored_path.exists():
document=json.loads(stored_path.read_text())
numbers=document['Numbers']
assert list(map(int,numbers))==list(range(1,31))
for n in range(1,31):
text=numbers[str(n)]
if isinstance(text,dict): text=text['number']
lo,hi=decimal_interval(text)
for ball in (values[n-1],independent[n-1]):
assert ball>lo and ball<hi
if n>1:
prev=numbers[str(n-1)]
if isinstance(prev,dict): prev=prev['number']
prev_lo,prev_hi=decimal_interval(prev)
assert prev_lo>hi
stored_checks+=1
skill_path=Path(__file__).with_name('skill-evidence.json')
skill=json.loads(skill_path.read_text()) if skill_path.exists() else None
report={'table':'T641','checks_passed':True,'count':30,
'controls':['integral of u^2/u^2 = 1','integral of u^2/(1+u)^2 = 3/2-2 log 2',
'integral A(1) = log(2 pi)-gamma','perturbed integral rejected',
'diagonal solve gives x=(1/2,1/4), distance 0'],
'gram_entries_checked':465,'distinct_independent_integrals':integral_A.cache_info().currsize,
'all_prefix_independent_solves':30,'all_coefficients_checked':465,
'positive_pivots':30,'positive_distances':30,'strict_decreases':29,
'quadratic_norms_checked':30,'schur_decreases_checked':29,
'stored_intervals_contain_both_enclosures':stored_checks,
'generator_working_bits':int(field.precision()),'independent_working_bits':PRECISION,
'tail_start':K,'tail_last_power':M,
'minimum_generator_accuracy_bits':min(int(v.accuracy()) for v in values),
'maximum_generator_radius':str(max(v.rad() for v in values)),
'maximum_independent_radius':str(max(v.rad() for v in independent)),
'maximum_normal_equation_residual_bound':str(max_residual),
'maximum_quadratic_radius':str(max_quadratic_radius),
'versions':{'sage':sage_version,'python':platform.python_version(),
'numberdb':getattr(numberdb,'__version__','unknown')},
'skills':skill,'rows':all_rows}
print(json.dumps(report,indent=2))
if __name__=='__main__': main()