"""Values of %s at rational arguments -- numberdb.org/TID_%s

Every x = a/b in lowest terms with b <= 6 and 0 <= x <= 10.

Run it with SageMath:

    $ sage -pip install numberdb          # once
    $ sage -python generate.py            # check the table against this code
    $ sage -python generate.py --publish  # send it, with NUMBERDB_API_KEY set

Values are computed as balls with arb.

One function per table, as T192 and T193 are the sine and cosine
integrals and T198 and T202 the two Airy functions. These six were one
table until it was split: how each is written in terms of the others
is in `Formulas`, and that they belong together is in `Similar
tables`, which is where a relation can be written down.

The grid is wider than it was, 0 <= x <= 10 rather than 0 < x <= 5.
Zero is the one argument where every one of the six is exact, and the
table that left it out left out the only entry a reader can check by
hand.
"""

import os
import sys

import numberdb.sage as numberdb
from sage.rings.complex_arb import ComplexBallField
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField


#: Which of the six this table holds.
FUNCTION = "erf"

#: Its value at x = 0, exact, and stored as the integer rather than as
#: a hundred places of one.
AT_ZERO = "0"

# Every a/b in lowest terms with b <= 6, and 0 <= x <= 10. A hundred and
# twenty-one of them.
#
# This reverses a decision, and the decision it reverses had a real argument
# in it, so here is that argument and the answer to it. The grid was bounded
# height with b <= 18, and the note that replaced it said: "that is counting
# rather than choosing: erf(1.96) is a number people arrive with and erf(17/18)
# is not, and a denominator bound admits the second to reach the first." Both
# halves of that are true. b <= 18 *was* chosen for how many entries it made,
# and 17/18 *is* not a number anybody writes.
#
# The mistake was the remedy. A bound picked for its entry count is fixed by
# picking the bound for its arguments instead -- b <= 6 admits 1/6 and 5/6 and
# stops -- not by moving to a grid that is round in base ten. Two decimal
# places is not a property of erf; it is a property of how we write numbers.
#
# And the 1.96 argument answers itself once you ask who is searching. Search
# here is by value, not by parameter: nobody asks this database for erf(1.96),
# somebody has 0.99442... and wants to know what it is. A reader who evaluated
# erf at 1.96 already knows what it is and is not asking. The numbers that
# arrive unlabelled are the ones with structure behind them, and those sit at
# arguments a person would write down: a half, a third, two fifths.
#
# The cost of getting this wrong compounds: the two-decimal grid here was
# quoted as precedent by later critiques -- "T197 does the same with its
# two-decimal grid, so it is the site's behaviour and not this table's" -- and
# by 2026-09-23 fifty-one tables had followed it.
DENOMINATOR_BOUND = 6
MAX_ARGUMENT = 10

# Bits of working precision beyond what the written digits need.
#
# `verify` recomputes every entry and compares, so a guard too small
# for some argument fails there rather than quietly rounding; it is
# stated as a knob and checked as a result.
WORKING_GUARD = 64


def _key_from_stdin():
    if os.environ.get("NUMBERDB_KEY_FROM_STDIN") != "1":
        return
    token = sys.stdin.read().strip()
    if "=" in token and token.split("=", 1)[0].isupper():
        token = token.split("=", 1)[1].strip().strip("'\"")
    if token:
        os.environ["NUMBERDB_API_KEY"] = token


def _arguments(bound=DENOMINATOR_BOUND, maximum=MAX_ARGUMENT):
    #Every a/b in lowest terms with b <= bound, in increasing order. Zero
    #included: see the note at the top -- it is the one argument where every
    #member of the family is exact, and the only entry a reader can check by
    #hand.
    seen = set()
    for denominator in range(1, int(bound) + 1):
        for numerator in range(0, int(maximum) * denominator + 1):
            value = QQ(numerator) / QQ(denominator)
            if value > maximum:
                break
            seen.add(value)
    for value in sorted(seen):
        yield str(value)


def _real_field(digits):
    return RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))


def _complex_field(digits):
    return ComplexBallField(numberdb.bits(digits, losing=WORKING_GUARD))


def _fresnel(function, x_text, digits):
    field = _complex_field(digits)
    i = field.gen(0)
    one = field(1)
    z = field(QQ(x_text))
    argument = field.pi().sqrt() * (one - i) * z / field(2)
    value = (one + i) * argument.erf() / field(2)
    component = value.imag() if function == "fresnel-S" else value.real()
    if not component.is_finite():
        raise ArithmeticError("computed a non-finite ball for %s(%s)"
                              % (function, x_text))
    return component


def _value_ball(function, x_text, digits):
    if function in ("fresnel-S", "fresnel-C"):
        return _fresnel(function, x_text, digits)

    field = _real_field(digits)
    x = field(QQ(x_text))
    if function == "erf":
        value = x.erf()
    elif function == "erfc":
        value = field(1) - x.erf()
    elif function == "erfi":
        value = x.erfi()
    elif function == "dawson":
        value = field.pi().sqrt() * (-x * x).exp() * x.erfi() / field(2)
    else:
        raise ValueError("unknown error-function-family member %r"
                         % (function,))
    if not value.is_finite():
        raise ArithmeticError("computed a non-finite ball for %s(%s)"
                              % (function, x_text))
    return value


class ErrorFunctionValues(numberdb.Generator):

    table = os.environ.get("NUMBERDB_TABLE") or "T197"
    parameters = ("x",)
    type = "R"
    digits = 100
    rigour = "proven"

    def enumerate(self, bound=DENOMINATOR_BOUND, maximum=MAX_ARGUMENT):
        for x in _arguments(bound, maximum):
            yield {"x": x}

    def value(self, params, digits):
        x_text = str(params["x"])
        if QQ(x_text) == 0:
            #A theorem, not a measurement: every one of the six is 0 at the
            #origin except erfc, which is 1. Checked against the computed
            #ball all the same, so the exact value is verified rather than
            #asserted.
            ball = _value_ball(FUNCTION, x_text, digits)
            if not ball.contains_exact(QQ(AT_ZERO)):
                raise ValueError("%s(0) is written as %s and the computation "
                                 "does not agree" % (FUNCTION, AT_ZERO))
            return ZZ(AT_ZERO)
        return _value_ball(FUNCTION, x_text, digits)


if __name__ == "__main__":
    _key_from_stdin()
    generator = ErrorFunctionValues()
    if os.environ.get("NUMBERDB_PUBLISH") == "1" or "--publish" in sys.argv:
        print(generator.publish(
            message="values of erf(x) at rational arguments"))
    else:
        report = generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)