"""Elliptic trilogarithm L31 at CM torsion points -- numberdb.org/T702.

Install numberdb in an existing SageMath environment:
    sage -pip install numberdb
Run with SageMath or passagemath and numberdb installed:
    sage -python generate.py                 # verify every stored entry
    sage -python generate.py --fill          # send entries to this draft
The latter requires NUMBERDB_API_KEY. It cannot publish or offer a draft.

Samart, arXiv:1309.7730v2, section 1 fixes the single-valued normalization.
The finite selection is explained in argument-selection.md. For 0<t<1,
|L3(z)| <= t/(1-t)*(1+|log(t)|+log(t)^2/3), t=|z|.
Each of the two omitted q-tails is bounded by summing this quadratic times
a geometric series. Its bound is added to the Arb result, including all
rounding and argument uncertainty. The boundary x=1 is zeta(3); for any
unit-modulus x the zero logarithmic factors are omitted exactly.

128 extra working bits and a tail target 24 decimal places beyond the
requested digits give substantially more than 100 significant digits across
the default selection; checks.py records the measured minimum accuracy.
"""
import os
import sys
from math import gcd
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

DISCRIMINANTS = (-3, -4, -7, -8, -11, -12, -16, -19, -27, -28, -43, -67, -163)
WORKING_GUARD = 128
TAIL_GUARD = 24


def tail_bound(R, A, v):
    """Bound sum_{k>=0} |L3(z_k)|, where |z_k|=R^(v+k)."""
    t = (-A*v).exp()
    S0 = 1/(1-R)
    S1 = R/(1-R)**2
    S2 = R*(1+R)/(1-R)**3
    return t/(1-t)*((1+A*v+A*A*v*v/3)*S0
                   +(A+2*A*A*v/3)*S1+A*A*S2/3)


def compute(params, digits=100, diagnostics=False):
    D,N,r,s = (int(params[k]) for k in ('Delta','N','r','s'))
    if D >= 0 or D % 4 not in (0,1) or N < 1:
        raise ValueError('Expected a negative quadratic-order discriminant and N>=1')
    # Accept redundant coordinates for checking periodicity, but default
    # enumeration is reduced and of exact order N.
    r %= N; s %= N
    B=RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
    C=ComplexBallField(B.precision())
    b=int(D%4 != 0)
    A=B.pi()*B(-D).sqrt()
    R=(-A).exp()
    q=C((-1 if b else 1)*R)
    a=QQ(r)/N
    x=C(-A*a, B.pi()*(QQ(2*s-b*r)/N)).exp()
    def L(z, v):
        logmod=-A*v
        if v==0:
            if s==0: return B(3).zeta()
            return z.polylog(3).real()
        return (z.polylog(3)-logmod*z.polylog(2)
                -logmod**2*(1-z).log()/3).real()
    total=L(x,a)
    qpower=C(1)
    target=B(10)**(-digits-TAIL_GUARD)
    for K in range(0,10000):
        bound=tail_bound(R,A,K+1+a)+tail_bound(R,A,K+1-a)
        if total.is_finite() and not total.contains_zero() and bound < abs(total)*target:
            result=total.add_error(bound)
            if not result.is_finite() or result.accuracy() < numberdb.bits(digits):
                raise ArithmeticError('Insufficient relative enclosure accuracy')
            if diagnostics:
                return result, {'terms_each_side':K,'relative_accuracy_bits':int(result.accuracy()),'tail_bound':str(bound)}
            return result
        qpower*=q
        n=K+1
        total+=L(qpower*x,n+a)+L(qpower/x,n-a)
    raise ArithmeticError('Tail target not reached; investigate an exact zero or cancellation')


class Table(numberdb.Generator):
    table='T702'
    parameters=('Delta','N','r','s')
    type='R'
    digits=100
    rigour='proven'

    def enumerate(self, max_order=6, discriminants=DISCRIMINANTS):
        for D in discriminants:
            for N in range(1,max_order+1):
                for r in range(N):
                    for s in range(N):
                        if gcd(gcd(N,r),s)!=1:
                            continue
                        if (r,s)>((-r)%N,(-s)%N):
                            continue
                        yield dict(Delta=D,N=N,r=r,s=s)

    def value(self, params, digits):
        return compute(params,digits)


if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        numberdb.configure(api_key=sys.stdin.read().strip())
    g=Table()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
        print(g.publish(message='Fill L31 CM torsion values with certified tails and independent checks'))
    else:
        report=g.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)