nx_sim_verification_gate.nx source

1// nx_sim_verification_gate.nx -- climb the credibility ladder's VERIFICATION axis with a real SOLUTION-VERIFICATION 2// artifact required by NASA-STD-7009 / ASME V&V 10 & 20 / AIAA G-077: an OBSERVED ORDER-OF-ACCURACY convergence 3// study. Integrate the harmonic oscillator (a=-x, omega=1, exact x(t)=A cos t) to a fixed time T=1 with shrinking 4// step dt=1/D, and measure the error vs the EXACT analytic solution (our sovereign cos). For a 2nd-order method 5// (leapfrog) the error must drop ~4x each time dt halves (order p=log2(ratio)~2). The test also genuinely MEASURES 6// order: forward-Euler (1st order) drops only ~2x -- so a method falsely claimed "2nd order" is CAUGHT. This 7// converts Verification from BASIC(2) toward CREDIBLE(3). No float, deterministic. GREEN iff 6/6. license_tier: ORIGINAL 8import "nx_syscalls.nx" 9 10func g_w(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 } 11func g_n(v: i64) -> i64 { var m: i64=v; if m<0{g_w("-");m=0-m} let t:*u8=sys_mmap(24); var k:i64=0; if m==0{t[0]=48 as u8;k=1}; while m>0{t[k]=(48+(m%10)) as u8;m=m/10;k=k+1}; var i:i64=0; let o:*u8=sys_mmap(24); while i<k{o[i]=t[k-1-i];i=i+1}; sys_write(1,o,k); return 0 } 12func g_row(id: *u8, ok: i64, pass: *i64) -> i64 { g_w(" "); g_w(id); g_w(": "); if ok==1 { g_w("OK\n"); pass[0]=pass[0]+1 } else { g_w("FAIL\n") } return 0 } 13func iabs(x: i64) -> i64 { if x<0 { return 0-x } return x } 14// sovereign cos for the EXACT reference solution (Taylor + range reduction, SC=1e6). 15func sin_fixed(Xin: i64, SC: i64) -> i64 { 16 let PI: i64=3141593; let TWO_PI: i64=6283185; let HALF: i64=1570796 17 var sign: i64=1; var X: i64=Xin 18 if X<0 { X=0-X; sign=0-1 } 19 X = X % TWO_PI 20 if X > PI { X = TWO_PI - X; sign = 0-sign } 21 if X > HALF { X = PI - X } 22 let x2: i64 = X*X/SC 23 var term: i64=X; var sum: i64=X; var k: i64=1 24 while k<=8 { term = (0-term)*x2/SC/((2*k)*(2*k+1)); sum=sum+term; k=k+1 } 25 return sum*sign 26} 27func cos_fixed(X: i64, SC: i64) -> i64 { return sin_fixed(X+1570796, SC) } 28 29// leapfrog (a=-x): return x after N steps with dt=1/D. 30func lf_xT(A: i64, D: i64, N: i64) -> i64 { 31 var x: i64=A; var v: i64=0; var i: i64=0 32 while i<N { let vh: i64=v - x/(2*D); x = x + vh/D; v = vh - x/(2*D); i=i+1 } 33 return x 34} 35// forward-Euler (a=-x): return x after N steps. 36func eu_xT(A: i64, D: i64, N: i64) -> i64 { 37 var x: i64=A; var v: i64=0; var i: i64=0 38 while i<N { let nx: i64=x + v/D; let nv: i64=v - x/D; x=nx; v=nv; i=i+1 } 39 return x 40} 41 42func main() -> i64 { 43 let pass: *i64 = sys_mmap(8) as *i64; pass[0]=0 44 g_w("=== NX-SIM-VERIFICATION GATE (observed order-of-accuracy: solution verification per NASA-7009 / ASME V&V) ===\n") 45 let SC: i64=1000000; let A: i64=1000000000 // positions scaled to 1e9 so dt^2 error >> fixed-point rounding floor 46 // EXACT reference at T=1: x(1) = A*cos(1) = 1000 * (SC*cos(1)) since A = 1000*SC. 47 let xref: i64 = 1000*cos_fixed(1000000, SC) 48 g_w(" exact x(1)=A*cos(1)="); g_n(xref); g_w(" (cos(0)*1000="); g_n(1000*cos_fixed(0,SC)); g_w(")\n") 49 50 // errors at T=1 (N=D steps so time = N/D = 1) for refining dt=1/D. 51 let Ds: *i64=sys_mmap(8*8) as *i64; Ds[0]=16; Ds[1]=32; Ds[2]=64; Ds[3]=128 52 let lpe: *i64=sys_mmap(8*8) as *i64; let eue: *i64=sys_mmap(8*8) as *i64 53 var i: i64=0 54 while i<4 { let D: i64=Ds[i]; lpe[i]=iabs(lf_xT(A,D,D)-xref); eue[i]=iabs(eu_xT(A,D,D)-xref); i=i+1 } 55 56 g_w(" leapfrog err: D16="); g_n(lpe[0]); g_w(" D32="); g_n(lpe[1]); g_w(" D64="); g_n(lpe[2]); g_w(" D128="); g_n(lpe[3]); g_w("\n") 57 g_w(" euler err: D16="); g_n(eue[0]); g_w(" D32="); g_n(eue[1]); g_w(" D64="); g_n(eue[2]); g_w(" D128="); g_n(eue[3]); g_w("\n") 58 // ratios *100 (err(D)/err(2D)); ~400 => order 2, ~200 => order 1 59 let lr0: i64=lpe[0]*100/lpe[1]; let lr1: i64=lpe[1]*100/lpe[2]; let lr2: i64=lpe[2]*100/lpe[3] 60 let er0: i64=eue[0]*100/eue[1]; let er1: i64=eue[1]*100/eue[2]; let er2: i64=eue[2]*100/eue[3] 61 g_w(" leapfrog ratio*100 (->~400=order2): "); g_n(lr0); g_w(" "); g_n(lr1); g_w(" "); g_n(lr2); g_w("\n") 62 g_w(" euler ratio*100 (->~200=order1): "); g_n(er0); g_w(" "); g_n(er1); g_w(" "); g_n(er2); g_w("\n") 63 64 g_row("EXACT REFERENCE: sovereign cos gives the analytic solution (cos(0)=A, x(1)=A cos1 ~0.5403A)" as *u8, ((iabs(cos_fixed(0,SC)-SC)<1500) as i64)*((iabs(xref-540303000)<2000000) as i64), pass) 65 var conv: i64=0; if lpe[1]<lpe[0] { if lpe[2]<lpe[1] { if lpe[3]<lpe[2] { conv=1 } } } 66 g_row("CONVERGENCE (leapfrog): error decreases monotonically under dt refinement (the sim converges)" as *u8, conv, pass) 67 var ord2: i64=0; if lr0>300 { if lr0<500 { if lr1>300 { if lr1<500 { ord2=1 } } } } 68 g_row("ORDER-OF-ACCURACY (leapfrog): err ratio ~4x per dt-halving => observed 2nd order (matches theory)" as *u8, ord2, pass) 69 var ord1: i64=0; if er0>150 { if er0<280 { if er1>150 { if er1<280 { ord1=1 } } } } 70 g_row("ORDER-OF-ACCURACY (euler): err ratio ~2x => observed 1st order (matches theory)" as *u8, ord1, pass) 71 g_row("DISCRIMINATION: the order test SEPARATES 2nd-order from 1st-order (leapfrog ratio > euler ratio)" as *u8, (lr1 > er1*13/10) as i64, pass) 72 var liar: i64=0; if ord2==1 { if er1<280 { liar=1 } } 73 g_row("LIAR-KILL: a method falsely claimed 2nd-order (Euler) is CAUGHT -- its ratio is ~2, not ~4" as *u8, liar, pass) 74 75 g_w("SIM-VERIFICATION-GATE rows=6 pass="); g_n(pass[0]) 76 if pass[0]==6 { g_w(" verdict=GREEN\n"); sys_exit(0); return 0 } 77 g_w(" verdict=RED\n"); sys_exit(1); return 1 78}