"""Spectral determinants of Laplacians on spheres -- numberdb.org/T595.

Run with SageMath and the numberdb package:
    sage -pip install numberdb
    sage -python generate.py                    # verify every stored entry
    sage -python generate.py --fill             # fill the existing draft
The --fill operation requires NUMBERDB_API_KEY. It does not publish a draft.

Cunha and Freitas, arXiv:2404.12114v1, Theorems 2.7 and 2.8, give finite
formulas for the scalar unit-sphere Laplacian with its zero mode omitted.
All combinatorial coefficients and harmonic sums are exact rationals.
Arb encloses zeta derivatives, logarithms and exponentials throughout.

The default uses 150 working decimal digits for 100 output digits. The
attached checks compare this with 200 digits and measure retained accuracy.
The dimension cutoff and the invertible conventions are justified in
argument-selection.md. checks.py independently continues the spectral sum.
"""
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

MAX_DIMENSION = 20
WORKING_EXTRA_DIGITS = 50


@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 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()


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

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

    def value(self, params, digits):
        L = logdet(int(params['d']), max(150, digits+WORKING_EXTRA_DIGITS))
        if params['convention'] == 'det':
            result = L.exp()
        elif params['convention'] == 'logdet':
            result = L
        else:
            raise ValueError('Unknown determinant convention.')
        assert result.is_finite() and result.accuracy() >= numberdb.bits(digits+10)
        return result


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
    generator = SphereDeterminants()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Enclose both sphere determinant conventions in dimensions 1 through 20'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)