_f64_igammaq_gate.nx source

1// _f64_igammaq_gate.nx -- ULP gate for integer-order Q(n,x) vs a sovereign 120-bit 2// finite-sum oracle (same closed form computed in bigfloat; bfc_exp reused from the erfc 3// organ). Grid n=1..12, x=k/2 (0.5..12). Plus the closed-form rail Q(1,x)==e^(-x) 4// (bit-exact: the n=1 kernel IS e^(-x)). RUNG bar: max ULP <= 4. 5// Durable: IGAMMAQ-GATE -> knowledge/status/math_engine.log. 6// license_tier: ORIGINAL 7import "nx_syscalls.nx" 8import "nx_bigfloat120.nx" 9import "nx_bigfloat120_div.nx" 10import "nx_bigfloat120_erfc.nx" 11import "nx_f64_exp.nx" 12import "_pe_f64igammaq.nx" 13 14func iq_p(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 } 15func iq_f(fd: i64, s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(fd,s,n); return 0 } 16func iq_n(fd: i64, v: i64) -> i64 { let bb: *u8=sys_mmap(28); var m: i64=v; if m<0{m=0-m; sys_write(fd,"-" as *u8,1)}; let t: *u8=sys_mmap(28); 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; while i<k{bb[i]=t[k-1-i];i=i+1}; sys_write(fd,bb,k); return 0 } 17func iq_ulp(a: i64, b: i64) -> i64 { if a >= b { return a - b } return b - a } 18 19// bigfloat Q(n,x) = e^(-x) * sum_{k=0}^{n-1} x^k/k! (the oracle) 20func iq_bf(out: *i64, n: i64, xb: *i64) -> i64 { 21 let term: *i64 = bf_new() 22 bf_set_int(term, 1) 23 let sum: *i64 = bf_new() 24 bf_copy(sum, term) 25 let t: *i64 = bf_new() 26 let snew: *i64 = bf_new() 27 var k: i64 = 1 28 while k < n { 29 bf_mul(t, term, xb) 30 bf_div_small(term, t, k) 31 bf_add(snew, sum, term) 32 bf_copy(sum, snew) 33 k = k + 1 34 } 35 let ex: *i64 = bf_new() 36 bfc_exp(ex, xb) 37 let one: *i64 = bf_new() 38 bf_set_int(one, 1) 39 let emx: *i64 = bf_new() 40 bf_div(emx, one, ex) 41 bf_mul(out, emx, sum) 42 return 0 43} 44 45func main() -> i64 { 46 iq_p("=== F64-IGAMMAQ GATE: integer-order Q(n,x) vs 120-bit finite-sum oracle ===\n" as *u8) 47 let kb: *i64 = bf_new() 48 let xb: *i64 = bf_new() 49 let oc: *i64 = bf_new() 50 var maxulp: i64 = 0 51 var worst_n: i64 = 0 52 var worst_x: i64 = 0 53 var npts: i64 = 0 54 var n: i64 = 1 55 while n <= 12 { 56 var kx: i64 = 1 57 while kx <= 24 { 58 bf_set_int(kb, kx) 59 bf_div_small(xb, kb, 2) // x = kx/2 60 let xbits: i64 = bf_to_f64(xb, 0, 0) 61 iq_bf(oc, n, xb) 62 let want: i64 = bf_to_f64(oc, 0, 0) 63 let got: i64 = nx_f64_igammaq(n, xbits) 64 let u: i64 = iq_ulp(got, want) 65 if u > maxulp { maxulp = u; worst_n = n; worst_x = xbits } 66 npts = npts + 1 67 kx = kx + 1 68 } 69 n = n + 1 70 } 71 // closed-form rail: Q(1,x) must equal e^(-x) bit-exactly at a few x 72 var rail_bad: i64 = 0 73 var rk: i64 = 1 74 while rk <= 10 { 75 bf_set_int(kb, rk) 76 bf_div_small(xb, kb, 2) 77 let xr: i64 = bf_to_f64(xb, 0, 0) 78 let q1: i64 = nx_f64_igammaq(1, xr) 79 let ex: i64 = nx_f64_exp(nx_f64_neg(xr)) 80 if q1 != ex { rail_bad = rail_bad + 1 } 81 rk = rk + 1 82 } 83 84 var ok: i64 = 0 85 if maxulp <= 4 { if rail_bad == 0 { ok = 1 } } 86 let lfd: i64 = sys_openat_append("knowledge/status/math_engine.log" as *u8, 0x1a4) 87 if lfd >= 0 { 88 iq_f(lfd, "IGAMMAQ-GATE epoch=" as *u8); iq_n(lfd, sys_now_realtime_sec()) 89 iq_f(lfd, " points=" as *u8); iq_n(lfd, npts) 90 iq_f(lfd, " max_ulp=" as *u8); iq_n(lfd, maxulp) 91 iq_f(lfd, " worst_n=" as *u8); iq_n(lfd, worst_n) 92 iq_f(lfd, " worst_x_bits=" as *u8); iq_n(lfd, worst_x) 93 iq_f(lfd, " q1_rail_bad=" as *u8); iq_n(lfd, rail_bad) 94 iq_f(lfd, " dlmf=ch8 form=closed-finite-sum oracle=120bit(SOVEREIGN)" as *u8) 95 if ok == 1 { iq_f(lfd, " verdict=GREEN\n" as *u8) } else { iq_f(lfd, " verdict=RED\n" as *u8) } 96 sys_close(lfd) 97 } 98 iq_p(" points=" as *u8); iq_n(1, npts) 99 iq_p(" max_ulp=" as *u8); iq_n(1, maxulp) 100 iq_p(" q1_rail_bad=" as *u8); iq_n(1, rail_bad) 101 if ok == 1 { iq_p(" IGAMMAQ-GATE: GREEN (integer-order incomplete gamma <=4 ulp + Q(1,x)=e^-x rail)\n" as *u8); sys_exit(0); return 0 } 102 iq_p(" IGAMMAQ-GATE: RED (measured gap named)\n" as *u8) 103 sys_exit(1) 104 return 1 105}