code wiki / _hdl_build / _gamma_probe.nx
_gamma_probe.nx source
↩ module page · 106 lines · 5154 B
1// _gamma_probe.nx -- DE-RISK the gamma kernel BEFORE building its emitter template
2// (design-intelligence law: validate the hard passage in isolation, don't rush a
3// multiday-bug door open). Computes Gamma(x) in plain f64 via the Stirling lnGamma
4// series with recurrence pull-up, using the EXISTING authored log + exp kernels
5// (kernel COMPOSITION proof), and compares against the known-good bigfloat oracle
6// (bf_gamma_f64) at integer + half-integer grid points. If max_ulp is small, the
7// GAMMA_STIRLING template is safe to author; if not, the probe NAMES the gap (which
8// term count / which domain) -- no emitter built on an unvalidated approach.
9// lnGamma(w) = (w-0.5)*ln(w) - w + 0.5*ln(2pi) + sum_k c_k / w^(2k-1) (w>=8)
10// Gamma(x) = exp(lnGamma(x+m)) / ((x)(x+1)...(x+m-1)) pull up to w=x+m>=8
11// Uses the spec table _pm_gamma_spec.nx (sovereign consts). Durable: GAMMAPROBE row
12// -> knowledge/status/math_engine.log. Exit 0 = approach validated (max_ulp<=8).
13// license_tier: ORIGINAL
14import "nx_syscalls.nx"
15import "nx_f64.nx"
16import "nx_f64_div.nx"
17import "nx_f64_cvt.nx"
18import "nx_bigfloat120.nx"
19import "nx_bigfloat120_div.nx"
20import "nx_bigfloat120_gamma.nx"
21import "_pm_gamma_spec.nx"
22import "_pe_f64log.nx"
23import "_pe_f64sinhcosh.nx"
24func gp_p(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 }
25func gp_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 }
26func gp_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 }
27// exp via the authored sinhcosh exp path: e^v = cosh(v)+sinh(v)
28func gp_exp(v: i64) -> i64 { return nx_f64_add(nx_f64_cosh(v), nx_f64_sinh(v)) }
29// Gamma(x) for x>0 finite, via Stirling lnGamma + pull-up to w>=8
30func gp_gamma(x: i64, c: *i64) -> i64 {
31 // pull up: multiply down-factors prod = x*(x+1)*...*(x+m-1) until w=x+m >= 8
32 let eight: i64 = c[3]
33 let one: i64 = c[0]
34 var w: i64 = x
35 var prod: i64 = one
36 var guard: i64 = 0
37 while guard < 64 {
38 if nx_f64_sub(w, eight) >> 63 == 0 { guard = 64 } else {
39 prod = nx_f64_mul(prod, w)
40 w = nx_f64_add(w, one)
41 guard = guard + 1
42 }
43 }
44 // lnGamma(w) = (w-0.5)*ln(w) - w + 0.5*ln(2pi) + sum c_k / w^(2k-1)
45 let lw: i64 = nx_f64_log(w)
46 var lg: i64 = nx_f64_mul(nx_f64_sub(w, c[1]), lw)
47 lg = nx_f64_sub(lg, w)
48 lg = nx_f64_add(lg, c[2])
49 // series: sum_{k=1..8} c[4+k] * w^-(2k-1); c[5..12] = c_k hi-first (k=8..1)
50 let winv: i64 = nx_f64_div(one, w)
51 let w2inv: i64 = nx_f64_mul(winv, winv)
52 // accumulate with smallest term first: k=8 (c[5]) ... k=1 (c[12]) is hi-first,
53 // so c[12] is k=1 (largest); build w^-(2k-1) incrementally from k=1
54 var term: i64 = winv // w^-1 (k=1)
55 var s: i64 = nx_f64_mul(c[12], term)
56 var kk: i64 = 2
57 while kk <= 8 {
58 term = nx_f64_mul(term, w2inv) // w^-(2k-1)
59 s = nx_f64_add(s, nx_f64_mul(c[12 - (kk - 1)], term))
60 kk = kk + 1
61 }
62 lg = nx_f64_add(lg, s)
63 let gw: i64 = gp_exp(lg)
64 return nx_f64_div(gw, prod)
65}
66func gp_ulp(a: i64, b: i64) -> i64 { if a >= b { return a - b } return b - a }
67func main() -> i64 {
68 gp_p("=== GAMMA PROBE: validate dd-Stirling approach vs bf oracle (de-risk template) ===\n" as *u8)
69 let c: *i64 = sys_mmap(8 * 40) as *i64
70 pm_gamma_spec_fill(c)
71 let two: *i64 = bf_new()
72 var maxulp: i64 = 0
73 var worst: i64 = 0
74 let xs: *i64 = sys_mmap(8 * 16) as *i64
75 // integer + half grid as f64 bits: 2,2.5,3,4,5,7,10
76 xs[0] = 0x4000000000000000 // 2.0
77 xs[1] = 0x4004000000000000 // 2.5
78 xs[2] = 0x4008000000000000 // 3.0
79 xs[3] = 0x4010000000000000 // 4.0
80 xs[4] = 0x4014000000000000 // 5.0
81 xs[5] = 0x401C000000000000 // 7.0
82 xs[6] = 0x4024000000000000 // 10.0
83 var npts: i64 = 7
84 var i: i64 = 0
85 while i < npts {
86 let want: i64 = bf_gamma_f64(xs[i])
87 let got: i64 = gp_gamma(xs[i], c)
88 let u: i64 = gp_ulp(got, want)
89 if u > maxulp { maxulp = u; worst = xs[i] }
90 i = i + 1
91 }
92 let lfd: i64 = sys_openat_append("knowledge/status/math_engine.log" as *u8, 0x1a4)
93 if lfd >= 0 {
94 gp_f(lfd, "GAMMAPROBE epoch=" as *u8); gp_n(lfd, sys_now_realtime_sec())
95 gp_f(lfd, " points=" as *u8); gp_n(lfd, npts)
96 gp_f(lfd, " max_ulp=" as *u8); gp_n(lfd, maxulp)
97 gp_f(lfd, " worst_x=" as *u8); gp_n(lfd, worst)
98 if maxulp <= 8 { gp_f(lfd, " verdict=APPROACH-VALIDATED\n" as *u8) } else { gp_f(lfd, " verdict=GAP-NAMED\n" as *u8) }
99 sys_close(lfd)
100 }
101 gp_p(" max_ulp=" as *u8); gp_n(1, maxulp)
102 if maxulp <= 8 { gp_p(" -- GAMMA PROBE: APPROACH VALIDATED (template safe to author)\n" as *u8); sys_exit(0); return 0 }
103 gp_p(" -- GAMMA PROBE: GAP (term count or domain -- refine before emitter)\n" as *u8)
104 sys_exit(1)
105 return 1
106}