"""Lambda_g Hodge integrals -- numberdb.org/T462.

Run with SageMath (or a suitable passagemath installation):
    sage -pip install numberdb
    sage -python generate.py              # verify every stored entry
    sage -python checks.py                # independent checks, no API key
    sage -python generate.py --fill       # write entries to the existing draft

The last command requires NUMBERDB_API_KEY for the draft's owner. It writes
numbers and attachments; it does not publish or offer the draft for review.

The table uses the Deligne-Mumford stack with labelled markings. Profiles are
sorted representatives, without division by a permutation multiplicity.
Their total degree is 2*g-3+n. The first five positive genera and at most five
markings form a complete low-genus reference rectangle. Genus zero, where
lambda_0=1, is described in prose and already covered by the psi-class table.

All arithmetic is exact in QQ. The coefficient of t^(2*g) in
(t/2)/sin(t/2) is obtained by triangular coefficient inversion, then multiplied
by a multinomial coefficient. No approximate numerical arithmetic is used.
"""
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from functools import lru_cache
from math import comb, factorial
import os
import sys

MAX_GENUS = 5
MAX_MARKS = 5


def profiles(total, length, minimum=0):
    """All non-decreasing nonnegative profiles, including zero exponents."""
    if length == 0:
        if total == 0:
            yield ()
        return
    for first in range(minimum, total // length + 1):
        for rest in profiles(total-first, length-1, first):
            yield (first,) + rest


@lru_cache(None)
def one_point(g):
    """Coefficient of t^(2*g) in (t/2)/sin(t/2), including b_0=1."""
    if g == 0:
        return QQ(1)
    return -sum((QQ((-1)**k) / QQ(2**(2*k)*factorial(2*k+1)))
                * one_point(g-k) for k in range(1, g+1))


def integral(g, ds):
    """A stable top-degree lambda_g integral; ds need not be sorted."""
    g = int(g)
    ds = tuple(int(d) for d in ds)
    n = len(ds)
    if g < 0 or any(d < 0 for d in ds) or 2*g-2+n <= 0:
        raise ValueError('A stable pair and nonnegative exponents are required')
    if sum(ds) != 2*g-3+n:
        raise ValueError('The exponent profile must have degree 2*g-3+n')
    coefficient, remaining = 1, sum(ds)
    for d in ds:
        coefficient *= comb(remaining, d)
        remaining -= d
    return one_point(g) * coefficient


class LambdaGHodge(numberdb.Generator):
    table = 'T462'
    parameters = ('g', 'n', 'd')
    type = 'Q'
    rigour = 'exact'
    files = ('generate.py', 'checks.py', 'check.py', 'evidence.json',
             'bernoulli.json', 'T458.json')

    def enumerate(self, max_genus=MAX_GENUS, max_marks=MAX_MARKS):
        for g in range(1, max_genus+1):
            for n in range(1, max_marks+1):
                assert 2*g-2+n > 0
                for ds in profiles(2*g-3+n, n):
                    assert len(ds) == n and sum(ds) == 2*g-3+n
                    yield {'g': g, 'n': n, 'd': ','.join(map(str, ds))}

    def value(self, params, digits=100):
        g, n = int(params['g']), int(params['n'])
        ds = tuple(map(int, str(params['d']).split(',')))
        if len(ds) != n or tuple(sorted(ds)) != ds:
            raise ValueError('The profile must be sorted and have n entries')
        result = integral(g, ds)
        if (g, ds) == (1, (0,)):
            return {'number': result,
                    'comment': r'The integral $\int_{\overline{\mathcal M}_{1,1}}\lambda_1=1/24$ includes the generic elliptic involution.'}
        return result


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
    generator = LambdaGHodge()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Fill exact lambda_g integrals after independent checks'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)