"""Independent angular-integral and serialization checks for numberdb.org/T565.

Run: sage -python checks.py
For checks of the stored draft, set NUMBERDB_CHECK_STORED=1 and provide a key.
The Fourier determinant and Pfaffians use no Gamma factors and no generator
helpers. Their derivation and known controls are in independent-checks.md.
"""
import os
import sys
import json
from math import factorial, comb
from decimal import Decimal
from fractions import Fraction
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from numberdb._write import to_text


def binomial(n, k):
    return QQ(comb(n, k)) if 0 <= k <= n else QQ(0)


def determinant(rows):
    a = [list(map(QQ, row)) for row in rows]
    result = QQ(1)
    for i in range(len(a)):
        pivot = next((j for j in range(i, len(a)) if a[j][i]), None)
        if pivot is None:
            return QQ(0)
        if pivot != i:
            a[i], a[pivot] = a[pivot], a[i]
            result = -result
        v = a[i][i]
        result *= v
        for r in range(i + 1, len(a)):
            ratio = a[r][i] / v
            for s in range(i + 1, len(a)):
                a[r][s] -= ratio*a[i][s]
    return result


def pfaffian(rows):
    a = [list(map(QQ, row)) for row in rows]
    assert len(a) % 2 == 0
    result = QQ(1)
    for i in range(0, len(a), 2):
        pivot = next((j for j in range(i + 1, len(a)) if a[i][j]), None)
        if pivot is None:
            return QQ(0)
        if pivot != i + 1:
            a[i+1], a[pivot] = a[pivot], a[i+1]
            for row in a:
                row[i+1], row[pivot] = row[pivot], row[i+1]
            result = -result
        v = a[i][i+1]
        result *= v
        for r in range(i+2, len(a)):
            for s in range(r+1, len(a)):
                a[r][s] -= (a[i][r]*a[i+1][s] - a[i][s]*a[i+1][r])/v
                a[s][r] = -a[r][s]
    return result


def sign_power(integer):
    return QQ(-1 if int(integer) % 2 else 1)


def skew_mode(a, b):
    """Integral sign(y-x) exp(i*a*x+i*b*y) dxdy, divided by i*pi.

    a,b belong to the same integer or half-integer lattice.
    """
    if a and b == -a:
        return -QQ(4)/a
    if a == 0 and b:
        assert b.denominator() == 1
        return -4*sign_power(b)/b
    if b == 0 and a:
        assert a.denominator() == 1
        return 4*sign_power(a)/a
    return QQ(0)


