"""Dirichlet spectral determinants of regular polygons -- numberdb.org/T620.

Run with SageMath or passagemath and numberdb installed:
    sage -pip install numberdb           # once
    sage -python generate.py             # verify every stored entry
    sage -python generate.py --fill      # fill the existing draft, API key required

The operator is -d_x^2-d_y^2 with Dirichlet boundary conditions on the
regular n-gon of circumradius one. Entries are det and its natural log,
in that order for n=3,...,30. See argument-selection.md for the selection.

Aurell--Salomonson (hep-th/9304031), (23)--(24), gives the polygon formula.
Kalvin (1910.00104v2), (3.15) and Lemma A.1, gives a finite corner formula.
All divisions of exact coefficients use QQ. Arb encloses the exact rational
arguments and carries errors through log-Gamma, zeta derivative, and exp.
The 180 guard bits are checked by checks.py, including the decimal intervals.
The independent integral and spectral controls are in checks.py.

This script writes entries only when --fill is explicit. It cannot offer,
review, or make a table public. Authentication may be piped on stdin with
NUMBERDB_KEY_FROM_STDIN=1; otherwise use NUMBERDB_API_KEY.
"""
import os
import sys
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 sage.rings.complex_arb import ComplexBallField

WORKING_GUARD = 180


@lru_cache(maxsize=None)
def ingredients(prec):
    field = RealBallField(prec)
    z = ComplexBallField(prec)(-1).zetaderiv(1)
    assert z.real().is_finite() and z.imag().is_finite() and z.imag().contains_zero()
    return field, z.real()


def corner(a, prec):
    """F(a), with a positive exact rational; no floating-point coefficients."""
    a = QQ(a)
    p, q = a.numerator(), a.denominator()
    field, zp = ingredients(prec)
    half = QQ(1) / 2
    dedekind = sum(((QQ(j) / p - half) * (QQ((j*q) % p) / p - half)
                   for j in range(1, p)), QQ(0))
    b = zp / (p*q) - field(q).log() / (12*p*q)
    b += (QQ(1)/4 + dedekind) * field(QQ(q)/p).log()
    for k in range(1, p):
        b += (half - QQ(k)/p) * field(QQ((k*q) % p)/p).log_gamma()
    for j in range(1, q):
        b += (half - QQ(j)/q) * field(QQ((j*p) % q)/q).log_gamma()
    f = b - a*zp + (a-1/a)/12 * field(2).log()
    f -= (a-1)/4 * (2*field.pi()).log()
    assert f.is_finite()
    return f


@lru_cache(maxsize=None)
def polygon_logdet(n, prec):
    n = int(n)
    if n < 3:
        raise ValueError('The number of sides must be at least three')
    field, _ = ingredients(prec)
    h = QQ(1)/n
    a = 1-2*h
    c = QQ(n-1)/(6*(n-2))
    log_scale = (field(1+h).log_gamma() + field(a).log_gamma()
                 - field(1-h).log_gamma())
    result = 2*c*log_scale + field(n).log()/(3*(n-2)) - n*corner(a, prec)
    assert result.is_finite()
    return result


class PolygonDirichletDeterminants(numberdb.Generator):
    table = 'T620'
    parameters = ('n', 'convention')
    type = 'R'
    digits = 100
    rigour = 'proven'

    def enumerate(self, maximum=30):
        for n in range(3, maximum+1):
            for convention in ('det', 'logdet'):
                yield {'n': str(n), 'convention': convention}

    def value(self, params, digits):
        value = polygon_logdet(int(params['n']), numberdb.bits(digits, losing=WORKING_GUARD))
        if params['convention'] == 'det':
            return value.exp()
        if params['convention'] == 'logdet':
            return value
        raise ValueError('Unknown convention')


def authenticated_client():
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        key = sys.stdin.read().strip()
        if '=' in key and key.split('=', 1)[0].replace('_', '').isupper():
            key = key.split('=', 1)[1].strip().strip('\"\'')
        return numberdb.Client(api_key=key)
    return numberdb.Client()


if __name__ == '__main__':
    gen = PolygonDirichletDeterminants()
    client = authenticated_client()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(gen.publish(client=client, message='Fill 56 certified polygon determinants and logarithms'))
    else:
        report = gen.verify(sample=None, client=client)
        print(report)
        sys.exit(0 if report.ok else 1)