"""Spectral determinants of the Laplacian on real projective spaces -- numberdb.org/T600.

Run with SageMath (install numberdb once with sage -pip install numberdb):
    sage -python generate.py             # verify every stored entry
    sage -python generate.py --fill      # fill this private draft, API key required

The radius-one quotient metric, untwisted scalar Laplacian, and omission of
its constant zero mode are fixed. This evaluates Cunha--Freitas,
arXiv:2404.12114v1, Theorem 3.3. Their s(n,k) are SIGNED falling-factorial
coefficients; v(n,i) come from x*product(x-(2j-1)^2, j=1..n-1).
Dimension one is evaluated directly: L_1=2 log(pi), D_1=pi^2.
All coefficients are exact; transcendental operations use Arb balls.
Working precision is 150 decimal digits for the 100-digit output, with
comparison at 200 digits recorded in numerical-evidence.json.
The finite dimension cutoff is explained in argument-selection.md.
No action in this file offers, reviews, or makes a table public.
"""
import os
import sys
from functools import lru_cache
from math import factorial
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_arb import RealBallField

WORKING_DECIMAL_GUARD = 50

@lru_cache(None)
def signed_stirling(n):
    row = [ZZ(1)]
    for j in range(n):
        nxt = [ZZ(0)] * (len(row)+1)
        for i, c in enumerate(row):
            nxt[i] -= j*c
            nxt[i+1] += c
        row = nxt
    return tuple(row)

@lru_cache(None)
def central_odd(n):
    row = [ZZ(0), ZZ(1)]
    for j in range(1,n):
        nxt = [ZZ(0)] * (len(row)+1)
        for i,c in enumerate(row):
            nxt[i] -= (2*j-1)**2*c
            nxt[i+1] += c
        row = nxt
    return tuple(row)

def harmonic(n):
    return sum((QQ(1)/j for j in range(1,n+1)),QQ(0))

@lru_cache(None)
def zeta_derivative(k,working_digits):
    C = ComplexBallField(numberdb.bits(working_digits))
    z = C(-k).zetaderiv(1)
    assert z.real().is_finite() and z.imag().is_finite() and z.imag().contains_zero()
    return z.real()

@lru_cache(None)
def logdet(d,working_digits=150):
    d=int(d)
    if d<1:
        raise ValueError('Dimension must be a positive integer')
    R=RealBallField(numberdb.bits(working_digits))
    if d==1:
        return 2*R.pi().log()
    s1=signed_stirling(d); s0=signed_stirling(d-1)
    if d%2:
        n=(d+1)//2
        zp=sum((R(QQ(2**(2*k+1)*(s1[2*k+1]+s0[2*k]))/factorial(d-1))
                *zeta_derivative(2*k,working_digits) for k in range(n)),R(0))
        zp+=R(4*(n-1)).log()
    else:
        n=d//2
        zp=sum((R(QQ((1-(k%2==0)*2**(k+1))*(s1[k+1]+s0[k]))/factorial(d-1))
                *zeta_derivative(k,working_digits) for k in range(d)),R(0))
        zp+=R(4*n-2).log()
        v=central_odd(n)
        anomaly=sum((QQ(v[i])/ (4**(n-1)*factorial(2*n-1))
                     *QQ((2*n-1)**(2*i))/(4*i)
                     *(harmonic(i-1)/2-harmonic(2*i-1))
                     for i in range(1,n+1)),QQ(0))
        zp+=R(anomaly)
    answer=-zp
    assert answer.is_finite()
    return answer

class ProjectiveDeterminants(numberdb.Generator):
    table='T600'
    parameters=('d','convention')
    type='R'
    digits=100
    rigour='proven'
    files=('generate.py','check.py','checks.py','argument-selection.md','continuation.md',
           'numerical-evidence.json','skill-evidence.json')

    def enumerate(self):
        for d in range(1,21):
            for convention in ('det','logdet'):
                yield {'d':str(d),'convention':convention}

    def value(self,params,digits):
        L=logdet(int(params['d']),working_digits=digits+WORKING_DECIMAL_GUARD)
        value=L.exp() if params['convention']=='det' else L
        result={'number':value}
        if params['d']=='1':
            result['comment']=(r'$D_1=\pi^2$.' if params['convention']=='det'
                               else r'$L_1=2\log\pi$.')
        return result

if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        raw=sys.stdin.read().strip()
        if '=' in raw and raw.split('=',1)[0].strip().replace('_','').isupper():
            raw=raw.split('=',1)[1].strip().strip('\"\'')
        os.environ['NUMBERDB_API_KEY']=raw
    g=ProjectiveDeterminants()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
        print(g.publish(message='Compute 40 projective-space determinant values with certified Arb enclosures'))
    else:
        report=g.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)