"""Unit roots of hypergeometric Calabi–Yau motives -- numberdb.org/T675.

Run with SageMath, with the attached frobenius-input.json beside this file:
    sage -pip install numberdb
    sage -python generate.py                 # verify all stored entries
    sage -python checks.py                   # independent offline checks
    sage -python generate.py --fill          # fill the draft, with an API key

This never offers, reviews or publishes a table. The exact polynomial input is
an archived API read of T674. checks.py recomputes all 390 factors with PARI
and compares all selected roots with a power-sum ratio. Extending the prime
range requires a newly checked polynomial input and argument-selection review.
See argument-selection.md for all three parameter choices.

Newton iteration is performed with integers modulo successively squared prime
powers. Q(r)=0 mod p^N and Q'(r) !=0 mod p give a Hensel certificate for the
returned p-adic enclosure. N is the least integer with p^N >=10^digits.
The print mode uses one integer residue plus O(p^N), avoiding long digit sums.
"""
import json
import os
import sys
from functools import lru_cache
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.padics.factory import Qp

TEMPLATES = ('2.2.2.2','2.2.3','2.2.4','2.2.6','3.3','3.4','3.6',
             '4.4','4.6','6.6','5','8','10','12')
FIBRES = ('1/2','-1','2')
PRIME_BOUND = 31
PARENT_SLUG = 'Frobenius_polynomials_of_hypergeometric_Calabi_Yau_motives'
R = PolynomialRing(ZZ, 'T')

@lru_cache(None)
def factors():
    data = json.loads(Path(__file__).with_name('frobenius-input.json').read_text())
    assert data['source_table'] == 'T674'
    result = {}
    for template, fibres in data['Numbers'].items():
        for z, primes in fibres.items():
            for p, text in primes.items():
                if isinstance(text, dict): text = text['number']
                poly = R(text)
                assert poly.degree() == 4 and poly[0] == 1
                result[(template,z,int(p))] = tuple(int(c) for c in poly.list())
    expected = {(key,z,p) for key in TEMPLATES for z in FIBRES
                for p in range(3,PRIME_BOUND+1)
                if ZZ(p).is_prime(proof=True) and all(int(d)%p for d in key.split('.'))}
    assert set(result) == expected
    return result

def exponent(p, digits):
    n, modulus = 0, 1
    while modulus < 10**int(digits):
        n += 1
        modulus *= p
    return n

def evaluate(coefficients, x):
    value = 0
    for c in reversed(coefficients):
        value = value*x+c
    return value

def hensel(coefficients, p, n, seed):
    derivative = [i*c for i,c in enumerate(coefficients)][1:]
    r, power = int(seed)%p, 1
    assert evaluate(coefficients,r)%p == 0
    assert evaluate(derivative,r)%p != 0
    while power < n:
        power = min(2*power,n)
        modulus = p**power
        r = (r-evaluate(coefficients,r)*pow(evaluate(derivative,r),-1,modulus))%modulus
        assert evaluate(coefficients,r)%modulus == 0
    assert r%p != 0
    return r

class CalabiYauUnitRoots(numberdb.Generator):
    table = 'T675'
    parameters = ('template','z','p')
    type = 'Qp'
    rigour = 'proven'
    digits = 100
    files = ('generate.py','checks.py','frobenius-input.json',
             'argument-selection.md','numerical-evidence.md','skill-evidence.json','check.py')

    def enumerate(self):
        data = factors()
        for template in TEMPLATES:
            for z in FIBRES:
                for p in range(3,PRIME_BOUND+1):
                    coeff = data.get((template,z,p))
                    if coeff is not None:
                        assert coeff[2]%p == 0
                        if coeff[1]%p:
                            yield {'template':template,'z':z,'p':str(p)}

    def value(self, params, digits=100):
        key = (params['template'],str(params['z']),int(params['p']))
        coeff = factors()[key]
        p, a = key[2], -coeff[1]
        assert a%p and coeff == (1,-a,coeff[2],-p**3*a,p**6)
        n = exponent(p,digits)
        r = hensel(tuple(reversed(coeff)),p,n,a)
        field = Qp(p,prec=n,type='capped-rel',print_mode='terse')
        value = field(r,absprec=n)
        assert value.precision_absolute() == n and value.valuation() == 0
        return {'number':value,'comment':
                'The corresponding local factor is HREF{%s#%s,%s,%s}[$P_{A,z,p}(T)$].'
                % (PARENT_SLUG, key[0],key[1],key[2])}

if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
    generator = CalabiYauUnitRoots()
    if '--fill' in sys.argv:
        print(generator.publish(message='Fill certified Calabi-Yau unit eigenvalues after independent checks'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)