"""Growth constants of the classes of planar graphs -- numberdb.org/T481.

Run with SageMath or passagemath and the numberdb package:
    sage -pip install numberdb
    sage -python generate.py                 # verify every stored entry
    sage -python checks.py                   # independent local checks
    sage -python generate.py --fill          # fill the existing draft

A private draft requires NUMBERDB_API_KEY even for reads. Alternatively set
NUMBERDB_KEY_FROM_STDIN=1 and pipe the key on stdin. NUMBERDB_PUBLISH=1 has
same effect as --fill in environments that do not forward script arguments.
Filling entries does not publish a draft or offer it for review.

All graphs are simple and labelled; size counts vertices. The implicit root
is the unique t in (0,1) with Y(t)=1, in Giménez--Noy, JAMS 22 (2009),
Theorem 1 and Lemma 4. Rational endpoints bracket it, with every sign
certified in Arb. Analytic expressions are evaluated on the enclosing ball.
The 3-connected constant is (272+112*sqrt(7))/27; the maximal constant is
exactly 256/27. The latter follows from Tutte's triangulation enumeration.

Working precision is bits(digits)+96. Bisection uses bits(digits)+40 steps,
leaving guard bits for cancellation and interval dependency. The 100-digit
run and the independent checks measure the resulting errors, rather than
assuming the working precision survives. No fixed-precision approximation
is wrapped in an interval.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField

WORKING_GUARD = 96
BISECTION_GUARD = 40
CLASSES = ('all', 'connected', '2-connected', '3-connected', 'maximal')


def y_of_t(t):
    exponent = -t**2 * (1-t) * (18+36*t+5*t**2)
    exponent /= 2*(3+t)*(1+2*t)*(1+3*t)**2
    return (1+2*t)/((1+3*t)*(1-t)) * exponent.exp() - 1


def root_bracket(digits):
    field = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
    lo, hi = QQ(3)/5, QQ(2)/3
    assert y_of_t(field(lo)) < 1 < y_of_t(field(hi))
    for _ in range(numberdb.bits(digits, losing=BISECTION_GUARD)):
        mid = (lo+hi)/2
        val = y_of_t(field(mid)) - 1
        if val < 0:
            lo = mid
        elif val > 0:
            hi = mid
        else:
            raise ArithmeticError('Working precision cannot certify root sign')
    assert y_of_t(field(lo)) < 1 < y_of_t(field(hi))
    return field, lo, hi


def radii(t):
    two_connected = (1+3*t)*(1-t)**3/(16*t**3)
    a = (1+t).log()*(3*t-1)*(1+t)**3/(16*t**3)
    a += (1+2*t).log()*(1+3*t)*(1-t)**3/(32*t**3)
    a += (1-t)*(185*t**4+698*t**3-217*t**2-160*t+6)/(64*t*(1+3*t)**2*(3+t))
    connected = (1+3*t).sqrt()*(1-t)**3/(16*t**3)*a.exp()
    return two_connected, connected


@lru_cache(maxsize=8)
def constants(digits):
    field, lo, hi = root_bracket(digits)
    t = field(lo).union(field(hi))
    radius_two, radius_all = radii(t)
    return {'all': 1/radius_all, 'connected': 1/radius_all,
            '2-connected': 1/radius_two,
            '3-connected': (272+112*field(7).sqrt())/27,
            'maximal': QQ(256)/27}


class PlanarGrowth(numberdb.Generator):
    table = 'T481'
    parameters = ('class', 'size', 'labelling')
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'check.py', 'convention.md', 'evidence.md')

    def enumerate(self):
        for graph_class in CLASSES:
            yield {'class': graph_class, 'size': 'vertices', 'labelling': 'labelled'}

    def value(self, params, digits):
        if params['size'] != 'vertices' or params['labelling'] != 'labelled':
            raise ValueError('Only labelled graphs counted by vertices are implemented')
        result = {'number': constants(digits)[params['class']]}
        comments = {
            'all': 'Disconnected graphs are included. CITE{GN}',
            'connected': 'The growth constant equals that for all planar graphs. CITE{GN}',
            '2-connected': 'Connectivity is vertex connectivity. CITE{GN}',
            '3-connected': r'The exact constant is $(272+112\sqrt{7})/27$. CITE{Survey}',
            'maximal': r'The exact constant is $256/27$, from the enumeration of simple triangulations. CITE{GN} CITE{Tutte}'
        }
        result['comment'] = comments[params['class']]
        return result


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
    gen = PlanarGrowth()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(gen.publish(message='Fill five labelled planar graph growth constants with certified enclosures', assisted_by='codex-cli'))
    else:
        report = gen.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)