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