nx_dpendulum_gate.nx source

1// nx_dpendulum_gate.nx -- the DOUBLE PENDULUM, the continuous chaotic system named as the sim 2// workstream's flagship "no-float determinism under chaos" vehicle (companion to the discrete 3// nx_chaos_logistic_gate). Pure integer/no-float fixed point (all quantities x1e6 = "q"); reuses the 4// Taylor sin/cos from nx_pendulum. Equal masses + lengths (m1=m2=1, l1=l2=1, g=9.8); full nonlinear 5// COUPLED equations of motion (myphysicslab form), denom = 3 - cos(2(th1-th2)) in [2,4] => no 6// singularity. Semi-implicit (symplectic-style) Euler: update omega from current theta, then theta 7// from the new omega -- energy-bounded enough to integrate the chaotic regime over a short horizon. 8// 9// THREE measured properties + liar-kill (GREEN means it discriminates): 10// (1) SENSITIVE DEPENDENCE: at LARGE amplitude (th1=th2=2.5 rad) two runs whose th2 differ by 1 11// micro-rad DIVERGE to O(1 rad) -- deterministic chaos. 12// (2) DETERMINISM: the SAME integer run recomputed gives a BIT-IDENTICAL max-separation (no-float 13// => exactly reproducible even under chaos = the Reproducibility-EXCEEDS property). 14// (3) NEGATIVE CONTROL: the SAME integrator at SMALL amplitude (th1=th2=0.05 rad) is quasi-linear 15// (two normal modes, NOT chaotic) -> the 1-micro-rad perturbation stays bounded/tiny. So the 16// gate cannot be fooled into calling regular motion "chaos". GREEN needs all three. 17// Sovereign: imports only nx_syscalls. Additive. license_tier: ORIGINAL 18import "nx_syscalls.nx" 19 20const SC: i64 = 1000000 // x1e6 fixed-point scale 21const GQ: i64 = 9800000 // g = 9.8 in q 22const DTDIV: i64 = 1000 // dt = 1ms: x += rate/1000 23 24func dp_len(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } return n } 25func dp_p(s: *u8) -> i64 { let n: i64 = dp_len(s); sys_write(1, s, n); return 0 } 26func dp_pn(v: i64) -> i64 { 27 let bb: *u8 = sys_mmap(28); var m: i64 = v 28 if m < 0 { sys_write(1, "-" as *u8, 1); m = 0 - m } 29 let t: *u8 = sys_mmap(28); var k: i64 = 0 30 if m == 0 { t[0] = 48 as u8; k = 1 } 31 while m > 0 { t[k] = (48 + (m % 10)) as u8; m = m / 10; k = k + 1 } 32 var i: i64 = 0; while i < k { bb[i] = t[k - 1 - i]; i = i + 1 } 33 sys_write(1, bb, k); return 0 34} 35 36// Taylor sine; input micro-rad, output sin x1e6 (the proven nx_pendulum sin_fixed). 37func sin_fixed(Xin: i64) -> i64 { 38 let PI: i64 = 3141593; let TWO_PI: i64 = 6283185; let HALF: i64 = 1570796 39 var sign: i64 = 1; var X: i64 = Xin 40 if X < 0 { X = 0 - X; sign = 0 - 1 } 41 X = X % TWO_PI 42 if X > PI { X = TWO_PI - X; sign = 0 - sign } 43 if X > HALF { X = PI - X } 44 let x2: i64 = X * X / SC 45 var term: i64 = X; var sum: i64 = X; var k: i64 = 1 46 while k <= 8 { term = (0 - term) * x2 / SC / ((2 * k) * (2 * k + 1)); sum = sum + term; k = k + 1 } 47 return sum * sign 48} 49func cos_fixed(X: i64) -> i64 { return sin_fixed(X + 1570796) } 50 51// fixed-point multiply of two q values -> q. 52func mulq(a: i64, b: i64) -> i64 { return a * b / SC } 53 54// one semi-implicit step of the double pendulum. s[0..3] = th1,th2,om1,om2 (all q). 55func dp_step(s: *i64) -> i64 { 56 let th1: i64 = s[0]; let th2: i64 = s[1]; let om1: i64 = s[2]; let om2: i64 = s[3] 57 let d: i64 = th1 - th2 58 let sind: i64 = sin_fixed(d) 59 let cosd: i64 = cos_fixed(d) 60 let cos2d: i64 = cos_fixed(2 * d) 61 let denom: i64 = 3 * SC - cos2d 62 let w1s: i64 = mulq(om1, om1) 63 let w2s: i64 = mulq(om2, om2) 64 let s1: i64 = sin_fixed(th1) 65 let s12: i64 = sin_fixed(th1 - 2 * th2) 66 let c1: i64 = cos_fixed(th1) 67 // alpha1 = [-3g sin th1 - g sin(th1-2th2) - 2 sin d (om2^2 + om1^2 cos d)] / denom 68 let inner1: i64 = w2s + mulq(w1s, cosd) 69 let num1: i64 = (0 - 3) * mulq(GQ, s1) - mulq(GQ, s12) - 2 * mulq(sind, inner1) 70 let a1: i64 = num1 * SC / denom 71 // alpha2 = [2 sin d (2 om1^2 + 2g cos th1 + om2^2 cos d)] / denom 72 let inner2: i64 = 2 * w1s + 2 * mulq(GQ, c1) + mulq(w2s, cosd) 73 let num2: i64 = 2 * mulq(sind, inner2) 74 let a2: i64 = num2 * SC / denom 75 // semi-implicit (symplectic-style): omega first, then theta with NEW omega 76 let no1: i64 = om1 + a1 / DTDIV 77 let no2: i64 = om2 + a2 / DTDIV 78 s[2] = no1; s[3] = no2 79 s[0] = th1 + no1 / DTDIV 80 s[1] = th2 + no2 / DTDIV 81 return 0 82} 83 84// run two pendulums from the same start except B's th2 is perturbed by `pert` micro-rad; both from 85// rest. Return MAX (|d th1| + |d th2|) over all N steps. Deterministic -> reproducible. 86func dp_run_sep(th0: i64, pert: i64, n: i64) -> i64 { 87 let A: *i64 = sys_mmap(8 * 4) as *i64 88 let B: *i64 = sys_mmap(8 * 4) as *i64 89 A[0] = th0; A[1] = th0; A[2] = 0; A[3] = 0 90 B[0] = th0; B[1] = th0 + pert; B[2] = 0; B[3] = 0 91 var mx: i64 = 0 92 var i: i64 = 0 93 while i < n { 94 dp_step(A); dp_step(B) 95 var d1: i64 = A[0] - B[0]; if d1 < 0 { d1 = 0 - d1 } 96 var d2: i64 = A[1] - B[1]; if d2 < 0 { d2 = 0 - d2 } 97 let sep: i64 = d1 + d2 98 if sep > mx { mx = sep } 99 i = i + 1 100 } 101 return mx 102} 103 104func main() -> i64 { 105 dp_p("=== nx_dpendulum_gate: no-float determinism UNDER CHAOS (double pendulum) ===\n" as *u8) 106 let BIG: i64 = 100000 // chaos threshold: > 0.1 rad of divergence 107 let SMALL: i64 = 10000 // control threshold: < 0.01 rad 108 let N: i64 = 14000 // 14 s at dt=1ms -- enough for the (observed) Lyapunov growth to saturate 109 let LARGE_TH: i64 = 2500000 // 2.5 rad (chaotic) 110 let SMALL_TH: i64 = 50000 // 0.05 rad (quasi-linear, NOT chaotic) 111 112 let chaos1: i64 = dp_run_sep(LARGE_TH, 1, N) 113 let chaos2: i64 = dp_run_sep(LARGE_TH, 1, N) // recompute -> must be identical 114 let order1: i64 = dp_run_sep(SMALL_TH, 1, N) 115 116 dp_p(" m1=m2=1 l1=l2=1 g=9.8 dt=1ms N=" as *u8); dp_pn(N) 117 dp_p("; perturbation = 1 micro-rad in th2\n" as *u8) 118 dp_p(" (1) sensitive-dependence th=2.5rad max sep = " as *u8); dp_pn(chaos1) 119 dp_p(" (need > " as *u8); dp_pn(BIG); dp_p(")\n" as *u8) 120 dp_p(" (2) determinism re-run th=2.5rad max sep = " as *u8); dp_pn(chaos2) 121 dp_p(" (must == run 1)\n" as *u8) 122 dp_p(" (3) neg-control (order) th=0.05rad max sep = " as *u8); dp_pn(order1) 123 dp_p(" (must < " as *u8); dp_pn(SMALL); dp_p(")\n" as *u8) 124 125 var pass: i64 = 0 126 var ok: i64 = 1 127 if chaos1 > BIG { pass = pass + 1 } else { ok = 0 } 128 if chaos1 == chaos2 { pass = pass + 1 } else { ok = 0 } 129 if order1 < SMALL { pass = pass + 1 } else { ok = 0 } 130 131 dp_p(" RESULT pass=" as *u8); dp_pn(pass); dp_p("/3 verdict=" as *u8) 132 if ok == 1 { 133 dp_p("GREEN (chaotic divergence, bit-reproducible no-float, regular motion stays regular)\n" as *u8) 134 sys_exit(0); return 0 135 } 136 dp_p("RED\n" as *u8) 137 sys_exit(1); return 1 138}