@lru_cache(None)
def angular_integral(beta, n, k):
    """Exact (r,p) from the defining integral, with value r*pi**p."""
    if beta == 2:
        a = [[binomial(2*k, k+i-j) for j in range(n)] for i in range(n)]
        return QQ(factorial(n))*determinant(a), 0
    if beta == 4:
        a = [[(j-i)*binomial(4*k, 2*k+i+j-(2*n-1))
              for j in range(2*n)] for i in range(2*n)]
        return QQ(factorial(n))*pfaffian(a), 0
    assert beta == 1
    modes = [[(QQ(2*j - (n-1) + k - 2*t)/2, binomial(k,t))
              for t in range(k+1)] for j in range(n)]
    size = n + (n % 2)
    a = [[QQ(0) for j in range(size)] for i in range(size)]
    for i in range(n):
        for j in range(i+1,n):
            a[i][j] = sum((ca*cb*skew_mode(fa,fb)
                           for fa,ca in modes[i] for fb,cb in modes[j]), QQ(0))
            a[j][i] = -a[i][j]
    border_pi = 0
    if n % 2:
        border_pi = int(k % 2 == 0)
        for i in range(n):
            value = QQ(0)
            for f,c in modes[i]:
                if border_pi:
                    if f == 0:
                        value += 2*c
                else:
                    assert f.denominator() == 2
                    sine = sign_power((int(2*f)-1)//2)
                    value += c*2*sine/f
            a[i][n], a[n][i] = value, -value
    m = n//2
    phase = sign_power(m) if n % 2 else QQ(1)
    return phase*QQ(factorial(n))/QQ(2**n)*pfaffian(a), m+border_pi-n


def controls():
    assert determinant([[2,3],[5,7]]) == -1
    assert pfaffian([[0,2,3,5],[-2,0,7,11],[-3,-7,0,13],[-5,-11,-13,0]]) == 28
    assert pfaffian([[0,0,2,0],[0,0,0,3],[-2,0,0,0],[0,-3,0,0]]) == -6
    assert angular_integral(1,1,0) == (QQ(1),0)
    assert angular_integral(1,1,1) == (QQ(4),-1)
    assert angular_integral(1,2,0) == (QQ(4),-1)
    assert angular_integral(1,3,0) == (QQ(6),-1)
    assert angular_integral(2,2,0) == (QQ(2),0)
    assert angular_integral(4,2,0) == (QQ(6),0)
    return 9


def decimal_interval(text):
    decimal = Decimal(text)
    center = Fraction(decimal)
    unit = Fraction(10)**decimal.as_tuple().exponent
    return QQ(center-unit), QQ(center+unit)


def contains_written(text, coefficient, power):
    if '.' not in text and 'e' not in text.lower():
        assert power == 0 and QQ(text) == coefficient
        return
    rb = RealBallField(640)
    ball = rb(coefficient)*rb.pi()**power
    assert ball.is_finite() and ball.accuracy() >= 620
    low, high = decimal_interval(text)
    assert low < QQ(ball.lower().exact_rational()) and QQ(ball.upper().exact_rational()) < high, text


def barnes_integer(m):
    result = QQ(1)
    for j in range(1,m-1):
        result *= factorial(j)
    return result


def main():
    from generate import MorrisIntegrals, exact_form
    total_controls = controls()
    generator = MorrisIntegrals()
    entries, precisions = [], []
    for params in generator.enumerate():
        beta,n,k = (int(params[a]) for a in generator.parameters)
        independent = angular_integral(beta,n,k)
        generated = exact_form(beta,n,k)
        assert generated == independent, (params,generated,independent)
        value = generator.value(params,100)
        text = to_text(value,100)
        contains_written(text,*independent)
        if independent[1]:
            assert value.is_finite()
            precisions.append(int(value.accuracy()))
        if beta == 2:
            g=barnes_integer
            assert independent[0]/QQ(factorial(n)) == g(n+1)*g(n+2*k+1)*g(k+1)**2/(g(n+k+1)**2*g(2*k+1))
            if k == 1:
                assert independent == (QQ(factorial(n+1)),0)
        if k == 1:
            assert independent == angular_integral(beta,n+1,0)
        entries.append({'params':params,'number':text,'rational_coefficient':str(independent[0]),'pi_power':independent[1]})
    if os.environ.get('NUMBERDB_CHECK_STORED') == '1':
        document = numberdb.table('T565')
        stored = document['Numbers']
        count=0
        for beta, ns in stored.items():
            for n, ks in ns.items():
                for k, text in ks.items():
                    if isinstance(text,dict): text=text['number']
                    q,p=angular_integral(int(beta),int(n),int(k))
                    contains_written(str(text),q,p)
                    count+=1
        assert count == 150
        report=generator.verify(sample=None)
        assert report.ok
        print('Stored entries:',count,'independent interval inclusions; generator verify:',report)
    result={'table':'T565','entries':len(entries),'exact_rows':sum(e['pi_power']==0 for e in entries),
            'irrational_rows':sum(e['pi_power']!=0 for e in entries),'independent_integral_checks':len(entries),
            'known_controls':total_controls,'working_bits':numberdb.bits(100,losing=96),
            'minimum_retained_bits':min(precisions),'independent_enclosure_bits':640,
            'longest_serialized':max(len(e['number']) for e in entries),
            'longest_at':max(entries,key=lambda e:len(e['number']))['params'],
            'serialized_values_bytes':sum(len(e['number']) for e in entries)}
    print('EVIDENCE_JSON='+json.dumps(result,sort_keys=True))
    print('ENTRIES_JSON='+json.dumps(entries))


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key=sys.stdin.read().strip()
        if '=' in key:key=key.split('=',1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY']=key
    main()