code wiki / _hdl_build / _f64_gamma_refl_neg_gate.nx

_f64_gamma_refl_neg_gate.nx source

↩ module page · 56 lines · 3064 B

1// _f64_gamma_refl_neg_gate.nx -- proves the SAME reflection kernel (_pe_f64gamma_refl) 2// is correct for NEGATIVE x, so gamma is full-domain across the reals minus the 3// non-positive-integer poles. Grid x = -(2j+1)/4 (quarter-points -0.25..-9.25, far from 4// the integer poles where sin(pi*x)=0). Oracle = the full-domain bf_gamma_f64 (Stirling 5// +shift+reflection, 120-bit). Bar <= 16 ulp (composes Gamma(1-x) <=7ulp + sin + div). 6// Durable: GAMMAREFLNEG-GATE -> knowledge/status/math_engine.log. 7// license_tier: ORIGINAL 8import "nx_syscalls.nx" 9import "nx_bigfloat120.nx" 10import "nx_bigfloat120_div.nx" 11import "nx_bigfloat120_gamma.nx" 12import "_pe_f64gamma_refl.nx" 13 14func gn_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 gn_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 gn_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 gn_ulp(a: i64, b: i64) -> i64 { if a >= b { return a - b } return b - a } 18 19func main() -> i64 { 20 gn_p("=== F64-GAMMA-REFL-NEG GATE: reflection kernel at NEGATIVE x vs full-domain bf_gamma ===\n" as *u8) 21 let kb: *i64 = bf_new() 22 let xb: *i64 = bf_new() 23 var maxulp: i64 = 0 24 var worst: i64 = 0 25 var npts: i64 = 0 26 var j: i64 = 0 27 while j <= 18 { 28 let kk: i64 = 2 * j + 1 // odd numerator -> never an integer 29 bf_set_int(kb, kk) 30 bf_div_small(xb, kb, 4) // magnitude (2j+1)/4 31 let mag: i64 = bf_to_f64(xb, 0, 0) 32 let xneg: i64 = mag | 0x8000000000000000 // make it negative 33 let want: i64 = bf_gamma_f64(xneg) 34 let got: i64 = nx_f64_gamma_refl(xneg) 35 let u: i64 = gn_ulp(got, want) 36 if u > maxulp { maxulp = u; worst = xneg } 37 npts = npts + 1 38 j = j + 1 39 } 40 let lfd: i64 = sys_openat_append("knowledge/status/math_engine.log" as *u8, 0x1a4) 41 if lfd >= 0 { 42 gn_f(lfd, "GAMMAREFLNEG-GATE epoch=" as *u8); gn_n(lfd, sys_now_realtime_sec()) 43 gn_f(lfd, " points=" as *u8); gn_n(lfd, npts) 44 gn_f(lfd, " max_ulp=" as *u8); gn_n(lfd, maxulp) 45 gn_f(lfd, " worst_x_bits=" as *u8); gn_n(lfd, worst) 46 gn_f(lfd, " domain=-9.25..-0.25(quarter-pts) method=euler-reflection oracle=bf_gamma_f64(full-domain)" as *u8) 47 if maxulp <= 16 { gn_f(lfd, " verdict=GREEN\n" as *u8) } else { gn_f(lfd, " verdict=RED\n" as *u8) } 48 sys_close(lfd) 49 } 50 gn_p(" points=" as *u8); gn_n(1, npts) 51 gn_p(" max_ulp=" as *u8); gn_n(1, maxulp) 52 if maxulp <= 16 { gn_p(" GAMMAREFLNEG-GATE: GREEN (gamma full-domain across reals minus integer poles)\n" as *u8); sys_exit(0); return 0 } 53 gn_p(" GAMMAREFLNEG-GATE: RED (measured gap named)\n" as *u8) 54 sys_exit(1) 55 return 1 56}