"""Independent arithmetic and embedding checks for numberdb.org/T657.

    sage -python checks.py
    sage -python checks.py --stored
For a private draft, --stored requires NUMBERDB_API_KEY. The first command
checks generated values without writing. The second reads all stored values
and applies the same independent tests. Both print a JSON evidence record.

Numerical checking uses real Euler products, not Arb's modular_eta. With
q=exp(-2*pi*sqrt(r)), truncate each Euler product after N=160 factors. Its tail
is in [1-q^(N+1)/(1-q),1]. Widening 1 by that bound is conservative and all
arithmetic is outward-rounded real ball arithmetic at 768 bits.

For B, independently substitute m=C/(x-b) into the canonical polynomial,
clear its denominator and eliminate j using PARI polresultant. For A,
eliminate B from its already checked minimal polynomial and the quadratic
A=B+C/(B-b). Require exact divisibility, irreducibility, vanishing on the
finite-product enclosure, and nonvanishing of the quotient after removing
all copies of the selected factor. Thus no other irreducible factor can
contain the specified value. Controls with known and deliberately wrong
answers precede these tests. Input H polynomials are checked with PARI polclass
and FLINT's certified Hilbert routine; canonical equations with PARI's data.
"""
import json
import os
import sys
from fractions import Fraction
from pathlib import Path

import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.real_arb import RealBallField
from sage.libs.pari import pari
from sage.schemes.elliptic_curves.cm import hilbert_class_polynomial
from sage.env import SAGE_VERSION

import generate as G

BITS = 768
TERMS = 160
RB = RealBallField(BITS)
ZX = PolynomialRing(ZZ, 'x')
x = ZX.gen()
INPUTS = json.loads((Path(__file__).resolve().parent / 'cm-inputs.json').read_text())


def precise(v, minimum=600):
    assert v.is_finite(), 'nonfinite independent enclosure'
    assert v.accuracy() >= minimum, 'independent enclosure too wide'
    return v


def euler(q):
    partial = RB(1)
    power = q
    for n in range(1, TERMS+1):
        partial *= 1-power
        power *= q
    omitted_bound = q**(TERMS+1)/(1-q)
    return partial * RB(1).add_error(omitted_bound)


