"""Independent checks for numberdb.org/T571; run: sage -python checks.py.

Use --stored to repeat the checks on values read from NumberDB (a private draft
requires NUMBERDB_API_KEY). No mutation is performed. This checker evaluates
Laguerre integrals by exact moments: a determinant for beta=2, or a de Bruijn
Pfaffian for beta=1. It does not use the generator's gamma product.

For p,q >= 0, I(p,q)=integral_{0<x<y} x^p y^q exp(-x-y) dx dy
= q! sum_{k=0}^q (p+k)!/(k! 2^(p+k+1)). The skew moment is I(p,q)-I(q,p).
Odd ranks are bordered by the single moments p!. See de Bruijn (1955),
On some multiple integrals involving determinants, (1.1), (4.3)-(4.5).
The trace constraint divides the Laguerre integral by (A-1)!, where
A=r*(1+beta*(N-r))+beta*r*(r-1)/2. The orbit volume is independently
assembled from successive real or complex projective-space volumes.

Controls run first: exact determinant/Pfaffian examples and eigenvalue
integrals; all four N=2 volumes come from a disk/ball of radius 1/sqrt(2).
The checks also cover pure-state geometry and the boundary/interior ratio.
"""
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from math import factorial
from functools import lru_cache
from decimal import Decimal
from fractions import Fraction
from pathlib import Path
import json
import os
import sys
import time
from check import exactness, measure, names_its_rings
from numberdb._write import to_text
from generate import HilbertSchmidtVolumes


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


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


@lru_cache(None)
def ordered_moment(p,q):
    return QQ(factorial(q))*sum((QQ(factorial(p+k))/QQ(factorial(k)*2**(p+k+1))
                                for k in range(q+1)),QQ(0))


@lru_cache(None)
def laguerre_integral(beta,d,r):
    if beta==2:
        a=[[factorial(2*d+i+j) for j in range(r)] for i in range(r)]
        return QQ(factorial(r))*determinant(a)
    size=r+(r%2)
    a=[[QQ(0) for j in range(size)] for i in range(size)]
    for i in range(r):
        for j in range(i+1,r):
            a[i][j]=ordered_moment(d+i,d+j)-ordered_moment(d+j,d+i)
            a[j][i]=-a[i][j]
        if r%2:a[i][r]=QQ(factorial(d+i));a[r][i]=-a[i][r]
    return QQ(factorial(r))*pfaffian(a)


def independent_volume(field,N,r,bits=560):
    R=RealBallField(bits);beta=1 if field=='real' else 2;d=N-r
    A=r*(1+beta*d)+beta*r*(r-1)//2
    integral=laguerre_integral(beta,d,r)/QQ(factorial(A-1))
    orbit=R(1)
    for k in range(d+1,N+1):
        if beta==2:orbit*=(2*R.pi())**(k-1)/R(factorial(k-1))
        else:orbit*=R(2)**(QQ(k-1)/2)*R.pi()**(QQ(k)/2)/R(QQ(k)/2).gamma()
    return R(r).sqrt()*R(integral)*orbit/R(factorial(r))


def endpoints(text):
    # Decimal text asserts a radius of one unit in its last place.
    x=Decimal(str(text).rstrip('?'));f=Fraction(x)
    c=QQ(f.numerator)/QQ(f.denominator)
    exponent=x.as_tuple().exponent
    ulp=QQ(10)**exponent
    return c-ulp,c+ulp


def narrow(value,bits=330):
    assert value.is_finite() and value>0
    assert value.accuracy()>=bits,(value,value.accuracy())


def contains_ball(text,ball):
    lo,hi=endpoints(text)
    assert ball.lower()>=lo and ball.upper()<=hi,(text,str(ball))


