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}