"""Independent numerical and algebraic checks for numberdb.org/T656. Run: sage -python checks.py [--stored] No source from the companion generator is used. The reciprocal equation is solved directly with exact rational signs. Controls run before comparisons. The only computation imported from generate.py is the output being tested. """ import json import sys from decimal import Decimal from fractions import Fraction from pathlib import Path import numberdb.sage as numberdb from numberdb._write import to_text from numberdb._generate import _stored_entries from sage.rings.rational_field import QQ from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing from sage.rings.real_arb import RealBallField from sage.rings.real_mpfi import RealIntervalField from sage.version import version as sage_version import mpmath INDEPENDENT_STEPS = 480 INDEPENDENT_BITS = 560 B = RealBallField(INDEPENDENT_BITS) I = RealIntervalField(INDEPENDENT_BITS) def rational(text): q = Fraction(str(text)) return QQ(q.numerator) / QQ(q.denominator) def written_bounds(text): text = str(text) center = rational(text) if '.' not in text and 'e' not in text.lower(): return center, center quantum = Decimal(1).scaleb(Decimal(text).as_tuple().exponent) ulp = rational(str(quantum)) return center - ulp, center + ulp def bounds_to_ball(lo, hi): return B((lo + hi) / QQ(2)).add_error(B((hi - lo) / QQ(2))) def assert_narrow(ball, tolerance=QQ(1)/QQ(10)**130): assert ball.is_finite() interval = I(ball) assert interval.absolute_diameter() < tolerance def encloses(text, lo, hi): a, b = written_bounds(text) assert a <= lo <= hi <= b, (text, 'does not contain certified bounds') def contains_ball(text, ball): assert_narrow(ball) interval = I(ball) encloses(text, interval.lower().exact_rational(), interval.upper().exact_rational()) def reciprocal_equation(n, g): # Clear the strictly positive denominators for g>1 in CVZ equation (31). return (3 - (n - 1)*(g - 1))*(g**n - 1) - 3*n*(g - 1) def independent_gamma(n): lo, hi = QQ(1) + QQ(1)/QQ(2*n), QQ(3) assert reciprocal_equation(n, lo) > 0 > reciprocal_equation(n, hi) for _ in range(INDEPENDENT_STEPS): mid = (lo + hi)/QQ(2) f = reciprocal_equation(n, mid) if f == 0: return mid, mid if f > 0: lo = mid else: hi = mid assert reciprocal_equation(n, lo) > 0 > reciprocal_equation(n, hi) return lo, hi def controls(): a, b = independent_gamma(2) assert a <= 2 <= b and b-a < QQ(1)/QQ(10)**140 assert reciprocal_equation(2, QQ(3)/2) > 0 assert reciprocal_equation(2, QQ(5)/2) < 0 p = B(QQ(1)/3) entropy = -p*p.log() - (1-p)*(1-p).log() known = B(3).log() - (QQ(2)/3)*B(2).log() assert_narrow(entropy); assert_narrow(known) assert entropy.overlaps(known) and entropy > B(QQ(6)/10) a, b = independent_gamma(3) gamma3 = (B(33).sqrt()+1)/4 assert_narrow(gamma3) assert bounds_to_ball(a,b).overlaps(gamma3) def algebraic_checks(n): R = PolynomialRing(QQ,'t'); t = R.gen() S = sum(t**j for j in range(n)) P = sum((3*j-n+1)*t**j for j in range(n)) assert P == 3*t*S.derivative()-(n-1)*S assert P(0) < 0 < P(1) rev = sum((3*j-n+1)*t**(n-1-j) for j in range(n)) assert (t-1)**2*rev == (3-(n-1)*(t-1))*(t**n-1)-3*n*(t-1) # All possible positive rational roots, by the rational root theorem. candidates = {QQ(a)/QQ(b) for a in range(1,n) if (n-1)%a==0 for b in range(1,2*n-1) if (2*(n-1))%b==0 and a < b} rational_roots = sorted(q for q in candidates if P(q)==0) assert rational_roots == ([QQ(1)/2] if n==2 else []) return P def main(stored=False): controls() from generate import MinimumPoint, configured_client client = configured_client() gen = MinimumPoint() raw = {f"{p['n']},{p['convention']}":gen.value(p,100) for p in gen.enumerate()} values = {k:to_text(v,100) for k,v in raw.items()} if stored: values = _stored_entries(numberdb.table('T656',client=client)) expected_keys = [f'{n},{c}' for n in range(2,33) for c in ('rho','gamma')] assert list(values)==expected_keys companion = _stored_entries(numberdb.table('T655',client=client)) entropy_table = _stored_entries(numberdb.table('T339',client=client)) algebraic_table = _stored_entries(numberdb.table('T35',client=client)) records=[]; min_bits=10000 for n in range(2,33): P = algebraic_checks(n) gl,gu = independent_gamma(n) if n == 2: assert gl <= 2 <= gu and reciprocal_equation(2, QQ(2)) == 0 gl = gu = QQ(2) rl,ru = QQ(1)/gu,QQ(1)/gl for convention,lo,hi in [('rho',rl,ru),('gamma',gl,gu)]: key=f'{n},{convention}'; text=values[key] encloses(text,lo,hi) a,b = written_bounds(text) value=raw[key] if n==2: assert a==b==(QQ(1)/2 if convention=='rho' else QQ(2)) else: vi = I(value) ratio = vi.lower().exact_rational() / (vi.upper().exact_rational()-vi.lower().exact_rational()) min_bits = min(min_bits,int(ratio.numerator().nbits()-ratio.denominator().nbits()-1)) assert value.is_finite() and I(value).absolute_diameter()0>reciprocal_equation(n,b) records.append({'n':n,'convention':convention,'number':text, 'independent_lower':str(lo),'independent_upper':str(hi)}) r=bounds_to_ball(rl,ru); g=bounds_to_ball(gl,gu) S=sum(r**j for j in range(n)); p=[r**j/S for j in range(n)] norm=sum(p); mean=sum(j*p[j] for j in range(n)) entropy=-sum(x*x.log() for x in p) variance=sum(j*j*p[j] for j in range(n))-mean**2 objective=S*r**(-QQ(n-1)/3) cvz=(g**n-1)/(g-1)*g**(-QQ(2*(n-1))/3) for value in [norm,mean,entropy,variance,objective,cvz,r*g]: assert_narrow(value) assert (1 in norm) and (QQ(n-1)/3 in mean) and variance>0 assert (1 in r*g) assert objective.overlaps(cvz) and entropy.overlaps(objective.log()) contains_ball(companion[f'{n},subrank'],objective) contains_ball(companion[f'{n},per-dimension'],objective/n) # The written reciprocal intervals must also overlap each other's inverse. a,b=written_bounds(values[f'{n},rho']); c,d=written_bounds(values[f'{n},gamma']) assert a*c <= 1 <= b*d if n==2: contains_ball(entropy_table['bernoulli,1/3,nats'],entropy) contains_ball(entropy_table['bernoulli,1/3,bits'],entropy/B(2).log()) if n==3: contains_ball(values['3,rho'],(B(33).sqrt()-1)/8) contains_ball(values['3,gamma'],(B(33).sqrt()+1)/4) contains_ball(algebraic_table['4,1,-2,2'],r) contains_ball(algebraic_table['2,-1,-4,2'],g) # Run precisely the approximation snippet intended for Programs. mp=mpmath.mp; mp.dps=110 errors=[] for n in range(2,33): stationary=lambda t:sum((3*j-n+1)*t**j for j in range(n)) rho=mp.findroot(stationary,(mp.mpf(0),mp.mpf(1)),solver='bisect', tol=mp.mpf('1e-105'),maxsteps=400) for c,v in [('rho',rho),('gamma',1/rho)]: a,b=written_bounds(values[f'{n},{c}']) q=rational(mp.nstr(v,108)); assert a <= q <= b or (n==2 and abs(q-a)