def run(stored=False):
    start=time.monotonic()
    # Nontrivial, positive, negative, singular and swapped-pivot controls.
    assert determinant([[1,2],[3,4]])==-2
    assert determinant([[0,1],[2,3]])==-2
    assert determinant([[1,2],[2,4]])==0
    a=[[0,2,3,5],[-2,0,7,11],[-3,-7,0,13],[-5,-11,-13,0]]
    assert pfaffian(a)==28 and determinant(a)==28**2
    assert pfaffian([[0,0,1,0],[0,0,0,1],[-1,0,0,0],[0,-1,0,0]])==-1
    assert laguerre_integral(1,0,1)==1
    assert laguerre_integral(1,0,2)==1
    assert laguerre_integral(2,0,2)==2
    print('PASS controls: determinant, Pfaffian, and small eigenvalue integrals')
    g=HilbertSchmidtVolumes();params=list(g.enumerate())
    assert len(params)==270 and len({tuple(p.values()) for p in params})==270
    assert not names_its_rings(str(Path(__file__).with_name('generate.py')))
    values={tuple(p.values()):g.value(p,g.digits) for p in params}
    assert not exactness(values)
    if stored:
        doc=numberdb.table(g.table)
        rows=doc['Numbers']
        assert isinstance(rows,list) and len(rows)==270
        texts={(str(x['params']['field']),int(x['params']['N']),int(x['params']['rank'])):x['number'] for x in rows}
        assert list(texts)==list(values)
    else:texts={key:to_text(value,digits=100) for key,value in values.items()}
    accuracy=[]
    for key,value in values.items():
        narrow(value);accuracy.append(int(value.accuracy()))
        oracle=independent_volume(*key)
        narrow(oracle,480)
        assert value.overlaps(oracle),key
        contains_ball(texts[key],value)
        contains_ball(texts[key],oracle)
    print('PASS all 270 exact-moment computations and serialized-interval containments')
    R=RealBallField(560)
    controls={('real',2,1):R(2).sqrt()*R.pi(),('real',2,2):R.pi()/2,
              ('complex',2,1):2*R.pi(),('complex',2,2):R.pi()*R(2).sqrt()/3}
    for key,v in controls.items():contains_ball(texts[key],v)
    for field in ('real','complex'):
        for N in range(2,17):
            if field=='complex':pure=(2*R.pi())**(N-1)/R(factorial(N-1))
            else:pure=R(2)**(QQ(N-1)/2)*R.pi()**(QQ(N)/2)/R(QQ(N)/2).gamma()
            contains_ball(texts[(field,N,1)],pure)
            D=N*N-1 if field=='complex' else N*(N+1)//2-1
            ratio=R(D)*R(N*(N-1)).sqrt()
            # Interval arithmetic applied to database decimal enclosures.
            plo,phi=endpoints(texts[(field,N,N-1)])
            vlo,vhi=endpoints(texts[(field,N,N)])
            assert plo/vhi<=ratio.lower() and ratio.upper()<=phi/vlo
            if field=='complex':
                interior=R(N).sqrt()*(2*R.pi())**(N*(N-1)//2)
                interior*=R(QQ(__import__('functools').reduce(lambda x,k:x*factorial(k),range(N),1))/QQ(factorial(N*N-1)))
                contains_ball(texts[(field,N,N)],interior)
    print('PASS four Bloch disk/ball controls, 30 pure-state volumes, 30 boundary ratios, 15 complex full-rank formulas')
    report={'table':g.table,'stored_readback':stored,'entries':len(values),'min_accuracy_bits':min(accuracy),
            'working_bits':numberdb.bits(100,losing=64),'independent_bits':560,
            'serialized_max_characters':max(map(len,texts.values())),
            'measure':measure(values),'seconds':round(time.monotonic()-start,3),
            'skill_manifest':'agents/runs/20261028T042605Z-build-skills/manifest.json',
            'skill_sha256':'e37efb7e2ededf6c1d7134c673f152f648f26e0267c556aebdf3ec65765d8b4c'}
    print('EVIDENCE '+json.dumps(report,sort_keys=True))
    for key in controls:print('CONTROL',key,texts[key])
    return report


if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
    run('--stored' in sys.argv)