#!/usr/bin/env python3 """Exact algebra checks for OMEGA--Navier--Stokes geometric transport. Requires Python 3.9+ and SymPy (tested with 1.14.0). Run: python verify_omega_ns.py Acceptance uses explicit exceptions, including under python -O. This checks symbolic identities, not a PDE discretization, a Lean proof, or the external Navier--Stokes blowup manuscript. Smooth localization and global existence for the explicitly prescribed forced fields are proved in the accompanying PDF. Version 2.0 adds initial-value and localization scope checks. Version 3.0 adds cutoff vorticity, energy, nonuniqueness, perturbation and tetrahedral checks; the frozen source packages are not modified. No numerical solution is claimed. """ import sympy as s def require(condition, name): if not condition: raise ValueError(name) def zero(expr, name): entries = list(expr) if isinstance(expr, s.MatrixBase) else [expr] require(all(s.simplify(v) == 0 for v in entries), name) def audit(): names = [] def check(expr, name): zero(expr, name) names.append(name) x, y, z = s.symbols('x y z', real=True) gamma, beta, dg, db = s.symbols('gamma beta dg db', real=True) X = s.Matrix([x, y, z]) H = s.Matrix([[1, s.Rational(1, 2)], [s.Rational(1, 2), -1]]) M = s.Matrix([[0, 1], [1, 1]]) K = (2*M-s.eye(2))/s.sqrt(5) phi = (1+s.sqrt(5))/2 phibar = (1-s.sqrt(5))/2 tau = s.symbols('tau', positive=True) A = tau*M check(K*K-s.eye(2), 'K squared equals identity') check(K.T-K, 'K is symmetric') check(s.trace(K), 'K is traceless') check(K.T*H+H*K, 'K preserves the golden quadratic form infinitesimally') check(M.T*H*M+H, 'Fibonacci matrix reverses the golden norm') check(A.T*H*A+tau**2*H, 'normalized matrix norm identity') ep = (s.eye(2)+K)/2 em = (s.eye(2)-K)/2 check(ep*em, 'spectral projectors are orthogonal') check(ep*ep-ep, 'positive spectral projector is idempotent') check(em*em-em, 'negative spectral projector is idempotent') lp = s.Matrix([[1, phi]]) lm = s.Matrix([[1, phibar]]) check(lp*K-lp, 'positive golden observation has eigenvalue one') check(lm*K+lm, 'conjugate observation has eigenvalue minus one') eta, theta = s.symbols('eta theta', real=True) c, d = s.symbols('c d', positive=True) def E(a): return s.exp(a)*ep+s.exp(-a)*em def flow(scale, shape): ans = s.zeros(3) ans[:2, :2] = scale*E(shape) ans[2, 2] = scale**-2 return ans P = flow(c, eta) check(E(eta).T*H*E(eta)-H, 'golden hyperbolic flow preserves norm') check(P.det()-1, 'three dimensional flow preserves volume') check(flow(c, eta)*flow(d, theta)-flow(c*d, eta+theta), 'two parameter group law') check(lp*P[:2, :2]-c*s.exp(eta)*lp, 'scalar observation of the deformation') N = x*x+x*y-y*y J = z*N transported = P*X check(J.subs(dict(zip((x, y, z), transported)), simultaneous=True)-J, 'cubic material invariant under the finite deformation') B = s.zeros(3) B[:2, :2] = gamma*s.eye(2)+beta*K B[2, 2] = -2*gamma Bd = B.diff(gamma)*dg+B.diff(beta)*db velocity = B*X pressure = -(X.T*(Bd+B*B)*X)[0]/2 gradient = lambda f: s.Matrix([s.diff(f, xi) for xi in X]) divergence = lambda v: sum(s.diff(v[i], X[i]) for i in range(3)) laplace = lambda v: s.Matrix([sum(s.diff(vi, xi, 2) for xi in X) for vi in v]) def curl(v): return s.Matrix([s.diff(v[2], y)-s.diff(v[1], z), s.diff(v[0], z)-s.diff(v[2], x), s.diff(v[1], x)-s.diff(v[0], y)]) nu = s.symbols('nu', positive=True) check(divergence(velocity), 'affine incompressibility') check(laplace(velocity), 'affine viscous Laplacian vanishes') check(Bd*X+velocity.jacobian(X)*velocity-nu*laplace(velocity)+gradient(pressure), 'full affine Navier Stokes momentum residual is zero') check(sum(s.diff(J, xi)*velocity[i] for i, xi in enumerate(X)), 'material derivative of cubic invariant vanishes') check(c*gamma*P.diff(c)+beta*P.diff(eta)-B*P, 'finite deformation solves particle flow equation') potential = velocity.cross(X)/3 check(curl(potential)-velocity, 'vector potential recovers affine velocity') chi = s.Function('chi')(x, y, z) localized = curl(chi*potential) check(divergence(localized), 'arbitrary smooth spatial cutoff preserves divergence') check(localized-curl(potential)*chi-gradient(chi).cross(potential), 'localized velocity product rule') check(phi**2*ep+phi**-2*em-M*M, 'continuous golden flow equals M squared') check(tau**2*(phi**2*ep+phi**-2*em)-A*A, 'normalized golden flow equals A squared') omega, dw = s.symbols('omega dw', real=True) Bsw = s.Matrix([[gamma, -omega, 0], [omega, gamma, 0], [0, 0, -2*gamma]]) Csw = Bsw.diff(gamma)*dg+Bsw.diff(omega)*dw+Bsw*Bsw sym = (Csw+Csw.T)/2 usw = Bsw*X psw = -(X.T*sym*X)[0]/2 expected_force = (dw+2*gamma*omega)*s.Matrix([-y, x, 0]) check(Csw*X+gradient(psw)-expected_force, 'swirling branch force identity') check(curl(usw)-s.Matrix([0, 0, 2*omega]), 'swirling branch vorticity') v1, v2, w1, w2 = s.symbols('v1 v2 w1 w2', real=True) V, W = s.Matrix([v1, v2]), s.Matrix([w1, w2]) qform = lambda q: (q.T*H*q)[0] check((qform(V+W)+qform(V-W))/2-qform(V)-s.trace(H*(W*W.T)), 'golden quadratic mean defect equals stress contraction') check(H*H-s.Rational(5, 4)*s.eye(2), 'golden stress bound spectral constant') vertices = [] for a in (-1, 1): for b in (-1, 1): vertices.extend((s.Matrix([0, a, b*phi]), s.Matrix([a, b*phi, 0]), s.Matrix([b*phi, 0, a]))) radius2 = 1+phi**2 check(sum(vertices, s.zeros(3, 1)), 'icosahedral first moment') second = sum((v*v.T/radius2 for v in vertices), s.zeros(3))/12 check(second-s.eye(3)/3, 'icosahedral second moment') discrete_energy = sum((v.T*B.T*B*v)[0]/radius2 for v in vertices)/12 check(discrete_energy-s.trace(B.T*B)/3, 'exact cubature of affine squared speed') Q, S = 339949771344778, 10**15 require((Q*(Q*26//S))//S == 2 and Q*Q*26//(S*S) == 3, 'iterated quantization must remain distinct from one final floor') names.append('exact counterexample separating iterated quantization and final rounding') c0 = s.symbols('c0', positive=True) eta0 = s.symbols('eta0', real=True) relative = flow(c/c0, eta-eta0) check(relative*flow(c0, eta0)-P, 'relative initial-to-current deformation') check(relative.subs({c:c0, eta:eta0})-s.eye(3), 'relative deformation starts at identity') check(c*gamma*relative.diff(c)+beta*relative.diff(eta)-B*relative, 'relative deformation solves particle flow equation') check(lp*relative[:2,:2]-(c/c0)*s.exp(eta-eta0)*lp, 'observation uses relative scale and shape') check(relative[:2,:2].T*H*relative[:2,:2]-(c/c0)**2*H, 'norm law uses relative initial scale') lam = phi/(phi+s.pi) require((1-lam**2).is_positive is True, 'nonzero constant schedule refutes unnormalized initial-value formula') names.append('constant-schedule counterexample to the old unnormalized formula') a_ref = s.symbols('a_ref', positive=True) observed = 1+(lp*s.Matrix([x,y]))[0]/a_ref check(1+(lp*(c*s.Matrix([x,y])))[0]/a_ref-(1+c*(observed-1)), 'dimensionless scalar observation of planar similarity') radial = s.Function('h')(x*x+y*y+z*z) radial_velocity = curl(radial*potential) hprime = s.Subs(s.Derivative(s.Function('h')(s.Symbol('v')),s.Symbol('v')), s.Symbol('v'),x*x+y*y+z*z) check(gradient(J).dot(radial_velocity)+2*hprime*(X.dot(velocity))*J, 'cutoff transition material-invariant defect') rr, rt, tt = s.symbols('rr rt tt', real=True) covariance = s.Matrix([[rr,rt],[rt,tt]]) check(s.trace(H*A*covariance*A.T)+tau**2*s.trace(H*covariance), 'general covariance transformation under normalized Fibonacci matrix') discrete_stress = sum((B*v*v.T*B.T/radius2 for v in vertices),s.zeros(3))/12 check(discrete_stress-B*B.T/3, 'exact cubature of affine stress tensor') ell0, r, dr, dt_angle, dz = s.symbols('ell0 r dr dt_angle dz', positive=True) physical_metric = dr**2+r**2*(dt_angle-2*s.pi*dr/ell0)**2+dz**2 check(s.expand(physical_metric)-(dr**2+r**2*dt_angle**2 -4*s.pi*r**2*dt_angle*dr/ell0+4*s.pi**2*r**2*dr**2/ell0**2+dz**2), 'dimensionally scaled spiral pullback metric') # Version 3.0: identities added after the adversarial review. radius_squared=X.dot(X) qrad=s.symbols('qrad') hp=s.Subs(s.diff(s.Function('h')(qrad),qrad),qrad,radius_squared) hpp=s.Subs(s.diff(s.Function('h')(qrad),qrad,2),qrad,radius_squared) expected_radial=(radial+s.Rational(2,3)*radius_squared*hp)*velocity-s.Rational(2,3)*hp*X.dot(velocity)*X check(radial_velocity-expected_radial,'explicit radial localized velocity') check(curl(radial_velocity)-(14*hp+4*radius_squared*hpp)/3*X.cross(velocity), 'explicit vorticity generated by the cutoff') check(s.trace(B.T*B)-6*gamma**2-2*beta**2,'nonzero core strain density') check(sum(s.diff(velocity.dot(velocity)/2,xi,2) for xi in X)-s.trace(B.T*B), 'local energy Laplacian balances gradient density') check(gradient(J).dot(localized)+J*velocity.dot(gradient(chi)), 'general spatial-cutoff cubic invariant defect') perturb=s.symbols('perturb',real=True) check(gradient(J).dot((B+perturb*s.diag(1,-1,0))*X)-2*perturb*z*(x*x+y*y), 'arbitrarily small trace-free perturbation breaks cubic invariant') check(J.subs(z,0),'cubic invariant is identically zero on marked plane') isotropic=s.symbols('isotropic',nonnegative=True) check(s.trace(H*(isotropic*s.eye(2))),'golden stress observation misses isotropic energy') tetra=[s.Matrix(v)/s.sqrt(3) for v in [(1,1,1),(1,-1,-1),(-1,1,-1),(-1,-1,1)]] check(sum(tetra,s.zeros(3,1)),'tetrahedral first moment') check(sum((v*v.T for v in tetra),s.zeros(3))/4-s.eye(3)/3,'tetrahedral second moment') kap=s.symbols('kap',real=True) Mk=s.Matrix([[0,1],[1,kap]]) Hk=s.Matrix([[1,kap/2],[kap/2,-1]]) dk=s.sqrt(kap*kap+4) Kk=(2*Mk-kap*s.eye(2))/dk ak=(kap+dk)/2;am=(kap-dk)/2 check(Kk.T-Kk,'general real-parameter generator symmetry') check(s.trace(Kk),'general real-parameter generator trace') check(Kk*Kk-s.eye(2),'general real-parameter generator involution') check(Kk.T*Hk+Hk*Kk,'general real-parameter quadratic form preservation') check(Mk.T*Hk*Mk+Hk,'general real-parameter recurrence reverses its form') check(ak*am+1,'general real-parameter reciprocal eigenvalue relation') check(Mk*Mk-ak**2*(s.eye(2)+Kk)/2-am**2*(s.eye(2)-Kk)/2, 'general real-parameter two-step spectral identity') print('OMEGA_NS_GEOMETRIC_EXTENSION/3.0') print('SymPy', s.__version__) for name in names: print('PASS:', name) print('TOTAL:', len(names), 'exact checks') print('SCOPE: symbolic identities; smooth localization is proved in the paper; no new Lean build or Julia run.') # Independent integer/rational certificate, retained with its own scope. from fractions import Fraction from math import isqrt Q = 339949771344778 S = 10**15 DSTAR = 42189065 M = 10**30 LAMBDA_LO = 339949771344778929576702809321 LAMBDA_HI = 339949771344778929576702809322 def require(condition, message): if not condition: raise ValueError(message) def atan_bounds(reciprocal, terms=64): require(type(reciprocal) is int and reciprocal > 1, "atan reciprocal") require(type(terms) is int and terms > 0, "atan terms") partial = sum( (Fraction((-1)**k, (2*k + 1) * reciprocal**(2*k + 1)) for k in range(terms)), Fraction(0), ) adjacent = partial + Fraction( (-1)**terms, (2*terms + 1) * reciprocal**(2*terms + 1) ) return min(partial, adjacent), max(partial, adjacent) def lambda_bounds(): a_lo, a_hi = atan_bounds(5) b_lo, b_hi = atan_bounds(239) pi_lo, pi_hi = 16*a_lo - 4*b_hi, 16*a_hi - 4*b_lo require(3 < pi_lo < pi_hi < 4, "pi bounds") scale = 10**80 root = isqrt(5*scale*scale) require(root*root < 5*scale*scale < (root+1)**2, "sqrt(5) bounds") phi_lo = (1 + Fraction(root, scale))/2 phi_hi = (1 + Fraction(root+1, scale))/2 return phi_lo/(phi_lo+pi_hi), phi_hi/(phi_hi+pi_lo) def floor_sum(n, modulus, coefficient, offset=0): require(all(type(x) is int for x in (n, modulus, coefficient, offset)), "floor_sum requires integers") require(n >= 0 and modulus > 0 and coefficient >= 0 and offset >= 0, "floor_sum domain") total = 0 while True: total += (n-1)*n*(coefficient//modulus)//2 total += n*(offset//modulus) coefficient %= modulus offset %= modulus height = coefficient*n + offset if height < modulus: return total n, offset = divmod(height, modulus) modulus, coefficient = coefficient, modulus def F(value): require(type(value) is int and value >= 0, "nonnegative integer required") return Q*value//S def tick(symbol, steps): require(type(symbol) is int and symbol >= 0, "tick symbol") require(type(steps) is int and steps >= 0, "tick steps") value = (symbol+2)*S for _ in range(steps): value = F(value) return value def orbit_certificate(seed): require(type(seed) is int and seed >= 0, "orbit seed") orbit, remainders = [seed], [] while orbit[-1]: next_value, rho = divmod(Q*orbit[-1], S) require(next_value < orbit[-1], "orbit must decrease") require(0 <= rho < S, "orbit remainder bound") require(Q*orbit[-1] == S*next_value + rho, "orbit identity") orbit.append(next_value) remainders.append(rho) recovered = orbit[-1] for rho in reversed(remainders): recovered, nonzero = divmod(S*recovered + rho, Q) require(nonzero == 0, "reverse divisibility") require(recovered == seed, "reverse seed") return orbit, remainders def arithmetic_audit(): lo, hi = lambda_bounds() require(Fraction(Q, S) < Fraction(LAMBDA_LO, M) < lo < hi < Fraction(LAMBDA_HI, M), "lambda rational enclosure") rational_sum = floor_sum(DSTAR, S, Q) upper_sum = floor_sum(DSTAR, M, LAMBDA_HI) require(rational_sum == upper_sum == 302541195258184, "agreement for every 0 <= d < dstar") require(Q*DSTAR//S == 14342162, "rational dstar floor") require(LAMBDA_LO*DSTAR//M == LAMBDA_HI*DSTAR//M == 14342163, "real dstar floor") require(Q*DSTAR % S == 999999976452570, "dstar remainder") print("PASS: first nonnegative floor-map disagreement is", DSTAR) print("Exact common floor sum on 0..42189064:", rational_sum) print("lambda lies strictly between", LAMBDA_LO, "and", LAMBDA_HI, "divided by", M) for seed, expected in ( (32767, (11139, 134157654340726)), (32768, (11139, 474107425685504)), ): require(divmod(Q*seed, S) == expected, "shared cell certificate") orbit_certificate(seed) orbit, _ = orbit_certificate(DSTAR) require(len(orbit)-1 == 17, "dstar orbit length") print("PASS: shared next cell, distinct remainders, exact reconstruction") print("dstar orbit:", orbit) gap = S for _ in range(31): gap = F(gap) require(gap == 2, "global step-31 gap bound") require(tick(30, 31) == 95 and tick(31, 31) == 97, "collision precursors") require(tick(30, 32) == tick(31, 32) == 32, "global step-32 collision") expected_rows = { 23: (256, 33386, 50079, 4290149), 24: (256, 11349, 17024, 1458435), 31: (256, 5, 8, 764), 32: (252, 1, 2, 259), } for steps, expected in expected_rows.items(): table = [tick(b, steps) for b in range(256)] require((len(set(table)), table[0], table[1], table[-1]) == expected, "byte tick table") table = [tick(b, 32) for b in range(256)] collisions = [(a, b, table[a]) for a in range(256) for b in range(a+1, 256) if table[a] == table[b]] require(collisions == [(30, 31, 32), (112, 113, 115), (171, 172, 175), (253, 254, 258)], "all step-32 byte collisions") print("PASS: all nonnegative integer symbols injective through step 31") print("Step-32 byte collisions (symbol, symbol, tick):", collisions) data = b"OMEGA-NAUGHT/1.0\n" + bytes([0, 1, 127, 128, 255]) + b"END" require(len(data) == 25, "object length") y_be = int.from_bytes(data, "big") require(y_be == 497785707148424555849356544956251674240246940257526445592132, "object Y_be") require(y_be % 32768 == 20036, "object D15 seed") orbit, _ = orbit_certificate(y_be % 32768) require(orbit == [20036, 6811, 2315, 786, 267, 90, 30, 10, 3, 1, 0], "object residual orbit") for source in (b"", b"\x00", b"\x00\x00", b"OMEGA", b"OMEGA\x00\xffNAUGHT", data): y_le = int.from_bytes(source, "little") value, residues = y_le, [] for _ in range(len(source)): next_value, rho = divmod(value, 256) require(value == 256*next_value + rho and 0 <= rho < 256, "byte tower identity") value = next_value residues.append(rho) require(value == 0 and len(residues) == len(source), "tower apex/length") require(bytes(residues) == source, "tower byte reconstruction") require(sum(r*256**i for i, r in enumerate(residues)) == y_le, "tower telescoping identity") reconstructed_be = sum(r*256**(len(residues)-1-i) for i, r in enumerate(residues)) require(reconstructed_be == int.from_bytes(source, "big"), "separate big-endian readout") print("PASS: supplied 25-byte object, residual orbit, and byte tower identities") print("SCOPE: the external SRSD decoder implementation and its own mutation suite/run_suite were not audited.") if __name__ == "__main__": audit() arithmetic_audit()