"""Values of the Riesz function -- numberdb.org/T645.

Run with SageMath or passagemath and numberdb installed:
    sage -python generate.py                  # verify every stored entry
    sage -python generate.py --fill-draft     # write entries to the existing draft
The write needs NUMBERDB_API_KEY. It does not offer or publish the draft.

R(x) = x sum_{j>=0} (-x)^j/(j! zeta(2j+2)). The default selection is
x=0,1,...,100. See argument-selection.md for its two distinct roles.

For x>=0, 1/zeta(2j+2)<1. After M terms, if M+1>x, the absolute tail is at
most x^(M+1)/M! / (1-x/(M+1)). Choose M with this exact rational majorant
below 10^(-digits-35). Arb propagates the rounding error in the finite sum;
the returned ball includes the omitted tail. At x=0 the defining factor x
proves the exact value zero. No assertion about the Riemann hypothesis is used.

Working precision is bits(digits)+2*ceil(x)+80. The x-dependent guard covers
cancellation of order exp(x); accuracy is also checked before returning.
The default range was measured in checks-evidence.json. The value method
accepts other nonnegative rational arguments, without adding them to the
selection. The attached checks use an independent accelerated Mobius sum,
exact Bernoulli numbers, and a Machin-series enclosure of pi.
"""
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField

GUARD_FIXED = 80
TAIL_EXTRA_DIGITS = 35


def truncation(x, digits):
    """Return M and a rational bound for terms j>=M, including outer x."""
    x = QQ(x)
    target = QQ(1) / ZZ(10)**(digits + TAIL_EXTRA_DIGITS)
    M, term = 0, QQ(1)
    while True:
        if M + 1 > x:
            bound = x * term / (1 - x / (M + 1))
            if bound < target:
                return M, bound
        M += 1
        term *= x / M


@lru_cache(maxsize=8)
def reciprocal_zetas(prec, count):
    F = RealBallField(prec)
    return tuple(1 / F(2*j + 2).zeta() for j in range(count))


def riesz(x, digits=100):
    x = QQ(x)
    if x < 0:
        raise ValueError('This generator evaluates nonnegative rational arguments.')
    if x == 0:
        return ZZ(0)
    prec = numberdb.bits(digits, losing=2*int(x.ceil()) + GUARD_FIXED)
    F = RealBallField(prec)
    M, tail = truncation(x, digits)
    inverses = reciprocal_zetas(prec, M)
    term, total = F(1), F(0)
    for j, inverse in enumerate(inverses):
        total += term * inverse
        term *= -F(x) / (j + 1)
    result = F(x) * total + F(0, tail)
    if not result.is_finite() or result.accuracy() < numberdb.bits(digits) + 16:
        raise ArithmeticError('Insufficient certified relative accuracy.')
    return result


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

    def enumerate(self, limit=100):
        for x in range(limit + 1):
            yield {'x': str(x)}

    def value(self, params, digits):
        return riesz(QQ(params['x']), digits)


if __name__ == '__main__':
    g = RieszValues()
    if '--fill-draft' in sys.argv:
        print(g.publish(message='fill the Riesz draft with certified values'))
    else:
        report = g.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)