def independent_values(p, r):
    t = 2*RB.pi()*RB(r).sqrt()
    q = (-t).exp()
    ratio = t.exp()*(euler(q)/euler(q**p))**(24//(p-1))
    b, c = 24//(p-1), ZZ(p)**(12//(p-1))
    return precise(ratio+b+c/ratio), precise(ratio+b)


def certify_factor(annihilator, selected, value):
    precise(value)
    ann, selected = ZX(annihilator).change_ring(QQ), ZX(selected).change_ring(QQ)
    assert selected.is_monic()
    assert bool(pari(selected).polisirreducible())
    quotient, remainder = ann.quo_rem(selected)
    assert remainder == 0, 'candidate does not divide exact annihilator'
    multiplicity = 1
    while quotient.degree() >= selected.degree():
        nxt, rem = quotient.quo_rem(selected)
        if rem != 0:
            break
        quotient = nxt
        multiplicity += 1
    residual, exclusion = selected(value), quotient(value)
    assert residual.is_finite() and exclusion.is_finite()
    assert residual.contains_zero(), 'candidate misses the CM embedding'
    assert not exclusion.contains_zero(), 'another factor may contain the CM embedding'
    return multiplicity


def direct_b_annihilator(p, r):
    # PARI's polynomial variables have priorities; explicitly name elimination.
    b, c = 24//(p-1), p**(12//(p-1))
    h = INPUTS['hilbert'][str(r)]
    phi = INPUTS['canonical'][str(p)]
    expr = ('polresultant(subst(%s,x,j),'
            'subst(subst(%s,y,j),x,%d/(x-%d))*(x-%d)^%d,j)'
            % (h, phi, c, b, b, p+1))
    answer = pari(expr)
    return ZX(str(answer))


def direct_a_annihilator(pb, p):
    b, c = 24//(p-1), p**(12//(p-1))
    expr = ('polresultant(subst(%s,x,t),(t-%d)^2-(x-%d)*(t-%d)+%d,t)'
            % (str(pb).replace('**','^'), b, b, b, c))
    return ZX(str(pari(expr)))


def exact_decimal(text):
    f = Fraction(str(text))
    return QQ(f.numerator)/QQ(f.denominator)


def within_companion(value, text):
    text = str(text)
    if '.' not in text and 'e' not in text.lower():
        assert ZZ(text) in value
        return
    centre = exact_decimal(text)
    mantissa, _, exponent = text.lower().partition('e')
    places = len(mantissa.split('.')[1]) if '.' in mantissa else 0
    unit = QQ(10)**(int(exponent or '0')-places)
    # Positive differences certify containment of the full independent ball.
    assert value-(centre-unit) > 0 and (centre+unit)-value > 0


def controls():
    assert euler(RB(0)) == 1
    a, b = independent_values(2, 1)
    assert 544 in a and 536 in b
    assert 536 not in a and 544 not in b
    certify_factor((x-3)*(x+2), x-3, RB(3))
    rejected = 0
    for ann, poly, ball in [((x-3)*(x+2),x+2,RB(3)),
                           ((x-3)*(x+2),x-4,RB(3)),
                           ((x-3)*(x+2),x-3,RB(3).add_error(10)),
                           ((x-3)*(x+2),x-3,RB('nan'))]:
        try:
            certify_factor(ann,poly,ball)
        except AssertionError:
            rejected += 1
    assert rejected == 4
    return {'known_square_lattice':True,'Euler_empty_product':True,
            'known_linear_factor':True,'deliberately_wrong_or_wide_rejected':rejected}


def read_values(stored):
    if not stored:
        generator = G.MonstrousCMPolynomials()
        return {(d['g'],int(d['r'])):generator.value(d,100)['number'] for d in generator.enumerate()}
    doc = numberdb.table('T657')
    numbers = doc['Numbers']
    if isinstance(numbers,list):
        values = {(str(row['params']['g']),int(row['params']['r'])):ZX(row['number']) for row in numbers}
    else:
        values = {}
        for g, rows in numbers.items():
            for r, entry in rows.items():
                if r in ('comment','param-latex'):
                    continue
                text = entry['number'] if isinstance(entry,dict) else entry
                values[(g,int(r))] = ZX(text)
    expected = [(g,r) for g in G.CLASSES for r in range(1,11)]
    assert list(values) == expected, 'stored enumeration differs'
    return values


def main(stored=False):
    result = {'table':'T657','stored':stored,'controls':controls(),
              'skill_manifest':'agents/runs/20261101T181018Z-build-skills/manifest.json',
              'skill_sha256':'a2b28d9953264d7b0ca62663785c2952190c5be6d6d3597fd695d9a54d606278',
              'sage':SAGE_VERSION,'pari':str(pari('version()')),
              'numberdb':numberdb.__version__,'real_product_bits':BITS,'product_terms':TERMS}
    values = read_values(stored)
    assert len(values) == 100
    # Controls precede every full-range comparison.
    assert ZX(str(pari(-4).polclass())) == x-1728
    assert ZX(hilbert_class_polynomial(-4,algorithm='arb')) == x-1728
    for r in range(1,11):
        wanted = ZX(INPUTS['hilbert'][str(r)])
        assert ZX(str(pari(-4*r).polclass())) == wanted
        assert ZX(hilbert_class_polynomial(-4*r,algorithm='arb')) == wanted
    for p in G.PRIMES:
        canonical = pari('ellmodulareqn(%d)' % p)
        assert canonical[1] == 0
        assert G.QMJ(str(canonical[0])) == G.QMJ(INPUTS['canonical'][str(p)])
    result['hilbert_PARI_and_arb_matches'] = 10
    result['canonical_PARI_matches'] = 5
    details = []
    for p in G.PRIMES:
        for r in range(1,11):
            ba, bb = independent_values(p,r)
            pa, pb = values[(str(p)+'A',r)], values[(str(p)+'B',r)]
            assert pa.base_ring() is ZZ and pb.base_ring() is ZZ
            mb = certify_factor(direct_b_annihilator(p,r),pb,bb)
            ma = certify_factor(direct_a_annihilator(pb,p),pa,ba)
            # Exact A/B relation in the quotient Q[x]/(P_B).
            b, c = 24//(p-1), ZZ(p)**(12//(p-1))
            d = pa.degree()
            composed = sum((pa[i]*((x-b)*x+c)**i*(x-b)**(d-i) for i in range(d+1)),ZX.zero())
            assert composed.quo_rem(pb)[1] == 0
            for suffix,poly,ball,mult in [('A',pa,ba,ma),('B',pb,bb,mb)]:
                g = str(p)+suffix
                within_companion(ball,INPUTS['cm_values'][g][str(r)])
                production = G.quotient_ball(p,r)
                production = production+b+c/production if suffix=='A' else production+b
                assert production.overlaps(ball)
                details.append({'g':g,'r':r,'degree':int(poly.degree()),
                    'characters':len(str(poly)),'factor_multiplicity':mult,
                    'product_accuracy_bits':int(ball.accuracy()),
                    'generator_accuracy_bits':int(production.accuracy()),
                    'embedding_ball':str(ball),
                    'polynomial':str(poly)})
    assert values[('2A',1)] == x-544 and values[('2B',1)] == x-536
    result.update(entries=100,exact_AB_relations=50,certified_factor_checks=100,
                  companion_enclosures=100,minimum_product_accuracy_bits=min(d['product_accuracy_bits'] for d in details),
                  minimum_generator_accuracy_bits=min(d['generator_accuracy_bits'] for d in details),
                  longest=max(details,key=lambda d:d['characters']),details=details)
    print('EVIDENCE '+json.dumps(result,ensure_ascii=False),flush=True)
    return result


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN'):
        numberdb.configure(api_key=sys.stdin.read().strip())
    main(stored='--stored' in sys.argv)