code wiki / _hdl_build / _pe_f64igammaq.nx
_pe_f64igammaq.nx source
↩ module page · 39 lines · 1426 B
1// _pe_f64igammaq.nx -- integer-order regularized upper incomplete gamma (DLMF 8.4.7,
2// the Poisson/Erlang/chi-square CDF):
3// Q(n,x) = Gamma(n,x)/Gamma(n) = e^(-x) * SUM_{k=0}^{n-1} x^k / k! (n >= 1)
4// P(n,x) = 1 - Q(n,x)
5// CLOSED FORM (finite sum) for integer order -- no series-vs-CF regime split, all terms
6// positive, self-contained (only nx_f64_exp + f64 arithmetic). term_k = term_{k-1}*x/k.
7// f64 software-IEEE throughout.
8// license_tier: ORIGINAL
9//
10// module: nishi-core.math.f64_igammaq
11// depends: nishi-core.math.f64, nishi-core.math.f64_exp
12// capability: F64_INCOMPLETE_GAMMA_INTEGER
13import "nx_syscalls.nx"
14import "nx_f64.nx"
15import "nx_f64_div.nx"
16import "nx_f64_exp.nx"
17
18const IGQ_ONE: i64 = 0x3FF0000000000000 // 1.0
19
20// Q(n,x), n>=1 integer, x>=0 (f64 bits)
21func nx_f64_igammaq(n: i64, x: i64) -> i64 {
22 var term: i64 = IGQ_ONE // x^0/0! = 1
23 var sum: i64 = IGQ_ONE
24 var kf: i64 = 0 // +0.0, becomes k as f64
25 var k: i64 = 1
26 while k < n {
27 kf = nx_f64_add(kf, IGQ_ONE) // kf = k
28 term = nx_f64_div(nx_f64_mul(term, x), kf) // term *= x/k
29 sum = nx_f64_add(sum, term)
30 k = k + 1
31 }
32 let e: i64 = nx_f64_exp(nx_f64_neg(x)) // e^(-x)
33 return nx_f64_mul(e, sum)
34}
35
36// P(n,x) = 1 - Q(n,x)
37func nx_f64_igammap(n: i64, x: i64) -> i64 {
38 return nx_f64_sub(IGQ_ONE, nx_f64_igammaq(n, x))
39}