"""Intersection numbers of psi classes -- numberdb.org/T458.

Run with SageMath:
    sage -pip install numberdb
    sage -python generate.py                 # verify every stored entry
    sage -python checks.py                   # independent checks, no writes

For an unpublished draft, verification needs NUMBERDB_API_KEY, or set
NUMBERDB_KEY_FROM_STDIN=1 and pipe the key to the script. This script never
makes a table public. An authorized contributor can fill its entries with
PsiIntersections().publish(); this is distinct from site publication.

The numbers are exact rational stack integrals with labelled marked points.
Profiles are stored once, in non-decreasing order. DVV recursion, string,
and dilaton equations start from <tau_0^3>_0=1 and <tau_1>_1=1/24.
All divisions have a Sage rational operand. Liu-Xu's genus recursion is
implemented separately in checks.py with Python Fraction arithmetic.
"""
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from functools import lru_cache
import os
import sys

MAX_GENUS = 5
MAX_MARKS = 5


def odd_double_factorial(k):
    assert k >= -1 and k % 2 == 1
    answer = 1
    for j in range(1, k + 1, 2):
        answer *= j
    return answer


def intersection(g, degrees):
    return _intersection(int(g), tuple(sorted(int(d) for d in degrees)))


@lru_cache(None)
def _intersection(g, ds):
    n = len(ds)
    if g < 0 or 2 * g - 2 + n <= 0 or any(d < 0 for d in ds):
        return QQ(0)
    if sum(ds) != 3 * g - 3 + n:
        return QQ(0)
    if g == 0 and ds == (0, 0, 0):
        return QQ(1)
    if g == 1 and ds == (1,):
        return QQ(1) / QQ(24)
    if ds[0] == 0:
        rest = ds[1:]
        return sum((intersection(g, rest[:j] + (d - 1,) + rest[j + 1:])
                    for j, d in enumerate(rest) if d > 0), QQ(0))
    if ds[0] == 1:
        return (2 * g - 3 + n) * intersection(g, ds[1:])
    # Distinguish the smallest exponent to keep the double-factorial sum short.
    lead, rest = ds[0], ds[1:]
    answer = QQ(0)
    for j, d in enumerate(rest):
        coefficient = QQ(odd_double_factorial(2 * lead + 2 * d - 1)) / QQ(odd_double_factorial(2 * d - 1))
        answer += coefficient * intersection(g, rest[:j] + (lead + d - 1,) + rest[j + 1:])
    for a in range(lead - 1):
        b = lead - 2 - a
        bracket = intersection(g - 1, (a, b) + rest)
        for mask in range(1 << len(rest)):
            left = tuple(rest[j] for j in range(len(rest)) if mask & (1 << j))
            right = tuple(rest[j] for j in range(len(rest)) if not mask & (1 << j))
            # Dimension determines the genus of the left component.
            genus_numerator = a + sum(left) + 2 - len(left)
            if genus_numerator < 0 or genus_numerator % 3:
                continue
            h = genus_numerator // 3
            if h <= g:
                bracket += intersection(h, (a,) + left) * intersection(g - h, (b,) + right)
        answer += QQ(odd_double_factorial(2 * a + 1) * odd_double_factorial(2 * b + 1)) * bracket / QQ(2)
    return answer / QQ(odd_double_factorial(2 * lead + 1))


def profiles(total, length, minimum=0):
    """All non-decreasing nonnegative profiles of the given length and sum."""
    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


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

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

    def value(self, params, digits):
        g, n = int(params['g']), int(params['n'])
        ds = tuple(map(int, params['d'].split(',')))
        if len(ds) != n or ds != tuple(sorted(ds)) or min(ds) < 0:
            raise ValueError('The exponent profile must be a non-decreasing list of n nonnegative integers.')
        if g < 0 or 2 * g - 2 + n <= 0 or sum(ds) != 3 * g - 3 + n:
            raise ValueError('The parameters must satisfy stability and the dimension constraint.')
        return intersection(g, ds)


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        numberdb.configure(api_key=sys.stdin.read().strip())
    report = PsiIntersections().verify(sample=None)
    print(report)
    sys.exit(0 if report.ok else 1)