"""Values of the Brjuno function -- numberdb.org/T677.

Run with SageMath (install the numberdb Python package once):
    sage -pip install numberdb
    sage -python generate.py
The default command verifies every stored entry. Reading a private draft needs
NUMBERDB_API_KEY. This script neither creates nor offers nor publishes a table.
An authorized writer can use BrjunoValues().publish() to fill its existing draft.

The ordinary Gauss-map sum is evaluated by summing one exact quadratic cycle
and dividing by 1 minus its product. All logarithms and roots use Arb balls.
The 96 extra working bits are checked by checks.py against written last-place
intervals; they are not themselves an accuracy claim. See argument-selection.md
for the discriminant cutoff and checks.py for the independent convergent sum.
"""
import numberdb.sage as numberdb
from sage.rings.real_arb import RealBallField
from math import gcd, isqrt

DISCRIMINANT_BOUND = 100
WORKING_GUARD = 96


def triples(bound=DISCRIMINANT_BOUND):
    """Primitive A*t^2+B*t-K, with 0<x<1 and conjugate x'<-1."""
    out = []
    for A in range(1, isqrt(bound) + 1):
        for B in range(1, isqrt(bound) + 1):
            for K in range(1, bound // 4 + 1):
                D = B*B + 4*A*K
                if D > bound or isqrt(D)**2 == D:
                    continue
                if gcd(gcd(A, B), K) != 1:
                    continue
                if A+B-K > 0 and A-B-K < 0:
                    out.append((A, B, K))
    return sorted(out, key=lambda t: (t[1]**2+4*t[0]*t[2], t))


def identity(t):
    return 'quad_%d_%d_m%d' % t


def parse_identity(key):
    import re
    match = re.fullmatch(r'quad_(\d+)_(\d+)_m(\d+)', key)
    if not match:
        raise ValueError('Expected a primitive quadratic root identity')
    return tuple(map(int, match.groups()))


def step(t):
    A, B, K = t
    D = B*B + 4*A*K
    m = (B + isqrt(D)) // (2*K)
    return (K, 2*m*K-B, A+B*m-K*m*m), m


def cycle(t):
    orbit, digits, seen = [], [], set()
    cur = t
    while cur not in seen:
        seen.add(cur)
        orbit.append(cur)
        cur, m = step(cur)
        digits.append(m)
    if cur != t:
        raise ValueError('Argument has a preperiod; not part of this selection')
    return orbit, digits


def root(t, field):
    A, B, K = t
    return field(2*K) / (field(B) + field(B*B+4*A*K).sqrt())


def brjuno(t, digits=100):
    field = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
    product, total = field(1), field(0)
    orbit, _ = cycle(t)
    for cur in orbit:
        x = root(cur, field)
        total += product * (-x.log())
        product *= x
    result = total / (1-product)
    if not result.is_finite() or not result > 0:
        raise ArithmeticError('Brjuno enclosure is not finite and positive')
    return result


class BrjunoValues(numberdb.Generator):
    table = 'T677'
    parameters = ('x',)
    type = 'R'
    digits = 100
    rigour = 'proven'

    def enumerate(self, discriminant_bound=DISCRIMINANT_BOUND):
        for t in triples(discriminant_bound):
            yield {'x': identity(t)}

    def value(self, params, digits):
        t = parse_identity(params['x'])
        A, B, K = t
        D = B*B+4*A*K
        result = {'number': brjuno(t, digits),
                  'param-latex': r'$\frac{- %d+\sqrt{%d}}{%d}$' % (B,D,2*A)}
        if t == (1,1,1):
            result['comment'] = 'The argument is HREF{Golden_ratio#phi_inv}[the reciprocal of the golden ratio], with continued fraction $[0;\\overline{1}]$.'
        elif t == (1,2,1):
            result['comment'] = 'The argument is the reciprocal of HREF{Markov_quadratic_irrationals#1,1,2}[$1+\\sqrt2$], with continued fraction $[0;\\overline{2}]$.'
        return result


if __name__ == '__main__':
    import sys
    report = BrjunoValues().verify(sample=None)
    print(report)
    sys.exit(0 if report.ok else 1)