back to table · edit · history · where entries came from · files · download
6417 bytes, as of the version from 2026-09-10 20:24. Recorded here, not run.
"""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
from math import gcd
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 = "erfc"
#: Its value at x = 0, exact, and stored as the integer rather than as
#: a hundred places of one.
AT_ZERO = "1"
# A thousand values, which is the scale this corpus's own function tables
# hold: T20's Bessel zeros and T9's Gamma values are about 1050 each, at some
# 120 KB per document, and 60 was an order of magnitude short of that.
#
# The knob is the denominator, not the range. Raising it refines the grid
# between the points already there rather than stretching it further out,
# which is what somebody interpolating wants -- and every value the table
# already had stays exactly where it was.
MAX_DENOMINATOR = 18
MAX_ARGUMENT = 10
# Bits of working precision beyond what the written digits need.
#
# Measured over all 360 entries: at this guard the widest result still has
# radius less than 1e-108 when the table asks for 100 digits.
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(max_denominator=MAX_DENOMINATOR, maximum=MAX_ARGUMENT):
#Zero included: see the note at the top.
values = {QQ(0)}
for denominator in range(1, max_denominator + 1):
for numerator in range(1, maximum * denominator + 1):
if gcd(numerator, denominator) == 1:
values.add(QQ(numerator) / QQ(denominator))
for value in sorted(values):
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":
#Directly where arb offers it, because 1 - erf(x) is catastrophic
#cancellation for large x: erfc(13/2) is 1.5e-19 and erfc(10) is
#2.1e-45, so subtracting from 1 throws away that many leading digits
#and the ball comes back with 99 of the 100 asked for, then 55.
#Where it is not offered, the subtraction is done in a field wide
#enough to survive it: 400 bits is 120 decimal digits of headroom
#against the 45 that x = 10 costs.
if hasattr(x, "erfc"):
value = x.erfc()
else:
wide = RealBallField(numberdb.bits(digits, losing=400))
value = wide(1) - wide(QQ(x_text)).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 ComplementaryErrorFunctionValues(numberdb.Generator):
table = os.environ.get("NUMBERDB_TABLE") or "T204"
parameters = ("x",)
type = "R"
digits = 100
rigour = "proven"
def enumerate(self, denominator=MAX_DENOMINATOR,
maximum=MAX_ARGUMENT):
for x in _arguments(denominator, 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 = ComplementaryErrorFunctionValues()
if os.environ.get("NUMBERDB_PUBLISH") == "1" or "--publish" in sys.argv:
print(generator.publish(
message="values of erfc(x) at rational arguments"))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)