"""Minkowski-Siegel masses of unimodular genera -- numberdb.org/T497.

Run with SageMath (or a compatible numberdb.sage installation):

    sage -pip install numberdb
    sage -python generate.py                 # verify every stored entry
    sage -python checks.py                   # independent checks, offline
    sage -python checks.py --stored          # check database values
    sage -python generate.py --publish       # send entries, API key required

The last command fills the existing table; it does not publish the draft.
For a private draft, reads also require NUMBERDB_API_KEY. In the campaign,
NUMBERDB_KEY_FROM_STDIN=1 reads that key from stdin without writing it.

The exact formulas are Theorem 2 of Conway and Sloane, European Journal of
Combinatorics 3 (1982), 219-231, DOI 10.1016/S0195-6698(82)80034-6.
They use the full integral isometry group, including determinant -1 maps.
Type I enumerates odd genera, type II even genera. Dimension zero is the
single even lattice with trivial isometry group, so its mass is 1.

The range follows their reference Tables 2 and 4: all odd dimensions
1..32, and even dimensions 0,8,16,24,32. These include the E8, dimension-16
and Niemeier classifications and the growth of the masses at dimension 32.
The rank-one printed value 1 in Table 2 is replaced by 1/2, as required
by the full isometry group {+1,-1}; King (2003), Table 6, agrees.
Bernoulli and Euler numbers are computed by exact rational recurrences;
no floating-point arithmetic, numerical L-functions or rounded data enter.
"""

import os
import sys
from functools import lru_cache
from math import comb, factorial

import numberdb.sage as numberdb
from sage.rings.rational_field import QQ

LIMIT = 32


@lru_cache(None)
def bernoulli(n):
    """Coefficient convention t/(exp(t)-1) = sum B_n t^n/n!."""
    if n == 0:
        return QQ(1)
    return -sum((QQ(comb(n + 1, j)) * bernoulli(j)
                 for j in range(n)), QQ(0)) / QQ(n + 1)


@lru_cache(None)
def euler(n):
    """Signed Euler numbers, 1/cosh(t) = sum E_n t^n/n!."""
    if n == 0:
        return QQ(1)
    if n % 2:
        return QQ(0)
    return -sum((QQ(comb(n, j)) * euler(j)
                 for j in range(0, n, 2)), QQ(0))


def bernoulli_product(k):
    result = QQ(1)
    for j in range(1, k + 1):
        result *= bernoulli(2 * j)
    return abs(result)


def mass(parity, n):
    """Exact mass for a positive definite integral unimodular genus."""
    n = int(n)
    if parity == 'II':
        if n < 0 or n % 8:
            raise ValueError('an even unimodular genus requires n >= 0, 8 | n')
        if n == 0:
            return QQ(1)
        k = n // 2
        answer = abs(bernoulli(k)) / QQ(2 * k)
        for j in range(1, k):
            answer *= abs(bernoulli(2 * j)) / QQ(4 * j)
        return answer
    if parity != 'I' or n < 1:
        raise ValueError('an odd unimodular genus requires n >= 1')
    if n == 1:
        return QQ(1) / QQ(2)
    k = n // 2
    if n % 2:
        sign = 1 if n % 8 in (1, 7) else -1
        return (QQ(2**k + sign) * bernoulli_product(k)
                / QQ(factorial(k) * 2**(2*k + 1)))
    if n % 4 == 2:
        return (abs(euler(k - 1)) * bernoulli_product(k - 1)
                / QQ(factorial(k - 1) * 2**(2*k + 1)))
    sign = 1 if n % 8 == 0 else -1
    return ((QQ(1) - QQ(1) / QQ(2**k))
            * (QQ(1) + sign * QQ(1) / QQ(2**(k - 1)))
            * abs(bernoulli(k)) * bernoulli_product(k - 1)
            / QQ(2 * factorial(k)))


class UnimodularMasses(numberdb.Generator):
    table = 'T497'
    parameters = ('parity', 'n')
    type = 'Q'
    rigour = 'exact'
    files = ('generate.py', 'checks.py', 'source_values.json', 'check.py',
             'evidence.txt', 'niemeier_orders.json')

    def enumerate(self, limit=LIMIT):
        for n in range(1, limit + 1):
            yield {'parity': 'I', 'n': str(n)}
        for n in range(0, limit + 1, 8):
            yield {'parity': 'II', 'n': str(n)}

    def value(self, params, digits=None):
        parity, n = params['parity'], int(params['n'])
        comments = {
            ('I', 8): r'The only odd class in this dimension is $\mathbb{Z}^8$.',
            ('II', 0): 'The zero-dimensional lattice has trivial isometry group.',
            ('II', 8): r'The genus consists of the single isometry class $E_8$.',
            ('II', 16): r'The two classes are $E_8\oplus E_8$ and $D_{16}^{+}$.',
            ('II', 24): 'This is the sum of reciprocal full automorphism orders of the 24 Niemeier lattices.',
        }
        answer = mass(parity, n)
        if (parity, n) in comments:
            return {'number': answer, 'comment': comments[parity, n]}
        return answer


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
    gen = UnimodularMasses()
    if '--publish' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(gen.publish(message='Fill exact genus masses through dimension 32, independently checked.'))
    else:
        report = gen.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)