"""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)