"""Dimensionless squared angular momentum at axisymmetric bifurcations of Maclaurin spheroids -- numberdb.org/T630.

Run with SageMath (or passagemath) and numberdb:
    sage -pip install numberdb
    sage -python generate.py
    sage -python generate.py --fill-draft

The default verifies all stored entries. Filling requires a NumberDB API key
and sends entries only; it does not publish or offer a draft. The enumeration
is k=1,...,10. See argument-selection.md for the selection and certification.md
for the exact root-uniqueness proof and error propagation. The root certificate
is adapted from the public T627 generator, fetched from the website; this
standalone file recomputes all roots and does not use stored eccentricities.

At 100 digits the computation uses 500 bits with a 160-bit minimum guard and
422 bisections. Bounds, second-precision agreement and independent integral
checks are recorded in numerical-evidence.json. Each failed certificate stops
the calculation. The same algorithm can extend the index only while its exact
polynomial sign certificates continue to pass.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.real_arb import RealBallField

WORKING_GUARD = 160
MIN_BITS = 500

@lru_cache(None)
def polynomials(k):
    k = int(k)
    if k < 1:
        raise ValueError('The branch index must be positive.')
    ring = PolynomialRing(QQ, 'x')
    x = ring.gen()
    p = [ring(1), x]
    r = [ring(0), ring(1)]
    ell = 2*k+2
    for n in range(1, ell):
        p.append(((2*n+1)*x*p[n] + n*p[n-1])/QQ(n+1))
        r.append(((2*n+1)*x*r[n] + n*r[n-1])/QQ(n+1))
    p, r = p[ell], r[ell]
    # d(q/p)/dx = -1/((1+x^2)p^2), q=p atan(1/x)-r.
    assert (1+x*x)*(r*p.derivative()-r.derivative()*p)-p*p == -1
    assert all(c >= 0 for c in p.list()) and p[0] > 0
    a, b = p*p+x*x, p*r+x
    numerator = (1+x*x)*(b*a.derivative()-b.derivative()*a)-a*a
    # H'=numerator/((1+x^2)A^2). Its numerator has one positive zero.
    assert numerator[0] < 0
    assert all(c >= 0 for c in numerator.list()[1:])
    assert numerator.leading_coefficient() > 0
    assert b.degree() < a.degree()
    return p, r, a, b, numerator


def h(k, x):
    _, _, a, b, _ = polynomials(k)
    if x == 0:
        return x.parent().pi()/2
    return (1/x).arctan() - b(x)/a(x)


def enclose(k, digits=100, working_bits=None):
    field = RealBallField(working_bits or max(MIN_BITS, numberdb.bits(digits, losing=WORKING_GUARD)))
    lo, hi = QQ(0), QQ(1)
    assert h(k, field(lo)) > 0 and h(k, field(hi)) < 0
    # H(0)>0; H decreases then increases to zero. The root is unique on x>0.
    steps = numberdb.bits(digits+22)
    for _ in range(steps):
        mid = (lo+hi)/2
        v = h(k, field(mid))
        if v > 0:
            lo = mid
        elif v < 0:
            hi = mid
        else:
            raise ArithmeticError('Increase working precision: unresolved sign.')
    xi = field((lo+hi)/2).add_error(field((hi-lo)/2))
    e = 1/(1+xi*xi).sqrt()
    assert e.is_finite() and e > 0 and e < 1
    assert field(e.rad()) < field(10)**(-digits-18)
    theta = (1/xi).arctan()
    omega2 = (xi*(1+3*xi*xi)*theta-3*xi*xi)/2
    j2 = QQ(4)/25 * (3/(4*field.pi()))**(QQ(4)/3) * omega2 * ((1+xi*xi)/(xi*xi))**(QQ(2)/3)
    assert j2.is_finite() and j2 > 0
    assert field(j2.rad()) < field(10)**(-digits-18)
    return j2, lo, hi


class MaclaurinAngularMomentum(numberdb.Generator):
    table = 'T630'
    parameters = ('k',)
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'argument-selection.md', 'certification.md',
             'numerical-evidence.json', 'skill-manifest.json', 'source-provenance.json')

    def enumerate(self, limit=10):
        for k in range(1, int(limit)+1):
            yield {'k': k}

    def value(self, params, digits):
        k = int(params['k'])
        value = enclose(k, digits)[0]
        if k == 1:
            return {'number': value, 'comment': 'This branch continues to the Dyson-ring sequence CITE{AKM}.'}
        return value


def client_from_stdin():
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key = sys.stdin.read().strip()
        if '=' in key and key.split('=',1)[0].isupper():
            key = key.split('=',1)[1].strip().strip('\'"')
        return numberdb.Client(api_key=key)
    return None

if __name__ == '__main__':
    client = client_from_stdin()
    generator = MaclaurinAngularMomentum()
    if '--fill-draft' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(client=client, message='Certify squared angular momentum at the first ten axisymmetric Maclaurin bifurcations.'))
    else:
        report = generator.verify(sample=None, client=client)
        print(report)
        sys.exit(0 if report.ok else 1)