"""Eigenvalues of the clamped plate -- numberdb.org/T488.

Run with SageMath (or passagemath-symbolics) and numberdb installed:
    sage -pip install numberdb mpmath        # once
    sage -python generate.py                 # verify every stored entry
    sage -python generate.py --fill          # fill the existing draft only
    sage -python checks.py                   # independent checks, offline
    sage -python checks.py --stored          # check the database values
Set NUMBERDB_API_KEY for authenticated access to the unpublished draft.
--fill sends entries, and does not offer, review or publish the table.

The first 30 eigenvalues on the radius-one disk, counted with multiplicity,
are the fourth powers of cross-product roots below 51/4. This cutoff falls
between complete eigenspaces. All roots are certified with Arb signs at exact
rational endpoints. mpmath supplies guesses only. See certificates.md for the
root-count argument and checks.py for the independent radial ODE computation.
The 180-bit guard leaves more than 110 decimal digits in every returned ball.
"""
import os
import sys
from functools import lru_cache
from fractions import Fraction

import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
import mpmath as mp

WORKING_GUARD = 180
CUTOFF = QQ(51)/4
MAX_ORDER = 9


def rational_decimal(text):
    f = Fraction(str(text))
    return QQ(f.numerator)/f.denominator


def sign(x):
    assert x.is_finite(), 'nonfinite sign enclosure'
    if x > 0:
        return 1
    if x < 0:
        return -1
    raise ArithmeticError('sign not certified')


def bessel(C, kind, m, x):
    z = getattr(C(x), 'bessel_' + kind)(m)
    assert z.real().is_finite() and z.imag().is_finite() and z.imag().contains_zero()
    return z.real()


def cross(C, m, x):
    return (bessel(C, 'J', m, x)*bessel(C, 'I', m+1, x)
            + bessel(C, 'J', m+1, x)*bessel(C, 'I', m, x))


def bisect_certified(f, left, right, width):
    sl, sr = sign(f(left)), sign(f(right))
    assert sl != sr
    while right-left > width:
        mid = (left+right)/2
        sm = sign(f(mid))
        if sm == sl:
            left = mid
        else:
            right = mid
    return left, right


def j_brackets(C, m):
    # There are no J_m zeros in (0,1]. On [1,infinity), sqrt(x) J_m(x)
    # has zero spacing > 2 for every integer m>=0 (Sturm comparison).
    # Hence a unit mesh and certified signs count *all* the zeros to 18.
    brackets = []
    left = QQ(1)
    previous = sign(bessel(C, 'J', m, left))
    assert previous == 1
    for stop in range(2, 19):
        right = QQ(stop)
        current = sign(bessel(C, 'J', m, right))
        if current != previous:
            brackets.append(bisect_certified(
                lambda x: bessel(C, 'J', m, x), left, right, QQ(1)/10**14))
        left, previous = right, current
    assert brackets and brackets[-1][0] > CUTOFF
    return brackets


def approximate_root(m, left, right, working):
    with mp.workdps(working):
        a = mp.mpf(str(left.numerator()))/mp.mpf(str(left.denominator()))
        b = mp.mpf(str(right.numerator()))/mp.mpf(str(right.denominator()))
        def f(x):
            return mp.besselj(m+1,x)+mp.besselj(m,x)*mp.besseli(m+1,x)/mp.besseli(m,x)
        root = mp.findroot(f, (a+(b-a)/5, a+2*(b-a)/5))
        assert a < root < b
        return mp.nstr(root, working-10)


@lru_cache(maxsize=None)
def spectrum(digits=100):
    C = ComplexBallField(numberdb.bits(digits, losing=WORKING_GUARD))
    R = RealBallField(C.precision())
    modes = []
    evidence = []
    eps = QQ(1)/10**(digits+12)
    for m in range(MAX_ORDER+1):
        js = j_brackets(C, m)
        for k in range(1, len(js)):
            # A root is between the kth and (k+1)st J_m zeros. If the kth
            # zero already exceeds the cutoff, all later roots do too.
            if js[k-1][0] > CUTOFF and k > 1:
                break
            a, b = js[k-1][1], js[k][0]
            q = rational_decimal(approximate_root(m, a, b, digits+40))
            lo, hi = q-eps, q+eps
            assert a < lo < hi < b
            assert sign(cross(C,m,lo)) != sign(cross(C,m,hi))
            root = R((lo+hi)/2) + R(0,(hi-lo)/2)
            derivative = (2*bessel(C,'J',m,root)*bessel(C,'I',m,root)
                          - cross(C,m,root)/root)
            sign(derivative)
            assert root.is_finite() and root.rad() < abs(root.center())/2**350
            evidence.append((m,k,lo,hi,js[k-1],js[k]))
            if lo > CUTOFF:
                break
            assert hi < CUTOFF
            eigenvalue = root**4
            modes.append((eigenvalue,m,k))
        if m == MAX_ORDER:
            assert evidence[-1][0:2] == (MAX_ORDER,1)
            assert evidence[-1][2] > CUTOFF
    # gamma_{m,k} increases with m; m=9's first root excludes all m>=9.
    modes.sort(key=lambda row: row[0].center())
    for a,b in zip(modes,modes[1:]):
        assert a[0] < b[0], 'distinct mode enclosures overlap'
    rows=[]
    for value,m,k in modes:
        rows.extend([(value,m,k)]*(1 if m==0 else 2))
    assert len(rows)==30
    return rows,evidence


class ClampedPlate(numberdb.Generator):
    table = 'T488'
    files = ('generate.py', 'checks.py', 'certificates.md')
    parameters = ('domain','n')
    type = 'R'
    digits = 100
    rigour = 'proven'

    def enumerate(self):
        for n in range(1,31):
            yield {'domain':'unit-disk','n':n}

    def value(self,params,digits):
        assert params['domain']=='unit-disk'
        n=int(params['n'])
        if not 1<=n<=30:
            raise ValueError('This generator covers the first 30 modes.')
        value,m,k=spectrum(digits)[0][n-1]
        result={'number':value}
        if n==1:
            result['comment']='The fundamental mode is axisymmetric.'
        if n>1 and spectrum(digits)[0][n-2][1:]==(m,k):
            result['equals']='HREF{Eigenvalues_of_the_clamped_plate#unit-disk,%d}'%(n-1)
        return result


if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        numberdb.configure(api_key=sys.stdin.read().strip())
    generator=ClampedPlate()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
        print(generator.publish(message='Fill 30 certified clamped-disk eigenvalues'))
    else:
        report=generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)