"""Independent checks for Euler-Bernoulli beam roots -- numberdb.org/T486. Download generate.py and check.py beside this file, then run: sage -python checks.py # check locally generated values sage -python checks.py --stored # check database values; draft needs key Every nonzero decimal is tested at BOTH exact endpoints of the interval its last digit denotes. A 4x4 boundary determinant is evaluated by its 24 signed permutation products in Arb at 768 bits, independently of the scalar equations used by generate.py. The boundary rows come directly from derivatives of cos(z*x), sin(z*x), cosh(z*x), sinh(z*x) on [0,1]. Row factors z^d are removed for nonzero z, which cannot alter zeros or endpoint signs. A second computation solves this boundary determinant with mpmath at 150 and 200 decimal digits, checking 90 places by nearest-integer comparison. Controls include the exact pinned-pinned determinant -4*sin(z)*sinh(z), a known pinned-pinned root pi, a nonroot at z=1, exact cubic nullities at z=0, and several low roots from the proposal's source checks. This is a numerical verification program, not an act of human review. """ import itertools import json import os import sys from decimal import Decimal from fractions import Fraction import numberdb.sage as numberdb from numberdb._write import Entries from sage.rings.integer_ring import ZZ from sage.rings.rational_field import QQ from sage.rings.real_arb import RealBallField import mpmath as mp from check import exactness, measure, agrees_with, names_its_rings from generate import BeamRoots, BOUNDARIES, MAX_N ORDERS = {'clamped': (0, 1), 'pinned': (0, 2), 'free': (2, 3)} PERMS = [(p, (-1)**sum(p[i] > p[j] for i in range(4) for j in range(i+1,4))) for p in itertools.permutations(range(4))] def det4(rows): return sum(sign * rows[0][p[0]] * rows[1][p[1]] * rows[2][p[2]] * rows[3][p[3]] for p, sign in PERMS) def boundary_determinant(boundary, z, trig): # This code shares no characteristic equations or root brackets with generate.py. sine, cosine, sinh, cosh = trig left, right = boundary.split('-') rows = [] for endpoint, support in [(0, left), (1, right)]: t = z * endpoint c, s, ch, sh = cosine(t), sine(t), cosh(t), sinh(t) for d in ORDERS[support]: rows.append(((c, -s, -c, s)[d], (s, c, -s, -c)[d], ch if d % 2 == 0 else sh, sh if d % 2 == 0 else ch)) return det4(rows) def exact_nullity(boundary): # At beta=0 use polynomial solutions, since the trigonometric basis degenerates. rows = [] for endpoint, support in zip((0, 1), boundary.split('-')): for d in ORDERS[support]: row = [] for j in range(4): coeff = 1 for h in range(d): coeff *= j-h row.append(QQ(0) if j < d else QQ(coeff) * QQ(endpoint)**(j-d)) rows.append(row) rank = 0 for col in range(4): pivot = next((i for i in range(rank,4) if rows[i][col]), None) if pivot is None: continue rows[rank], rows[pivot] = rows[pivot], rows[rank] pivot_value = rows[rank][col] rows[rank] = [v / pivot_value for v in rows[rank]] for i in range(4): if i != rank: factor = rows[i][col] rows[i] = [a-factor*b for a,b in zip(rows[i],rows[rank])] rank += 1 return 4-rank def rational(text): f = Fraction(text) return QQ(f.numerator) / QQ(f.denominator) def decimal_interval(text): c = rational(text) if '.' not in text and 'e' not in text.lower(): return c, c exponent = Decimal(text).as_tuple().exponent u = QQ(10)**exponent return c-u, c+u def main(stored=False): g = BeamRoots() rb = RealBallField(768) arb_trig = (lambda x:x.sin(), lambda x:x.cos(), lambda x:x.sinh(),lambda x:x.cosh()) mp_trig = (mp.sin,mp.cos,mp.sinh,mp.cosh) mp.mp.dps = 200 # Controls first. In the pinned-pinned matrix, add row 1 to row 2. # Expanding the resulting [0,0,2,0] row gives -4*sin(z)*sinh(z). assert det4([[1,0,1,0],[-1,0,1,0],[1,2,3,4],[-1,-2,3,4]]) == -32 control = boundary_determinant('pinned-pinned', rb(1), arb_trig) closed = -4 * rb(1).sin() * rb(1).sinh() assert control.is_finite() and control < 0 assert (control-closed).contains_zero() assert QQ(control.rad()) < QQ(1)/10**200 root_control = boundary_determinant('pinned-pinned', rb.pi(), arb_trig) assert root_control.is_finite() and root_control.contains_zero() assert QQ(root_control.rad()) < QQ(1)/10**200 static = {b:exact_nullity(b) for b in BOUNDARIES} assert static == {'clamped-free':0,'clamped-clamped':0,'clamped-pinned':0, 'pinned-pinned':0,'free-free':2,'pinned-free':1} print('Controls passed: nonzero determinant, known pi root, exact zero-mode nullities.') if stored: doc = numberdb.table('T486') block = doc['Numbers'] if isinstance(block,list): records=block else: records=[] for boundary, group in block.items(): for n, node in group.items(): fields=node if isinstance(node,dict) else {'number':node} records.append({'params':{'boundary':boundary,'n':str(n)},**fields}) print('Read back',len(records),'records from T486.') else: values = {} entries = Entries('boundary','n') min_accuracy = 1000000 for p in g.enumerate(): v = g.value(p,g.digits) values[(p['boundary'],int(p['n']))] = v data = v if isinstance(v,dict) else {'number':v} x = data['number'] if hasattr(x,'rad'): relative = abs(QQ(x.mid())) / QQ(x.rad()) lower_bits = relative.numerator().nbits()-relative.denominator().nbits()-1 min_accuracy = min(min_accuracy,int(lower_bits)) entries.add(**p,**data,digits=g.digits) assert not exactness(values) assert not names_its_rings(os.path.join(os.path.dirname(__file__),'generate.py')) records=entries.as_list() print('Local measurements:',json.dumps(measure(values))) assert 366 <= min_accuracy < 1000000 print('Minimum relative accuracy lower bound:',min_accuracy,'bits.') expected = {(b,n) for b in BOUNDARIES for n in range(1,MAX_N+1)} lookup = {(r['params']['boundary'],int(r['params']['n'])):r for r in records} assert set(lookup) == expected and len(records)==len(expected)==72 det_checked=0 brackets={} for key,r in lookup.items(): lo,hi=decimal_interval(str(r['number'])) brackets[key]=(lo,hi) if hi == 0: assert key in {('free-free',1),('free-free',2),('pinned-free',1)} continue dl=boundary_determinant(key[0],rb(lo),arb_trig) dh=boundary_determinant(key[0],rb(hi),arb_trig) assert dl.is_finite() and dh.is_finite() assert (dl<0 and dh>0) or (dl>0 and dh<0), (key,'stored interval lacks a determinant sign change') det_checked+=1 for b in BOUNDARIES: for n in range(1,MAX_N): if b=='free-free' and n==1: assert brackets[b,n]==brackets[b,n+1]==(0,0) else: assert brackets[b,n][1] < brackets[b,n+1][0],(b,n,'unordered') for n in range(3,MAX_N+1): assert lookup['free-free',n]['number']==lookup['clamped-clamped',n-2]['number'] assert lookup['free-free',n]['equals']=='clamped-clamped,%d'%(n-2) for n in range(2,MAX_N+1): assert lookup['pinned-free',n]['number']==lookup['clamped-pinned',n-1]['number'] assert lookup['pinned-free',n]['equals']=='clamped-pinned,%d'%(n-1) print('Certified boundary determinant sign changes:',det_checked, '; exact zeros: 3; ordering and 21 isospectral equalities passed.') # Rounding to a common 90-place grid is solely a comparison, never stored data. expected_integers={} for key,r in lookup.items(): q=rational(str(r['number'])) * 10**90 expected_integers[key]=ZZ((q+QQ(1)/2).floor()) max_error=mp.mpf(0) def independent(key): nonlocal max_error b,n=key if (b=='free-free' and n<=2) or (b=='pinned-free' and n==1): return ZZ(0) k=n-(2 if b=='free-free' else 1 if b=='pinned-free' else 0) with mp.workdps(200): if b=='clamped-free':guess=(mp.mpf(k)-mp.mpf('0.5'))*mp.pi elif b in ('clamped-clamped','free-free'):guess=(mp.mpf(k)+mp.mpf('0.5'))*mp.pi elif b in ('clamped-pinned','pinned-free'):guess=(mp.mpf(k)+mp.mpf('0.25'))*mp.pi else:guess=k*mp.pi roots=[] for precision in (150,200): with mp.workdps(precision): f=lambda x:boundary_determinant(b,x,mp_trig) x=mp.findroot(f,(guess-mp.mpf('0.1'),guess+mp.mpf('0.1')), solver='secant',tol=mp.mpf(10)**(-precision+30),maxsteps=80) roots.append(x) assert abs(roots[0]-roots[1]) < mp.mpf('1e-120'), key lo,hi=brackets[key] mpl=mp.mpf(str(lo.numerator()))/int(lo.denominator()) mph=mp.mpf(str(hi.numerator()))/int(hi.denominator()) assert mpl