back to table · edit · history · where entries came from · files · download
2716 bytes, as of the version from 2026-10-27 02:15 (current). Recorded here, not run.
"""Legendre spectral helpers for PSWF zeros -- numberdb.org/T534.
Matrix and Sturm bisection adapted from the attached generator of T518,
https://numberdb.org/files/T518/generate.py?raw=1, downloaded 2026-10-27.
The differential eigenvalue convention is identical. No stored eigenvalues
are used: bandwidth and eigenvalues are recomputed at working precision.
"""
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_mpfr import RealField
BANDWIDTHS = (1, 2, 4, 8)
MODES = 12
SETTINGS = ((160, 80), (220, 80), (220, 120), (220, 180))
REQUIRED_AGREEMENT = 110
def field(dps):
return RealField(numberdb.bits(dps, losing=16))
@lru_cache(maxsize=None)
def rational_matrix(parity, size):
diagonal, off_squared = [], []
for j in range(size):
k = 2*j + parity
x2 = QQ((k+1)**2) / ((2*k+1)*(2*k+3))
x2 += QQ(k*k) / ((2*k-1)*(2*k+1))
diagonal.append((k*(k+1), x2))
if j < size-1:
off_squared.append(QQ((k+1)**2 * (k+2)**2) /
((2*k+1)*(2*k+3)**2*(2*k+5)))
return tuple(diagonal), tuple(off_squared)
def matrix(q, parity, size):
R = q.parent()
diagonal, off_squared = rational_matrix(parity, size)
return ([R(k) + q*R(a) for k, a in diagonal],
[q*q*R(b) for b in off_squared])
def sturm_count(diagonal, off_squared, x):
"""Number of eigenvalues strictly below x (zero pivots by a limit)."""
pivot = diagonal[0] - x
count = int(pivot < 0)
tiny = x.parent()(2)**(-x.parent().precision()*2)
if not pivot:
pivot = -tiny
count = 1
for i, square in enumerate(off_squared):
pivot = diagonal[i+1] - x - square/pivot
if not pivot:
pivot = -tiny
count += int(pivot < 0)
return count
def bisect_mode(q, n, size, dps):
R = q.parent()
diagonal, off_squared = matrix(q, n % 2, size)
index = n // 2
lo = R(n*(n+1))-1
hi = R(n*(n+1))+q+1
assert sturm_count(diagonal, off_squared, lo) <= index
assert sturm_count(diagonal, off_squared, hi) > index
tolerance = R(10)**(-(dps-10))
while hi-lo > tolerance:
mid = (lo+hi)/2
if sturm_count(diagonal, off_squared, mid) <= index:
lo = mid
else:
hi = mid
return (lo+hi)/2
@lru_cache(maxsize=None)
def spectrum(S, dps, size):
R = field(dps)
c = R.pi()*R(QQ(S))/2
return tuple(bisect_mode(c*c, n, size, dps) for n in range(MODES))
def decimal_string(x, digits):
# Fixed significant digits, retaining a decimal point, including final zeros.
return x.str(truncate=False, digits=digits, skip_zeroes=False)