back to table · edit · history · where entries came from · files · download
5466 bytes, as of the version from 2026-10-26 15:32 (current). Recorded here, not run.
"""Slepian concentration eigenvalues -- numberdb.org/T519.
Run with SageMath and mpmath:
sage -pip install numberdb mpmath
sage -python generate.py # verify all entries, authenticated for draft
sage -python checks.py # independent numerical checks
sage -python generate.py --fill # fill the existing draft; API key required
The values are the decreasing eigenvalues of sin(c(x-y))/(pi(x-y))
on [-1,1], where c=pi*S/2. We retain 100 significant digits after requiring
110 digits of relative agreement. The Legendre truncation is tested with
80, 120 and 180 terms in each parity block, and arithmetic with 160 and 220
decimal digits. No rigorous bound on infinite-dimensional truncation is claimed.
The independent Gauss-Legendre Nystrom computation is in checks.py.
"""
import functools
import os
import sys
import numberdb.sage as numberdb
import mpmath as mp
SHANNON_NUMBERS = (1, 2, 4, 8)
MODES = 12
SETTINGS = ((160, 80), (220, 80), (220, 120), (220, 180))
AGREEMENT_DIGITS = 110
def differential_block(ctx, c, parity, size):
"""Matrix of -((1-x^2)f')'+c^2*x^2*f in orthonormal Legendre basis."""
def a(k):
return ctx.mpf(k + 1) / ctx.sqrt((2*k + 1)*(2*k + 3))
diagonal, off = [], []
for j in range(size):
ell = 2*j + parity
diagonal.append(ell*(ell+1) + c*c*(a(ell)**2 + (a(ell-1)**2 if ell else 0)))
if j + 1 < size:
off.append(c*c*a(ell)*a(ell+1))
return diagonal, off
def sturm_count(ctx, diagonal, off, x):
pivot = diagonal[0] - x
count = int(pivot < 0)
tiny = ctx.power(10, -2*ctx.dps)
for j in range(1, len(diagonal)):
if not pivot:
pivot = -tiny
pivot = diagonal[j] - x - off[j-1]**2 / pivot
count += int(pivot < 0)
return count
def mode_vector(ctx, c, n, diagonal, off):
lower = ctx.mpf(n*(n+1))
upper = lower + c*c
target = n//2
tolerance = ctx.power(10, -(ctx.dps-15))
while upper - lower > tolerance:
middle = (lower + upper)/2
if sturm_count(ctx, diagonal, off, middle) <= target:
lower = middle
else:
upper = middle
chi = (lower+upper)/2
# Backward recurrence selects the decaying Legendre tail. The scale
# cancels from the finite-Fourier eigenvalue ratio at the origin.
size = len(diagonal)
beta = [ctx.zero]*size
beta[-1] = ctx.one
for j in range(size-1, 0, -1):
beta[j-1] = ((chi-diagonal[j])*beta[j]
- (off[j]*beta[j+1] if j+1 < size else 0))/off[j-1]
scale = max(abs(b) for b in beta)
return chi, [b/scale for b in beta]
@functools.lru_cache(maxsize=None)
def spectrum(S, dps, terms):
ctx = mp.mp.clone()
ctx.dps = dps
c = ctx.pi*ctx.mpf(S)/2
result = {}
for parity in (0, 1):
diagonal, off = differential_block(ctx, c, parity, terms)
# P_(2j)(0) and P'_(2j+1)(0)=(2j+1)P_(2j)(0).
origin = []
p0 = ctx.one
for j in range(terms):
ell = 2*j + parity
origin.append(ctx.sqrt(ctx.mpf(2*ell+1)/2)*p0*(ell if parity else 1))
p0 *= -ctx.mpf(2*j+1)/(2*j+2)
for n in range(parity, MODES, 2):
chi, beta = mode_vector(ctx, c, n, diagonal, off)
denominator = ctx.fsum(b*p for b, p in zip(beta, origin))
if parity:
mu_abs = abs(c*ctx.sqrt(ctx.mpf(2)/3)*beta[0]/denominator)
else:
mu_abs = abs(ctx.sqrt(2)*beta[0]/denominator)
lam = ctx.mpf(S)*mu_abs**2/4
assert ctx.isfinite(lam) and 0 < lam < 1
result[n] = ctx.nstr(lam, dps-20)
assert all(ctx.mpf(result[n]) > ctx.mpf(result[n+1]) for n in range(MODES-1))
return result
def agreement_digits(a, b):
ctx = mp.mp.clone()
ctx.dps = 240
a, b = ctx.mpf(a), ctx.mpf(b)
if a == b:
return ctx.inf
return -ctx.log10(abs(a-b)/max(abs(a),abs(b)))
def rounded(text, digits):
ctx = mp.mp.clone()
ctx.dps = 240
return ctx.nstr(ctx.mpf(text), digits, strip_zeros=False)
class SlepianConcentration(numberdb.Generator):
table = 'T519'
parameters = ('S', 'n')
type = 'R'
digits = 100
rigour = 'heuristic (agreement-checked)'
files = ('generate.py', 'checks.py', 'evidence.json', 'build-report.md')
def enumerate(self):
for S in SHANNON_NUMBERS:
for n in range(MODES):
yield {'S': str(S), 'n': str(n)}
def value(self, params, digits):
if digits > 100:
raise ValueError('This measured configuration supports at most 100 stored digits.')
S, n = int(params['S']), int(params['n'])
values = [spectrum(S, dps, terms)[n] for dps, terms in SETTINGS]
for a, b in zip(values, values[1:]):
assert agreement_digits(a, b) >= AGREEMENT_DIGITS, (S, n, 'precision/truncation disagreement')
return rounded(values[-1], digits)
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
generator = SlepianConcentration()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(generator.publish(message='fill 48 concentration eigenvalues with 100 agreement-checked digits'))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)