back to table · edit · history · where entries came from · files · download
3328 bytes, as of the version from 2026-11-10 20:59 (current). Recorded here, not run.
"""Integral kernel reused from T641 checks.py, with its derivation rechecked.
Original SHA-256: 3a6ea0c72d6588b2fdda2d91a0414232f1d64092a5e928cd4121cc4d05798678.
No imports from the generating cotangent implementation.
"""
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
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