_gamma_dd_probe.nx source

1// _gamma_dd_probe.nx -- the DD-precision gamma probe (ME2-GAMMA-DD). The plain-f64 2// probe found 13 ulp; the named requirement was double-double lnGamma. This composes 3// the authored kernels' dd primitives (_sh_expdd/_sh_mul2/_sh_2sum/_sh_scale from the 4// sinhcosh module + nx_f64_log) into a dd lnGamma, and re-tests vs bf_gamma_f64. If 5// max_ulp <= 4 the GAMMA_STIRLING template is safe to author at dd; else the gap is 6// refined further. No emitter built on an unvalidated approach (design-intelligence). 7// ln_dd(w) : hi = nx_f64_log(w); lo = (w - exp_dd(hi))/w (Newton correction) 8// lnGamma_dd : (w-0.5)*ln_dd(w) [dd product] - w + 0.5ln2pi + series, kept dd 9// Gamma : exp(lnGamma_hi)*(1+lnGamma_lo) / pull-up-product 10// license_tier: ORIGINAL 11import "nx_syscalls.nx" 12import "nx_f64.nx" 13import "nx_f64_div.nx" 14import "nx_f64_cvt.nx" 15import "nx_bigfloat120.nx" 16import "nx_bigfloat120_div.nx" 17import "nx_bigfloat120_gamma.nx" 18import "_pm_gamma_spec.nx" 19import "_pe_f64log.nx" 20import "_pe_f64sinhcosh.nx" 21func gd_p(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 } 22func gd_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 } 23func gd_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 } 24// exp(v) via the PROVEN authored path (cosh+sinh) -- the plain-f64 probe confirmed this 25// gives correct values; the _sh_expdd reconstruction was the structural bug the first 26// dd run caught (484e9 ulp). Root-cause fix, not a patch: use the verified exp. 27func gd_expf(v: i64) -> i64 { return nx_f64_add(nx_f64_cosh(v), nx_f64_sinh(v)) } 28// ln_dd(w) -> out[0]=hi out[1]=lo (Newton: lo = (w - exp(hi))/w) 29func gd_lndd(w: i64, out: *i64) -> i64 { 30 let hi: i64 = nx_f64_log(w) 31 let eh: i64 = gd_expf(hi) 32 let resid: i64 = nx_f64_sub(w, eh) 33 out[0] = hi 34 out[1] = nx_f64_div(resid, w) 35 return 0 36} 37// Gamma(x) for x>0 via dd Stirling lnGamma + pull-up to w>=8 38func gd_gamma(x: i64, c: *i64) -> i64 { 39 let eight: i64 = c[3] 40 let one: i64 = c[1] 41 let half: i64 = c[1] 42 var w: i64 = x 43 var prod: i64 = c[0] 44 var guard: i64 = 0 45 while guard < 64 { 46 if nx_f64_sub(w, eight) >> 63 == 0 { guard = 64 } else { 47 prod = nx_f64_mul(prod, w) 48 w = nx_f64_add(w, c[0]) 49 guard = guard + 1 50 } 51 } 52 // ln_dd(w) 53 let ln: *i64 = sys_mmap(32) as *i64 54 gd_lndd(w, ln) 55 // A = (w-0.5) * ln_dd(w), dd: _sh_mul2(w-0.5, ln_hi) + (w-0.5)*ln_lo into lo 56 let wm: i64 = nx_f64_sub(w, c[1]) 57 let mp: *i64 = sys_mmap(32) as *i64 58 _sh_mul2(wm, ln[0], mp) // mp[0]=hi mp[1]=lo of (w-0.5)*ln_hi 59 var a_hi: i64 = mp[0] 60 var a_lo: i64 = nx_f64_add(mp[1], nx_f64_mul(wm, ln[1])) 61 // lnGamma = A - w + 0.5ln2pi + series, accumulate dd via _sh_2sum on the hi chain 62 let acc: *i64 = sys_mmap(32) as *i64 63 _sh_2sum(a_hi, nx_f64_neg(w), acc) // A + (-w) 64 var g_hi: i64 = acc[0] 65 var g_lo: i64 = nx_f64_add(acc[1], a_lo) 66 _sh_2sum(g_hi, c[2], acc) // + 0.5*ln(2pi) 67 g_hi = acc[0]; g_lo = nx_f64_add(g_lo, acc[1]) 68 // series sum c_k / w^(2k-1), k=1..8 (c[12]=k1 .. c[5]=k8 hi-first), small -> into lo 69 let winv: i64 = nx_f64_div(c[0], w) 70 let w2inv: i64 = nx_f64_mul(winv, winv) 71 var term: i64 = winv 72 var ser: i64 = nx_f64_mul(c[12], term) 73 var kk: i64 = 2 74 while kk <= 8 { 75 term = nx_f64_mul(term, w2inv) 76 ser = nx_f64_add(ser, nx_f64_mul(c[12 - (kk - 1)], term)) 77 kk = kk + 1 78 } 79 // the Stirling series (~0.01 for w=8) is a SIGNIFICANT lnGamma term, NOT a tiny dd 80 // correction -- it belongs in g_hi (the exponent), not g_lo. Folding it into g_lo 81 // and linearizing exp*(1+g_lo) was the ~0.5-abs bug (g_lo must stay ~2^-52). 82 _sh_2sum(g_hi, ser, acc) 83 g_hi = acc[0]; g_lo = nx_f64_add(g_lo, acc[1]) 84 // Gamma = exp(g_hi) * (1 + g_lo) / prod (g_lo now genuinely tiny -> linearization ok) 85 let eg: i64 = gd_expf(g_hi) 86 let gw: i64 = nx_f64_add(eg, nx_f64_mul(eg, g_lo)) 87 return nx_f64_div(gw, prod) 88} 89func gd_ulp(a: i64, b: i64) -> i64 { if a >= b { return a - b } return b - a } 90func main() -> i64 { 91 gd_p("=== GAMMA DD PROBE: double-double lnGamma vs bf oracle (the named requirement) ===\n" as *u8) 92 let c: *i64 = sys_mmap(8 * 40) as *i64 93 pm_gamma_spec_fill(c) 94 let xs: *i64 = sys_mmap(8 * 16) as *i64 95 xs[0] = 0x4000000000000000 96 xs[1] = 0x4004000000000000 97 xs[2] = 0x4008000000000000 98 xs[3] = 0x4010000000000000 99 xs[4] = 0x4014000000000000 100 xs[5] = 0x401C000000000000 101 xs[6] = 0x4024000000000000 102 var npts: i64 = 7 103 var maxulp: i64 = 0 104 var worst: i64 = 0 105 var i: i64 = 0 106 while i < npts { 107 let want: i64 = bf_gamma_f64(xs[i]) 108 let got: i64 = gd_gamma(xs[i], c) 109 let u: i64 = gd_ulp(got, want) 110 if u > maxulp { maxulp = u; worst = xs[i] } 111 i = i + 1 112 } 113 let lfd: i64 = sys_openat_append("knowledge/status/math_engine.log" as *u8, 0x1a4) 114 if lfd >= 0 { 115 gd_f(lfd, "GAMMADDPROBE epoch=" as *u8); gd_n(lfd, sys_now_realtime_sec()) 116 gd_f(lfd, " points=" as *u8); gd_n(lfd, npts) 117 gd_f(lfd, " max_ulp=" as *u8); gd_n(lfd, maxulp) 118 gd_f(lfd, " worst_x=" as *u8); gd_n(lfd, worst) 119 if maxulp <= 4 { gd_f(lfd, " verdict=DD-VALIDATED\n" as *u8) } else { gd_f(lfd, " verdict=DD-GAP\n" as *u8) } 120 sys_close(lfd) 121 } 122 gd_p(" max_ulp=" as *u8); gd_n(1, maxulp) 123 if maxulp <= 4 { gd_p(" -- GAMMA DD PROBE: VALIDATED (author GAMMA_STIRLING template at dd)\n" as *u8); sys_exit(0); return 0 } 124 gd_p(" -- GAMMA DD PROBE: still gap (refine dd-ln correction or series)\n" as *u8) 125 sys_exit(1) 126 return 1 127}