"""Hecke-type Faber polynomials of monstrous McKay–Thompson series -- numberdb.org/T635.

Run with SageMath:
    sage -pip install numberdb
    sage -python generate.py                 # verify every stored entry
    sage -python generate.py --publish       # fill the draft, NUMBERDB_API_KEY set
    sage -python checks.py                   # independent checks, no network

The --publish option sends entries to the existing draft; it does not make
that draft public. The class index uses ATLAS labels, q=exp(2*pi*i*tau),
zero constant coefficient, and P_n(T_g)=q^(-n)+O(q), without division by n.
All arithmetic is exact over ZZ. See argument-selection.md for the range.
The input series are computed from eta products using their logarithmic
coefficients; the polynomials use the Faber recurrence. The separate checks
use direct product multiplication, Laurent cancellation, archived canonical
modular equations and j coefficients, and published polynomial examples.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

R = PolynomialRing(ZZ, 'X')
X = R.gen()
PRIMES = (2, 3, 5, 7, 13)
CLASSES = ('1A',) + tuple(str(p) + suffix for p in PRIMES for suffix in ('A', 'B'))
MAX_DEGREE = 20


def convolution(a, b, limit):
    return [sum((a[k] * b[n-k] for k in range(max(0, n-len(b)+1), min(n+1, len(a)))), ZZ(0))
            for n in range(limit+1)]


def sigma(n, power=1):
    return sum((ZZ(d)**power for d in range(1, n+1) if n % d == 0), ZZ(0))


def from_log_derivative(log_coeffs, limit):
    out = [ZZ(1)]
    for n in range(1, limit+1):
        total = sum((log_coeffs[k] * out[n-k] for k in range(1, n+1)), ZZ(0))
        assert total % n == 0
        out.append(total // n)
    return out


@lru_cache(None)
def coefficients(g, limit):
    """Return c_g(1),...,c_g(limit) as a list with c_g(0)=0."""
    if g not in CLASSES:
        raise ValueError('Unsupported Monster class: ' + g)
    length = limit+1                 # q*T requires one additional coefficient
    if g == '1A':
        e4 = [ZZ(1)] + [240*sigma(n, 3) for n in range(1, length+1)]
        reciprocal_delta = from_log_derivative([ZZ(0)]+[24*sigma(n) for n in range(1, length+1)], length)
        numerator = convolution(convolution(e4, e4, length), e4, length)
        qj = convolution(numerator, reciprocal_delta, length)
        assert qj[:2] == [1, 744]
        return tuple([ZZ(0)] + qj[2:])
    p = int(g[:-1]); e = 24 // (p-1)
    logs = [ZZ(0)] + [-e*(sigma(n) - (p*sigma(n//p) if n % p == 0 else 0))
                      for n in range(1, length+1)]
    product = from_log_derivative(logs, length)
    assert product[:2] == [1, -e]
    c = [ZZ(0)] + product[2:]
    if g[-1] == 'A':
        inverse = from_log_derivative([-v for v in logs], length)
        scale = ZZ(p)**(12 // (p-1))
        for n in range(1, limit+1):
            c[n] += scale*inverse[n-1]
    return tuple(c)


@lru_cache(None)
def faber(g, degree):
    if degree < 1:
        raise ValueError('The degree must be positive')
    c = coefficients(g, degree-1)
    polynomials = [R(1), X]
    for n in range(2, degree+1):
        polynomials.append(X*polynomials[n-1]
                           - sum((c[k]*polynomials[n-1-k] for k in range(1, n-1)), R(0))
                           - n*c[n-1])
    return polynomials[degree]


class MonstrousFaber(numberdb.Generator):
    table = 'T635'
    parameters = ('g', 'n')
    type = 'Z[]'
    rigour = 'exact'
    digits = 100
    files = ('generate.py', 'checks.py', 'check-data.json', 'argument-selection.md',
             'numerical-evidence.md', 'skill-manifest.json')

    def enumerate(self, max_degree=MAX_DEGREE):
        for g in CLASSES:
            for n in range(2, max_degree+1):
                yield {'g': g, 'n': n}

    def value(self, params, digits):
        return faber(str(params['g']), int(params['n']))


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key = sys.stdin.read().strip()
        if '=' in key and key.split('=', 1)[0].strip().isidentifier():
            key = key.split('=', 1)[1].strip().strip('\"\'')
        numberdb.configure(api_key=key)
    generator = MonstrousFaber()
    if '--publish' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Exact monic Faber polynomials for eleven Monster classes, degrees 2 through 20'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)