back to table · edit · history · where entries came from · files · download
4231 bytes, as of the version from 2026-10-31 22:58 (current). Recorded here, not run.
"""Selberg approximation errors in the Nyman–Beurling criterion -- numberdb.org/T642.
Run with SageMath (or compatible passagemath) and the numberdb package:
sage -pip install numberdb
sage -python generate.py # verify the stored draft/table
sage -python generate.py --fill # write entries to an existing draft
For an unpublished draft, set NUMBERDB_API_KEY to an authorized key.
--fill writes numbers and attaches this source; it does not publish or offer.
The sole parameter N is the integer truncation index, N >= 2. Default coverage
is every N from 2 through 30, in order. See argument-selection.md. The squared
L2((0,infinity), dx) residual uses a PLUS sign, since the Selberg approximation
approximates -chi. Weights at n=N and at nonsquarefree n vanish exactly.
The Gram kernel uses the finite Vasyunin-sum identity (BBLS, Proposition 89).
All elementary functions are evaluated as Arb enclosures at bits(digits)+128.
No stored decimal from a companion table is treated as an exact input.
checks.py independently integrates the fractional-part products and verifies
all errors and the least-squares comparison. See derivation.md for proofs.
"""
import functools
import os
import sys
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
MAX_N = 30
WORKING_GUARD = 128
@functools.lru_cache(None)
def mobius_values(limit):
"""Exact sieve; a separate prime-factorization implementation checks it."""
mu = [1] * (limit + 1)
prime = [True] * (limit + 1)
for p in range(2, limit + 1):
if prime[p]:
for n in range(p, limit + 1, p):
prime[n] = False
mu[n] *= -1
for n in range(p * p, limit + 1, p * p):
mu[n] = 0
return tuple(mu)
@functools.lru_cache(None)
def kernel(limit, bits):
R = RealBallField(bits)
C = (2*R.pi()).log() - R.euler_constant()
@functools.lru_cache(None)
def V(p,q):
return sum((R(QQ((m*p) % q)/QQ(q)) *
(R.pi()*R(QQ(m)/QQ(q))).cot()
for m in range(1,q)),R(0))
@functools.lru_cache(None)
def A(r):
p,q = int(r.numerator()),int(r.denominator())
x = R(r)
return ((1-x)*x.log()+(1+x)*C)/2 - R.pi()*(V(p,q)+V(q,p))/(2*q)
G = tuple(tuple(A(QQ(j)/QQ(k))/j for k in range(1,limit+1))
for j in range(1,limit+1))
b = tuple((R(n).log()+1-R.euler_constant())/n for n in range(1,limit+1))
return R,G,b
def weights(N,R):
mu = mobius_values(N)
logN = R(N).log()
return tuple(R(0) if n == N or mu[n] == 0 else
R(mu[n])*(1-R(n).log()/logN) for n in range(1,N+1))
def squared_error(N,digits=100):
if N < 2:
raise ValueError('The Selberg taper requires N >= 2')
R,G,b = kernel(max(N,MAX_N),numberdb.bits(digits,losing=WORKING_GUARD))
w = weights(N,R)
linear = sum((w[i]*b[i] for i in range(N)),R(0))
quadratic = sum((w[i]*w[i]*G[i][i] for i in range(N)),R(0))
quadratic += 2*sum((w[i]*w[j]*G[i][j] for i in range(N) for j in range(i)),R(0))
result = 1+2*linear+quadratic
assert result.is_finite() and result > 0
return result
class SelbergErrors(numberdb.Generator):
table = 'T642'
parameters = ('N',)
type = 'R'
digits = 100
rigour = 'proven'
def enumerate(self,maximum=MAX_N):
for N in range(2,maximum+1):
yield {'N':str(N)}
def value(self,params,digits):
return squared_error(int(params['N']),digits)
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
token = sys.stdin.read().strip()
if '=' in token and token.split('=',1)[0].isupper():
token = token.split('=',1)[1].strip().strip('\'"')
numberdb.configure(api_key=token)
generator = SelbergErrors()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(generator.publish(message='Fill all 29 Selberg squared errors with 100 certified digits'))
else:
report=generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)