"""Bloch elliptic dilogarithm at CM torsion points -- numberdb.org/T700.

Run with SageMath (or a compatible passagemath installation):
    sage -pip install numberdb
    sage -python generate.py             # verify all stored entries
    sage -python generate.py --publish   # submit entries to this draft only

Set NUMBERDB_API_KEY locally for writing or reading the private draft. The
campaign may supply it on stdin with NUMBERDB_KEY_FROM_STDIN=1. Submitting
entries does not offer, review or publish the draft.

D_q(u) = D(x) + sum_{n>=1} [D(x q^n) - D(q^n/x)], x=exp(2*pi*i*u).
The inversion law D(1/z)=-D(z) puts all arguments in the closed unit disk.
For |z|=t<1, |D(z)| <= t*(1+|log(t)|)/(1-t), by the power-series bounds
on Li_2 and log(1-z). With rho=|q|=exp(-L), a=r/N and c=M+1 +/- a,
the tail after n=M is at most the sum over the two c of
    rho^c/(1-rho^c) * ((1+L*c)/(1-rho)+L*rho/(1-rho)^2).
All evaluations and this tail bound use Arb balls; add_error expands the
result by an upper bound. The 128 extra working bits are checked by the
attached all-entry precision measurements, not used as a substitute for a
bound. No library approximation is wrapped in an artificial interval.

The sole selection source is this file. argument-selection.md explains its
arithmetic families and noncanonical finite torsion bound. checks.py provides
independent Fourier/Hurwitz, Mahler-integral and lattice-sign checks.
"""
import os
import sys
from functools import lru_cache
from math import gcd
import numberdb.sage as numberdb
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ

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


@lru_cache(maxsize=None)
def context(delta, digits):
    bits = numberdb.bits(digits, losing=WORKING_GUARD)
    R, C = RealBallField(bits), ComplexBallField(bits)
    L = R.pi()*R(-delta).sqrt()
    rho = (-L).exp()
    b = delta % 2
    q = rho if b == 0 else -rho
    return R, C, L, rho, b, q


def bloch_wigner(z, log_abs):
    """All calls have |z|<=1; z=1 is handled by its limiting value."""
    if z == 1:
        return z.real().parent()(0)
    ans = z.polylog(2).imag()
    if log_abs != 0:
        ans += log_abs*(1-z).log().imag()
    return ans


def elliptic_d(delta, N, r, s, digits=100, diagnostics=False):
    """Enclose D_q((r*tau+s)/N), for 0<=r<N, without using point oddness."""
    assert delta < 0 and delta % 4 in (0, 1)
    assert N > 0 and 0 <= r < N
    R, C, L, rho, b, q = context(delta, digits)
    a = QQ(r)/N
    theta = 2*R.pi()*R((QQ(s)-QQ(b*r)/2)/N)
    x = C(-L*a, theta).exp()
    if r == 0 and s % N == 0:
        x = C(1)
    total = bloch_wigner(x, -L*a)
    target = R(2)**(-R.precision()+24)
    M = 0
    qn = R(1)
    while True:
        M += 1
        qn *= q
        total += bloch_wigner(x*qn, -L*(M+a))
        total -= bloch_wigner(qn/x, -L*(M-a))
        bound = R(0)
        for c in (M+1+a, M+1-a):
            t = (-L*c).exp()
            bound += t/(1-t)*((1+L*c)/(1-rho)+L*rho/(1-rho)**2)
        if bound < target:
            break
        assert M < 1000, 'tail bound did not converge'
    value = total.add_error(bound)
    assert value.is_finite()
    if diagnostics:
        return value, {'terms_each_side':M, 'accuracy_bits':int(value.accuracy()),
                       'tail_bound':str(bound), 'working_bits':int(R.precision())}
    return value


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

    def enumerate(self):
        for delta in DISCRIMINANTS:
            b = delta % 2
            for N in range(3, 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
                        if (2*s-b*r) % N == 0:
                            continue
                        yield dict(zip(self.parameters, (delta,N,r,s)))

    def value(self, params, digits):
        delta,N,r,s = (int(params[k]) for k in self.parameters)
        return elliptic_d(delta,N,r,s,digits)


if __name__ == '__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
        os.environ['NUMBERDB_API_KEY'] = sys.stdin.read().strip()
    generator = Table()
    if '--publish' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
        print(generator.publish(message='Fill CM torsion values with certified q-series enclosures'))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)