"""Independent numerical checks for numberdb.org/T465. sage -python checks.py # local generator, writes nothing sage -python checks.py --stored # authenticated draft readback 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 5 through 12, along with zeta(2,1,1,1)=zeta(5) and zeta(2,2,2,2)=pi^8/9!. """ 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 LengthFourMZV, 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('T465') out = {} for a, row1 in t['Numbers'].items(): for b, row2 in row1.items(): for c, row3 in row2.items(): for d, value in row3.items(): out[tuple(map(int, (a, b, c, d)))] = 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) g = LengthFourMZV() 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, 10), range(1, 9), range(1, 9), range(1, 9)) if sum(s) <= 12} assert set(values) == expected and len(values) == 330 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 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(5, 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), B(5).zeta()), ((2,2,2,2), B.pi() ** 8 / factorial(9))]: 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':'T465','mode':'stored' if stored else 'local', 'entries':len(values), 'controls_passed':7, 'independent_enclosures_contained':len(values), 'sum_identities_passed':8,'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)