"""Ellenberg–Gijswijt constants and reduced polynomial multiplication subrank -- numberdb.org/T655.

Run with SageMath and the numberdb client:
    sage -pip install numberdb
    sage -python generate.py              # verify all stored entries
    sage -python generate.py --fill-draft # write numbers; never publishes a table

The write mode requires NUMBERDB_API_KEY and an existing writable draft.
For the campaign wrapper, NUMBERDB_KEY_FROM_STDIN=1 reads the key from stdin;
NUMBERDB_PUBLISH=1 selects the same fill-only operation.

For every integer n >= 2, isolate the unique positive zero of
P_n(t) = sum((3*j-n+1)*t**j, j=0..n-1) using exact rational signs.
Strict convexity of log(S_n(exp(s))) - (n-1)*s/3 identifies its minimum.
Arb evaluates the positive-real objective on the entire root enclosure.
The default enumeration is all n=2..32, with C_n and J_n=C_n/n at each n.

128 extra working bits and 32 additional root-bisection bits are used.
The accompanying checks measure retained precision and test containment of
all 100-digit serialized intervals. See argument-selection.md, checks.py,
and skill-evidence.json for conventions, evidence, and the archived skill.
"""
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
from sage.rings.real_mpfi import RealIntervalField

WORKING_GUARD = 128
ROOT_GUARD = 32
MAX_DIMENSION = 32
CONVENTIONS = ('subrank', 'per-dimension')


def stationary(n, t):
    """P_n(t), evaluated by Horner's rule without floating-point arithmetic."""
    value = QQ(0)
    for j in reversed(range(n)):
        value = value * t + (3*j - n + 1)
    return value


@lru_cache(maxsize=None)
def root_bracket(n, steps):
    """Exact rational enclosure of the global minimizing argument."""
    if n < 2:
        raise ValueError('n must be an integer at least two')
    lo, hi = QQ(0), QQ(1)
    assert stationary(n, lo) < 0 < stationary(n, hi)
    for unused in range(steps):
        middle = (lo + hi) / QQ(2)
        sign = stationary(n, middle)
        if sign == 0:
            return middle, middle
        if sign < 0:
            lo = middle
        else:
            hi = middle
    assert stationary(n, lo) < 0 < stationary(n, hi)
    return lo, hi


@lru_cache(maxsize=None)
def constant(n, digits=100):
    precision = numberdb.bits(digits, losing=WORKING_GUARD)
    lo, hi = root_bracket(n, precision + ROOT_GUARD)
    field = RealBallField(precision)
    rho = field(RealIntervalField(precision)(lo, hi))
    total = field(0)
    for unused in range(n):
        total = total * rho + 1
    value = total * (-field(QQ(n-1) / QQ(3)) * rho.log()).exp()
    assert value.is_finite() and value > 0
    return value


class EllenbergGijswijt(numberdb.Generator):
    table = 'T655'
    parameters = ('n', 'convention')
    type = 'R'
    digits = 100
    rigour = 'proven'
    files = ('generate.py', 'checks.py', 'argument-selection.md',
             'skill-evidence.json', 'numerical-evidence.json', 'verification.md',
             'entropy-reference.json', 'source-manifest.json')

    def enumerate(self, max_dimension=MAX_DIMENSION):
        for n in range(2, max_dimension + 1):
            for convention in CONVENTIONS:
                yield {'n': n, 'convention': convention}

    def value(self, params, digits):
        n = int(params['n'])
        if str(n) != str(params['n']) or n < 2:
            raise ValueError('n must be an integer at least two')
        convention = params['convention']
        if convention not in CONVENTIONS:
            raise ValueError('unknown normalization')
        value = constant(n, digits)
        if convention == 'per-dimension':
            value = value / n
        entry = {'number': value,
                 'param-latex': '$C_n$' if convention == 'subrank' else '$J_n$'}
        if n == 2 and convention == 'subrank':
            entry['comment'] = r'$C_2=3/2^{2/3}$.'
        if n == 3 and convention == 'subrank':
            entry['comment'] = r'$C_3=3(207+33\sqrt{33})^{1/3}/8$, the Ellenberg–Gijswijt cap-set bound constant.'
        return entry


def configure_from_stdin():
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key = sys.stdin.read().strip()
        if '=' in key and key.split('=', 1)[0].isidentifier():
            key = key.split('=', 1)[1].strip().strip('\"\'')
        numberdb.configure(api_key=key)


if __name__ == '__main__':
    configure_from_stdin()
    generator = EllenbergGijswijt()
    if '--fill-draft' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Fill 62 certified constants in both source normalizations'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)