back to table · edit · history · where entries came from · files · download
9459 bytes, as of the version from 2026-11-01 17:38 (current). Recorded here, not run.
"""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()<QQ(1)/10**120
vi=I(value)
encloses(text,vi.lower().exact_rational(),vi.upper().exact_rational())
assert len(Decimal(text).as_tuple().digits)==100
if convention=='rho':
assert P(a)<0<P(b)
else:
assert reciprocal_equation(n,a)>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)<QQ(1)/10**103)
report={'table':'T656','mode':'stored' if stored else 'local',
'skill_manifest':'agents/runs/20261101T171341Z-build-skills/manifest.json',
'skill_sha256':'a2b28d9953264d7b0ca62663785c2952190c5be6d6d3597fd695d9a54d606278',
'entries_checked':len(values),'exact_rationals':2,'certified_decimals':60,
'minimum_generator_relative_accuracy_bits':min_bits,
'independent_bisection_steps':INDEPENDENT_STEPS,'independent_ball_bits':INDEPENDENT_BITS,
'control':'n=2 gamma=2, known nonzero Bernoulli entropy, n=3 radical',
'checks':['exact stationary and reciprocal polynomial identities for n=2..32',
'exhaustive rational-root theorem candidates',
'exact endpoint signs for every printed irrational interval',
'independently computed entire root brackets inside all printed intervals',
'normalization, mean, entropy, strictly positive variance, reciprocals',
'all 62 T655 companion constants and normalizations',
'both T339 Bernoulli(1/3) entropy units',
'both n=3 quadratic-algebraic entries in T35',
'Programs snippet on every n=2..32'],
'versions':{'sage':sage_version,'numberdb':numberdb.__version__,
'python':sys.version.split()[0],'mpmath':mpmath.__version__},
'records':records}
print(json.dumps(report,indent=2))
if __name__=='__main__':
main(stored='--stored' in sys.argv)