"""Independent checks for numberdb.org/T500. sage -python checks.py # recompute, sends nothing sage -python checks.py --stored # also check all database decimal intervals A draft needs NUMBERDB_API_KEY; the campaign may instead pipe its key with NUMBERDB_KEY_FROM_STDIN=1. Download generate.py, check.py and euler-controls.json alongside this. The independent construction integrates the square wave sign(sin(pi*x)). Its antiperiodic primitives are E_r(x)/r! on [0,1]. All coefficients, Sturm root counts, extrema and C(n,k)^n are computed over QQ. This uses neither zeta values nor pi. Exact endpoint-power comparisons verify the interval claimed by every decimal stored in the database. """ import json import sys import math from pathlib import Path from fractions import Fraction import numberdb.sage as numberdb from sage.rings.rational_field import QQ from sage.rings.real_arb import RealBallField from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.version import version as sage_version from numberdb._write import to_text from check import exactness, measure, names_its_rings, agrees_with, exact_in_disguise, prose from generate import LandauKolmogorov, MAX_N, WORKING_GUARD, constant, favard, configure_from_stdin P = PolynomialRing(QQ, 'x') x = P.gen() def interior_roots(p): """Count distinct roots in (0,1) by an exact Sturm sequence.""" p = P(p) if not p: raise ValueError('The zero polynomial has no finite root count.') for a in (0, 1): while p.degree() > 0 and p(a) == 0: p, remainder = p.quo_rem(x-a) assert remainder == 0 if p.degree() == 0: return 0 seq = [p, p.derivative()] while seq[-1]: remainder = -(seq[-2].quo_rem(seq[-1])[1]) if not remainder: break seq.append(remainder) def variations(a): signs = [1 if q(a) > 0 else -1 for q in seq if q(a) != 0] return sum(a != b for a, b in zip(signs, signs[1:])) return variations(0) - variations(1) def spline_data(limit): f = P.one() polynomials = {0: f} norms = {0: QQ(1)} for n in range(1, limit+1): primitive = sum((c*x**(j+1)/QQ(j+1) for j, c in enumerate(f.list())), P.zero()) g = primitive - primitive(1)/2 assert g.derivative() == f assert g(0) == -g(1) assert g(1-x) == (-1)**n*g assert g(x+1)+g == 2*x**n/QQ(math.factorial(n)) assert interior_roots(g.derivative()) == (1 if n % 2 == 0 else 0) if n % 2 == 0: assert g.derivative()(QQ(1)/2) == 0 A = abs(g(QQ(1)/2 if n % 2 == 0 else 0)) assert A > 0 assert A == max(abs(g(a)) for a in (QQ(0), QQ(1)/2, QQ(1))) norms[n] = A polynomials[n] = g f = g return polynomials, norms def rational_power(norms, n, k): return norms[n-k]**n / norms[n]**(n-k) def rational_nth_root(q, n): a, a_exact = q.numerator().nth_root(n, truncate_mode=True) b, b_exact = q.denominator().nth_root(n, truncate_mode=True) return QQ(a)/b if a_exact and b_exact else None def decimal_bounds(text): """Exact NumberDB interval for a decimal; one unit in its last place.""" text = str(text) centre = Fraction(text) mantissa, sep, exponent = text.lower().partition('e') exponent = int(exponent) if sep else 0 places = len(mantissa.partition('.')[2]) if '.' in mantissa else 0 unit = Fraction(10)**(exponent-places) def rational(v): return QQ(v.numerator)/v.denominator return rational(centre-unit), rational(centre+unit) def validate_stored(document, norms): rows = document['Numbers'] if isinstance(rows, dict): rows = [dict({'params': {'n': n, 'k': k}}, **(entry if isinstance(entry, dict) else {'number': entry})) for n, group in rows.items() for k, entry in group.items()] expected = [(n,k) for n in range(2,MAX_N+1) for k in range(1,n)] actual = [(int(row['params']['n']), int(row['params']['k'])) for row in rows] assert actual == expected, 'Stored rows are missing, duplicated, or out of order.' widths = [] for row, (n,k) in zip(rows, actual): q = rational_power(norms,n,k) lo, hi = decimal_bounds(row['number']) assert lo > 0 and lo**n < q < hi**n, (n,k,'decimal interval excludes the root') assert hi-lo <= QQ(2)/10**99 widths.append(str(hi-lo)) return {'entries':len(rows),'ordered_complete_triangle':True, 'exact_endpoint_power_checks':len(rows),'largest_interval_width':max(widths,key=QQ)} def run(stored=False): # Controls are checked before trusting root counts or interval tests. assert interior_roots((x-QQ(1)/3)*(x-QQ(2)/3)) == 2 assert interior_roots(x*(x-1)) == 0 assert interior_roots(x-QQ(1)/2) == 1 assert interior_roots(x*x+1) == 0 lo,hi = decimal_bounds('1.4142135623730950488016887242096980785696718753769') assert lo**2 < 2 < hi**2 lo,hi = decimal_bounds('1.4142135623730950488016887242096980785696718753869') assert not lo**2 < 2 < hi**2 assert rational_nth_root(QQ(8)/27,3)==QQ(2)/3 assert rational_nth_root(QQ(2),2) is None polys,norms = spline_data(MAX_N) assert agrees_with({r:norms[r] for r in range(1,5)}, lambda r: {1:QQ(1)/2,2:QQ(1)/8,3:QQ(1)/24,4:QQ(5)/384}[r]) == [] source = json.loads(Path(__file__).with_name("euler-controls.json").read_text()) assert not agrees_with({r:polys[r]*math.factorial(r) for r in polys}, lambda r: P(source["Numbers"][str(r)])) R = RealBallField(512) precision = numberdb.bits(100, losing=WORKING_GUARD) values = {} balls = {} gen = LandauKolmogorov() for params in gen.enumerate(): n,k = params['n'],params['k'] v = gen.value(params,100) values[(n,k)] = v balls[(n,k)] = v['number'] if isinstance(v,dict) else v assert not exactness(values) assert not names_its_rings(__file__.replace('checks.py','generate.py')) assert not prose(values) for r in range(1,MAX_N+1): K = favard(r,precision) assert K.is_finite() and R(K.rad()) < R(10)**(-110) expected = R(norms[r]) * R.pi()**r assert expected.is_finite() and R(expected.rad()) < R(10)**(-140) assert K.overlaps(expected), ('Favard/spline disagreement',r) assert favard(0,precision) == 1 rational_rows=[] for (n,k),v in balls.items(): q=rational_power(norms,n,k) if rational_nth_root(q,n) is not None: rational_rows.append((n,k)) independent = R(q)**(QQ(1)/n) assert v.is_finite() and R(v.rad()) < R(10)**(-110) assert independent.is_finite() and R(independent.rad()) < R(10)**(-140) assert v.overlaps(independent), (n,k,'independent disagreement') assert (v**n-q).contains_zero(), (n,k,'power identity') assert v > 1 and v < R.pi()/2 assert rational_rows == [], rational_rows special={(2,1):QQ(2),(3,1):QQ(9)/8,(3,2):QQ(3), (4,1):QQ(512)/375,(4,2):QQ(36)/25,(4,3):QQ(24)/5, (5,1):QQ(1953125)/1572864,(5,2):QQ(125)/72} assert not agrees_with({key:rational_power(norms,*key) for key in special},lambda key:special[key]) # OEIS A245198's displayed 61 significant digits, not an unavailable b-file. oeis='1.0400419115259520572650284121789426931689026701866310548487955401' lo,hi=decimal_bounds(oeis) assert lo**3 < QQ(9)/8 < hi**3 out={'table':'T500','sage':sage_version,'numberdb':numberdb.__version__, 'range':{'n':[2,MAX_N],'k':'1..n-1'},'working_bits':precision, 'guard_bits':WORKING_GUARD,'measurement':measure(values), 'largest_ball_radius':str(max(v.rad() for v in balls.values())), 'exact_rational_rows':rational_rows,'independent_spline_comparisons':len(balls), 'corpus_euler_polynomial_comparisons':MAX_N+1,'favard_norm_comparisons':MAX_N,'published_small_case_checks':len(special), 'oeis_displayed_digits_check':True,'finite_and_narrow_controls':True, 'sturm_and_interval_negative_controls':True, 'exact_in_disguise_remarks':exact_in_disguise(values)} serialized = [{'params':{'n':str(n),'k':str(k)},'number':to_text(v,100,gen.format)} for (n,k),v in balls.items()] out['serialized_measurement'] = measure({(row['params']['n'],row['params']['k']):row['number'] for row in serialized}) out['local_serialized_intervals'] = validate_stored({'Numbers':serialized}, norms) out['computed_values'] = serialized if stored: configure_from_stdin() out['stored'] = validate_stored(numberdb.table('T500'), norms) verification=gen.verify(sample=None) assert verification.ok, str(verification) out['verify']=str(verification) return out if __name__ == '__main__': print(json.dumps(run('--stored' in sys.argv),indent=2,default=str))