back to table · edit · history · where entries came from · files · download
8884 bytes, as of the version from 2026-10-26 08:15 (current). Recorded here, not run.
"""Independent numerical checks for numberdb.org/T508.
sage -python checks.py # local generator, writes nothing
sage -python checks.py --stored # authenticated draft readback
Adapted from the published depth-four checker at
https://numberdb.org/files/T465/checks.py?raw=1 (CC BY-SA 4.0).
The checker does not call PARI for its reference values. It evaluates
iterated integrals split at 1/2 using convergent multiple-polylogarithm
series in Arb ball arithmetic. For a word w=0^(s1-1)1 ... 0^(sr-1)1,
let L(w) be its iterated integral from 0 to 1/2, with the leftmost letter
outermost. Reflection t -> 1-t and subdivision at 1/2 give
zeta(w) = sum_{w=uv} L(complement(reverse(u))) * L(v).
Every resulting nonempty word ends in 1, so L is a convergent series
sum_{n>=1} 2^-n n^-a1 H_{n-1}(a2,...,ar). Here H is the finite nested
harmonic sum with strictly descending indices. For all words used, r<=12.
The tail after N terms is bounded by sum_{n>N} n^11 2^-n. Its ratio is at
most q=((N+2)/(N+1))^11/2, giving the explicit bound
(N+1)^11 2^(-N-1)/(1-q). Every operation and this tail use 550-bit balls.
No result is accepted unless finite, positive and of relative radius <1e-125.
Controls compare the empty integral, L(1)=log(2), L(1,1)=log(2)^2/2,
zeta(2), zeta(3), zeta(2,1), and zeta(3,1) with independent Arb constants.
All table values must contain the narrow independent enclosures after
serialization to 100 digits. Sum identities are checked at every weight
6 through 12, along with zeta(2,1,1,1,1)=zeta(6) and zeta(2,2,2,2,2)=pi^10/11!.
"""
import json
import os
import sys
from decimal import Decimal
from functools import lru_cache
from math import factorial
from itertools import product
import numberdb.sage as numberdb
from sage.rings.real_arb import RealBallField
from sage.rings.real_mpfi import RealIntervalField
from sage.rings.real_mpfr import RealField
from sage.rings.rational_field import QQ
from numberdb._write import to_text
from generate import LengthFiveMZV, indices, working_values, WORKING_DIGITS
from check import exactness, exact_in_disguise, measure, names_its_rings, prose
BITS = 550
TERMS = 650
B = RealBallField(BITS)
I = RealIntervalField(BITS)
def word(index):
return tuple(bit for a in index for bit in [0] * (a - 1) + [1])
def index_of(w):
assert not w or w[-1] == 1
out, count = [], 1
for bit in w:
if bit:
out.append(count)
count = 1
else:
count += 1
return tuple(out)
def cuts(index):
w = word(index)
return [(index_of(tuple(1 - x for x in reversed(w[:i]))), index_of(w[i:]))
for i in range(len(w) + 1)]
class HalfIntegral:
def __init__(self, requested):
needed = {()}
for index in requested:
for u, v in cuts(index):
needed.update((u, v))
# Include explicit controls, not merely values derived from the table.
needed.update(((1,), (1, 1)))
suffixes = {()}
for a in needed:
for j in range(len(a)):
suffixes.add(a[j:])
harmonic = {a: B(0) for a in suffixes}
harmonic[()] = B(1)
polylog = {a: B(0) for a in needed}
polylog[()] = B(1)
ordered = sorted(suffixes - {()}, key=lambda a: -len(a))
powers = max((max(a) for a in suffixes if a), default=1)
zpower = B(1)
for n in range(1, TERMS + 1):
inverse = B(1) / B(n)
invpowers = [B(1)]
for a in range(1, powers + 1):
invpowers.append(invpowers[-1] * inverse)
zpower /= 2
for a in needed - {()}:
polylog[a] += zpower * invpowers[a[0]] * harmonic[a[1:]]
# Longer indices first: every suffix on the right is H_(n-1).
for a in ordered:
harmonic[a] += invpowers[a[0]] * harmonic[a[1:]]
q = (B(TERMS + 2) / B(TERMS + 1)) ** 11 / 2
assert q < 1
tail = B(TERMS + 1) ** 11 * B(2) ** (-TERMS - 1) / (1 - q)
for a in needed - {()}:
polylog[a] = polylog[a].add_error(tail.upper())
self.polylog = polylog
self.tail = tail
self.term_count = len(needed)
def __call__(self, index):
return sum((self.polylog[u] * self.polylog[v] for u, v in cuts(index)), B(0))
def narrow(x):
assert x.is_finite() and x > 0, str(x)
assert B(x.rad()) < abs(B(x.center())) * B('1e-125'), str(x)
def control(actual, expected):
narrow(actual)
narrow(expected)
assert actual.overlaps(expected), (str(actual), str(expected))
def decimal_claim(text):
"""Parse the stored decimal as center +/- exactly one last-place unit."""
d = Decimal(text)
n, den = d.as_integer_ratio()
c = QQ(n) / QQ(den)
ulp = QQ(10) ** int(d.as_tuple().exponent)
return I(c - ulp, c + ulp)
def stored_values():
# check.stored requires Django; this uses the authoritative public API.
t = numberdb.table('T508')
out = {}
for a, row1 in t['Numbers'].items():
for b, row2 in row1.items():
for c, row3 in row2.items():
for d, row4 in row3.items():
for e, value in row4.items():
out[tuple(map(int, (a, b, c, d, e)))] = value
return out
def run(stored=False):
controls = [(2,), (3,), (2, 1), (3, 1)]
h = HalfIntegral(list(indices()) + controls)
control(h.polylog[()], B(1))
control(h.polylog[(1,)], B(2).log())
control(h.polylog[(1, 1)], B(2).log() ** 2 / 2)
control(h((2,)), B(2).zeta())
control(h((3,)), B(3).zeta())
control(h((2, 1)), B(3).zeta())
control(h((3, 1)), B.pi() ** 4 / 360)
assert not h((2,)).overlaps(B(3).zeta())
assert not decimal_claim('1.000').overlaps(I(2))
print('Known-value and negative controls passed.', flush=True)
g = LengthFiveMZV()
values = stored_values() if stored else {
tuple(p[k] for k in g.parameters): g.value(p, 100) for p in g.enumerate()}
expected = {s for s in product(range(2, 9), range(1, 8), range(1, 8), range(1, 8), range(1, 8)) if sum(s) <= 12}
assert set(values) == expected and len(values) == 462
assert not names_its_rings(os.path.join(os.path.dirname(__file__), 'generate.py'))
assert not exactness(values)
assert not exact_in_disguise(values)
assert not prose(values)
evidence, serialized, intervals = [], {}, {}
for k, v in values.items():
text = v.get('number') if isinstance(v, dict) else v
text = to_text(text, digits=100)
ref = h(k)
narrow(ref)
claim = decimal_claim(text)
enclosure = I(ref)
assert claim.lower() < enclosure.lower() < enclosure.upper() < claim.upper(), k
ulp = QQ(10) ** int(Decimal(text).as_tuple().exponent)
assert not (claim + I(10 * ulp)).overlaps(enclosure), k
intervals[k] = claim
serialized[k] = {**v, 'number': text} if isinstance(v, dict) else text
evidence.append({'index': list(k), 'number': text,
'independent_lower': str(enclosure.lower()),
'independent_upper': str(enclosure.upper())})
# The equality is checked on actual serialized values, including readback.
for weight in range(6, 13):
total = sum((v for k, v in intervals.items() if sum(k) == weight), I(0))
assert I(B(weight).zeta()).overlaps(total), weight
for k, ref in [((2,1,1,1,1), B(6).zeta()),
((2,2,2,2,2), B.pi() ** 10 / factorial(11))]:
assert intervals[k].overlaps(I(ref)), k
R = RealField(800)
agreements = []
for k in indices():
lo, hi = [R(working_values(d)[k]) for d in WORKING_DIGITS]
relative = abs((lo - hi) / hi)
if relative:
agreements.append(-relative.log10())
report = {'tid':'T508','mode':'stored' if stored else 'local',
'entries':len(values), 'controls_passed':7,'negative_controls_passed':len(values)+2,
'independent_enclosures_contained':len(values),
'sum_identities_passed':7,'special_values_passed':2,
'working_digits':list(WORKING_DIGITS),
'minimum_agreement_digits':str(min(agreements)),
'independent_bits':BITS,'series_terms':TERMS,
'half_integrals':h.term_count,'tail_bound':str(h.tail),
'serialized_measurement':measure(serialized),
'values':evidence}
print('EVIDENCE_JSON='+json.dumps(report, sort_keys=True))
return report
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
key=sys.stdin.read().strip()
if '=' in key and key.split('=',1)[0].isupper():
key=key.split('=',1)[1].strip().strip('\'"')
numberdb.configure(api_key=key)
run(stored='--stored' in sys.argv)