"""Exact finite determinant formulas used by numberdb.org/T604.

Cunha and Freitas, arXiv:2404.12114v1, Theorems 2.7, 2.8 and 3.2.
Adapted from the attached generators of NumberDB T595 and T603, downloaded
in this run and checked against the cited formulas. These compute afresh;
no rounded database values enter the computation.
"""
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_DIGITS = 150

@lru_cache(None)
def stirling_row(n):
    """Signed first-kind Stirling numbers from an integer recurrence."""
    if n == 0:
        return (ZZ(1),)
    previous = stirling_row(n - 1)
    return tuple((previous[k-1] if k else ZZ(0))
                 - (n-1)*(previous[k] if k < n else ZZ(0))
                 for k in range(n+1))


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


@lru_cache(None)
def even_correction(n):
    # Coefficients of x product_{j=1}^{n-1}(x-(2j-1)^2).
    v = [ZZ(0), ZZ(1)]
    for j in range(1, n):
        new = [ZZ(0)]*(len(v)+1)
        for k, coefficient in enumerate(v):
            new[k] -= (2*j-1)**2*coefficient
            new[k+1] += coefficient
        v = new
    scale = QQ(1)/(ZZ(4)**(n-1)*factorial(2*n-1))
    return sum((scale*v[i]*ZZ(2*n-1)**(2*i)/QQ(2*i)
                * (harmonic(i-1)/2-harmonic(2*i-1))
                for i in range(1, n+1)), QQ(0))


@lru_cache(None)
def sphere_logdet(d, working_digits=150):
    if d < 1:
        raise ValueError('The dimension must be a positive integer.')
    C = ComplexBallField(numberdb.bits(working_digits))
    R = RealBallField(numberdb.bits(working_digits))
    if d == 1:
        return 2*(2*R.pi()).log()
    n = d//2
    if d % 2:
        a, b = stirling_row(2*n), stirling_row(2*n+1)
        derivative = sum((C(-2*k).zetaderiv(1)*(a[2*k]+b[2*k+1])
                          for k in range(1, n+1)), C(0))
        derivative *= QQ(2)/factorial(2*n)
        derivative += (R(n)/R.pi()).log()
    else:
        a, b = stirling_row(2*n-1), stirling_row(2*n)
        derivative = sum((C(1-2*k).zetaderiv(1)*(a[2*k-1]+b[2*k])
                          for k in range(1, n+1)), C(0))
        derivative *= QQ(2)/factorial(2*n-1)
        derivative += R(even_correction(n)) + R(2*n-1).log()
    assert derivative.real().is_finite() and derivative.imag().is_finite() and derivative.imag().contains_zero()
    return -derivative.real()


@lru_cache(None)
def stirling(n,k):
    if n == 0:
        return int(k == 0)
    if k < 1 or k > n:
        return 0
    return stirling(n-1,k-1)-(n-1)*stirling(n-1,k)


def roots_polynomial(roots):
    """Exact coefficient list of product (X-r); no symbolic ring needed."""
    c=[1]
    for r in roots:
        out=[0]*(len(c)+1)
        for k,v in enumerate(c):
            out[k]-=r*v
            out[k+1]+=v
        c=out
    return c


@lru_cache(None)
def dirichlet_logdet(d, working=WORKING_DIGITS):
    if not isinstance(d,int) or d<1:
        raise ValueError('dimension must be a positive integer')
    C=ComplexBallField(numberdb.bits(working))
    R=RealBallField(numberdb.bits(working))
    if d == 1:
        return (2*R.pi()).log()
    # The signed Stirling combination is valid for either parity.
    zp=sum((QQ(stirling(d,k+1)+(-1)**(d-1+k)*stirling(d-1,k))
            /QQ(factorial(d-1))*C(-k).zetaderiv(1)
            for k in range(d)), C(0))
    if d%2:
        n=(d+1)//2
        u=roots_polynomial([j*j for j in range(n-1)])
        correction=-sum((QQ(2*u[i])/QQ(factorial(2*n-2))
                         *QQ((n-1)**(2*i+1))/QQ(2*i)
                         *(harmonic(i-1)/2-harmonic(2*i-1))
                         for i in range(1,n)),QQ(0))
    else:
        n=d//2
        v=roots_polynomial([0]+[(2*j-1)**2 for j in range(1,n)])
        correction=sum((QQ(v[i])/QQ(4**(n-1)*factorial(2*n-1))
                        *QQ((2*n-1)**(2*i))/QQ(4*i)
                        *(harmonic(i-1)/2-harmonic(2*i-1))
                        for i in range(1,n+1)),QQ(0))
    result=-(zp+correction)
    assert result.real().is_finite() and result.imag().is_finite()
    assert result.real().accuracy() > 400 and result.imag().contains_zero()
    return result.real()