code wiki / _hdl_build / _pe_f64gamma.nx
_pe_f64gamma.nx source
↩ module page · 74 lines · 3053 B
1// AUTHORED BY THE NISHI BUILDER (pattern: MATH_KERNEL / GAMMA_STIRLING) -- no Claude core logic.
2// dd lnGamma (Stirling + pull-up), composed from authored log+sinhcosh kernels; spec=sovereign bigfloat.
3// Validated max_ulp=2 vs bf oracle. Rung-1 x>=0.5; x<0.5 reflection = LOUD NAN (named debt).
4import "nx_syscalls.nx"
5import "nx_tier.nx"
6import "nx_f64.nx"
7import "nx_f64_div.nx"
8import "nx_f64_cvt.nx"
9import "_pe_f64log.nx"
10import "_pe_f64sinhcosh.nx"
11func _gm_expf(v: i64) -> i64 { return nx_f64_add(nx_f64_cosh(v), nx_f64_sinh(v)) }
12func _gm_lndd(w: i64, out: *i64) -> i64 {
13 let hi: i64 = nx_f64_log(w)
14 let eh: i64 = _gm_expf(hi)
15 out[0] = hi
16 out[1] = nx_f64_div(nx_f64_sub(w, eh), w)
17 return 0
18}
19func nx_f64_gamma(x: i64) -> i64 {
20 let cls: nx_int = nx_f64_classify(x)
21 if cls == NX_F64_CLS_NAN { return NX_F64_NAN_RAW }
22 if nx_f64_sign(x) == 1 { return NX_F64_NAN_RAW }
23 if cls == NX_F64_CLS_INF { return x }
24 if (x & 0x7FFFFFFFFFFFFFFF) < 4602678819172646912 { return NX_F64_NAN_RAW }
25 if nx_f64_sub(x, 4640242534423986176) >> 63 == 0 { return 0x7FF0000000000000 }
26 var w: i64 = x
27 var prod: i64 = 4607182418800017408
28 var guard: i64 = 0
29 while guard < 64 {
30 if nx_f64_sub(w, 4620693217682128896) >> 63 == 0 { guard = 64 } else {
31 prod = nx_f64_mul(prod, w)
32 w = nx_f64_add(w, 4607182418800017408)
33 guard = guard + 1
34 }
35 }
36 let ln: *i64 = sys_mmap(32) as *i64
37 _gm_lndd(w, ln)
38 let wm: i64 = nx_f64_sub(w, 4602678819172646912)
39 let mp: *i64 = sys_mmap(32) as *i64
40 _sh_mul2(wm, ln[0], mp)
41 var a_hi: i64 = mp[0]
42 var a_lo: i64 = nx_f64_add(mp[1], nx_f64_mul(wm, ln[1]))
43 let acc: *i64 = sys_mmap(32) as *i64
44 _sh_2sum(a_hi, nx_f64_neg(w), acc)
45 var g_hi: i64 = acc[0]
46 var g_lo: i64 = nx_f64_add(acc[1], a_lo)
47 _sh_2sum(g_hi, 4606452282016710324, acc)
48 g_hi = acc[0]
49 g_lo = nx_f64_add(g_lo, acc[1])
50 let winv: i64 = nx_f64_div(4607182418800017408, w)
51 let w2inv: i64 = nx_f64_mul(winv, winv)
52 var term: i64 = winv
53 var ser: i64 = nx_f64_mul(4590669220166325589, term)
54 term = nx_f64_mul(term, w2inv)
55 ser = nx_f64_add(ser, nx_f64_mul(0 - 4654820494858425321, term))
56 term = nx_f64_mul(term, w2inv)
57 ser = nx_f64_add(ser, nx_f64_mul(4560459359808757786, term))
58 term = nx_f64_mul(term, w2inv)
59 ser = nx_f64_add(ser, nx_f64_mul(0 - 4664742711180314604, term))
60 term = nx_f64_mul(term, w2inv)
61 ser = nx_f64_add(ser, nx_f64_mul(4560903004447375139, term))
62 term = nx_f64_mul(term, w2inv)
63 ser = nx_f64_add(ser, nx_f64_mul(0 - 4656886181880316803, term))
64 term = nx_f64_mul(term, w2inv)
65 ser = nx_f64_add(ser, nx_f64_mul(4574040544619111450, term))
66 term = nx_f64_mul(term, w2inv)
67 ser = nx_f64_add(ser, nx_f64_mul(0 - 4639197419445202024, term))
68 _sh_2sum(g_hi, ser, acc)
69 g_hi = acc[0]
70 g_lo = nx_f64_add(g_lo, acc[1])
71 let eg: i64 = _gm_expf(g_hi)
72 let gw: i64 = nx_f64_add(eg, nx_f64_mul(eg, g_lo))
73 return nx_f64_div(gw, prod)
74}