nx_bigfloat120_gamma.nx source
↩ module page · 289 lines · 10481 B
1// nx_bigfloat120_gamma.nx -- SOVEREIGN gamma oracle on the 120-bit bigfloat
2// (ME2 rung 1: DLMF ch.5 enters the engine). Composition of blessed organs:
3// ln Gamma(z) = (z - 1/2) ln z - z + ln(2 pi)/2
4// + sum_{n=1..17} B_2n / (2n (2n-1) z^(2n-1)) [z >= 30]
5// with the INTEGER Bernoulli rationals B_2..B_34 (DLMF 24.2 -- mathematical
6// constants like k!; B_36's numerator exceeds i64 and is not needed: at z = 30
7// the n=17 tail is < 2^-138, 28 bits past the oracle's 2^-110 contract).
8// x < 30 shifts up by the recurrence Gamma(x) = Gamma(x+k)/(x(x+1)..(x+k-1));
9// negative non-integers reflect through pi/(sin(pi x) Gamma(1-x)) over the
10// blessed trig organ. exp from the pow organ; ln rebuilt bigfloat-in (the
11// atanh form; exponent split is exact).
12//
13// C99 tgamma specials: +-0 -> +-inf, negative integers and -inf -> NaN,
14// +inf -> +inf, x > 172 -> +inf (true overflow 171.624), x < -184 -> +-0
15// (|Gamma| below the subnormal floor; sign by interval parity).
16//
17// Self-anchor (see _bf_gamma_gate_authored): EXACT factorials to 18!,
18// the recurrence identity across different shift depths, Gamma(1/2)^2 = pi,
19// and the reflection product r(x) r(1-x) sin(pi x) = pi on (0,1).
20// license_tier: ORIGINAL
21
22import "nx_syscalls.nx"
23import "nx_tier.nx"
24import "nx_bigfloat120.nx"
25import "nx_bigfloat120_div.nx"
26import "nx_bigfloat120_trig.nx"
27import "nx_bigfloat120_exp.nx"
28import "nx_bigfloat120_ln.nx"
29import "nx_bigfloat120_pow.nx"
30const K_MAGIC_2730: i64 = 2730
31const K_MAGIC_3617: i64 = 3617
32const K_MAGIC_43867: i64 = 43867
33const K_MAGIC_174611: i64 = 174611
34const K_MAGIC_854513: i64 = 854513
35const K_MAGIC_236364091: i64 = 236364091
36const K_MAGIC_8553103: i64 = 8553103
37const K_MAGIC_23749461029: i64 = 23749461029
38const K_MAGIC_8615841276005: i64 = 8615841276005
39const K_MAGIC_14322: i64 = 14322
40const K_MAGIC_7709321041217: i64 = 7709321041217
41const K_MAGIC_2577687858367: i64 = 2577687858367
42const K_MAGIC_2047: i64 = 2047
43
44// ln of a positive bigfloat x >= 1: exponent split (exact) + atanh series on
45// the [1,2) mantissa (s = (m-1)/(m+1) <= 1/3 -> 38 terms < 2^-120 tail).
46func bfg_ln_big(out: *i64, xb: *i64) -> i64 {
47 let e: i64 = xb[0]
48 let m: *i64 = bf_new()
49 bf_copy(m, xb)
50 m[0] = 0 // m in [1, 2)
51 let one: *i64 = bf_new()
52 bf_set_int(one, 1)
53 let lnm: *i64 = bf_new() // |ln m| (m >= 1 -> ln m >= 0)
54 if bf_cmp(m, one) != 0 {
55 let num: *i64 = bf_new()
56 bf_sub(num, m, one)
57 let den: *i64 = bf_new()
58 bf_add(den, m, one)
59 let s: *i64 = bf_new()
60 bf_div(s, num, den)
61 let z: *i64 = bf_new()
62 bf_mul(z, s, s)
63 let term: *i64 = bf_new()
64 bf_copy(term, s)
65 let sum: *i64 = bf_new()
66 bf_copy(sum, s)
67 let tmul: *i64 = bf_new()
68 let tdiv: *i64 = bf_new()
69 let snew: *i64 = bf_new()
70 var j: i64 = 1
71 while j <= 38 {
72 bf_mul(tmul, term, z)
73 bf_copy(term, tmul)
74 bf_div_small(tdiv, term, 2 * j + 1)
75 bf_add(snew, sum, tdiv)
76 bf_copy(sum, snew)
77 j = j + 1
78 }
79 bf_copy(lnm, sum)
80 lnm[0] = lnm[0] + 1 // * 2
81 }
82 if e == 0 { return bf_copy(out, lnm) }
83 let ln2: *i64 = bf_new()
84 bf_ln2(ln2)
85 let kl: *i64 = bf_new()
86 bf_mul_small(kl, ln2, e)
87 return bf_add(out, kl, lnm) // e >= 0 for x >= 1
88}
89
90// sin of a positive bigfloat (|xb| < 2^20): mirrors bf_sin_f64's quadrant
91// logic with a bigfloat argument; returns the sign bit, |sin| in out.
92func bfg_sin_big(out: *i64, xb: *i64) -> i64 {
93 if bf_is_zero(xb) == 1 { return bf_copy(out, xb) }
94 let r: *i64 = bf_new()
95 let kq: *i64 = sys_mmap(16) as *i64
96 _bf_trig_reduce(r, kq, xb)
97 let sv: *i64 = bf_new()
98 let cv: *i64 = bf_new()
99 bf_sin_r(sv, r)
100 bf_cos_r(cv, r)
101 let kk: i64 = kq[0]
102 let rneg: i64 = kq[1]
103 var outsig: i64 = 0
104 var usecos: i64 = 0
105 if kk == 0 { outsig = rneg }
106 if kk == 1 { usecos = 1; outsig = 0 }
107 if kk == 2 { outsig = 1 - rneg }
108 if kk == 3 { usecos = 1; outsig = 1 }
109 if usecos == 1 { bf_copy(out, cv) } else { bf_copy(out, sv) }
110 return outsig
111}
112
113// ln Gamma(z) for bigfloat z >= 30 (Stirling + 17 Bernoulli terms)
114func bfg_lngamma_big(out: *i64, z: *i64) -> i64 {
115 let one: *i64 = bf_new()
116 bf_set_int(one, 1)
117 let half: *i64 = bf_new()
118 bf_copy(half, one)
119 half[0] = half[0] - 1
120 // A = (z - 1/2) * ln z (z >= 30: every piece positive)
121 let zm: *i64 = bf_new()
122 bf_sub(zm, z, half)
123 let lz: *i64 = bf_new()
124 bfg_ln_big(lz, z)
125 let a: *i64 = bf_new()
126 bf_mul(a, zm, lz)
127 // C = ln(2 pi)/2 = (ln 2 + ln pi)/2
128 let pi: *i64 = bf_new()
129 bf_pi(pi)
130 let lnpi: *i64 = bf_new()
131 bfg_ln_big(lnpi, pi)
132 let ln2: *i64 = bf_new()
133 bf_ln2(ln2)
134 let c: *i64 = bf_new()
135 bf_add(c, ln2, lnpi)
136 c[0] = c[0] - 1
137 // S = sum B_2n/(2n(2n-1) z^(2n-1)), signs (+,-,+,..): pos/neg accumulators.
138 // |B_2n| numerator/denominator pairs, DLMF 24.2 (all numerators < 2^53).
139 let bn: *i64 = sys_mmap(8 * 20) as *i64
140 let bd: *i64 = sys_mmap(8 * 20) as *i64
141 bn[1] = 1; bd[1] = 6
142 bn[2] = 1; bd[2] = 30
143 bn[3] = 1; bd[3] = 42
144 bn[4] = 1; bd[4] = 30
145 bn[5] = 5; bd[5] = 66
146 bn[6] = 691; bd[6] = K_MAGIC_2730
147 bn[7] = 7; bd[7] = 6
148 bn[8] = K_MAGIC_3617; bd[8] = 510
149 bn[9] = K_MAGIC_43867; bd[9] = 798
150 bn[10] = K_MAGIC_174611; bd[10] = 330
151 bn[11] = K_MAGIC_854513; bd[11] = 138
152 bn[12] = K_MAGIC_236364091; bd[12] = K_MAGIC_2730
153 bn[13] = K_MAGIC_8553103; bd[13] = 6
154 bn[14] = K_MAGIC_23749461029; bd[14] = 870
155 bn[15] = K_MAGIC_8615841276005; bd[15] = K_MAGIC_14322
156 bn[16] = K_MAGIC_7709321041217; bd[16] = 510
157 bn[17] = K_MAGIC_2577687858367; bd[17] = 6
158 let z2: *i64 = bf_new()
159 bf_mul(z2, z, z)
160 let zp: *i64 = bf_new()
161 bf_copy(zp, z) // z^(2n-1), n=1 -> z
162 let pos: *i64 = bf_new()
163 let neg: *i64 = bf_new()
164 let t1: *i64 = bf_new()
165 let t2: *i64 = bf_new()
166 let t3: *i64 = bf_new()
167 let snew: *i64 = bf_new()
168 var n: i64 = 1
169 while n <= 17 {
170 bf_set_int(t1, bn[n])
171 bf_div_small(t2, t1, bd[n])
172 bf_div_small(t1, t2, (2 * n) * (2 * n - 1))
173 bf_div(t3, t1, zp)
174 if (n & 1) == 1 {
175 bf_add(snew, pos, t3)
176 bf_copy(pos, snew)
177 } else {
178 bf_add(snew, neg, t3)
179 bf_copy(neg, snew)
180 }
181 if n < 17 {
182 bf_mul(t1, zp, z2)
183 bf_copy(zp, t1)
184 }
185 n = n + 1
186 }
187 let s: *i64 = bf_new()
188 bf_sub(s, pos, neg) // first term dominates
189 // out = A - z + C + S (A > z + C + S margin at z >= 30: A >= 100)
190 let acs: *i64 = bf_new()
191 bf_add(acs, a, c)
192 let acs2: *i64 = bf_new()
193 bf_add(acs2, acs, s)
194 return bf_sub(out, acs2, z)
195}
196
197// Gamma(x) for positive bigfloat x (any magnitude with lnGamma < ~1090)
198func bfg_gamma_pos(out: *i64, xb: *i64) -> i64 {
199 let thirty: *i64 = bf_new()
200 bf_set_int(thirty, 30)
201 let z: *i64 = bf_new()
202 bf_copy(z, xb)
203 let prod: *i64 = bf_new()
204 bf_set_int(prod, 1)
205 let one: *i64 = bf_new()
206 bf_set_int(one, 1)
207 let t: *i64 = bf_new()
208 var guard: i64 = 0
209 var go: i64 = 1
210 while go == 1 {
211 if guard >= 32 { go = 0 } // mathematical bound: 30 shifts
212 if bf_cmp(z, thirty) >= 0 { go = 0 }
213 if go == 1 {
214 bf_mul(t, prod, z)
215 bf_copy(prod, t)
216 bf_add(t, z, one)
217 bf_copy(z, t)
218 guard = guard + 1
219 }
220 }
221 let lg: *i64 = bf_new()
222 bfg_lngamma_big(lg, z)
223 let g: *i64 = bf_new()
224 bfp_exp_pm(g, lg, 0) // lnGamma > 0 at z >= 30
225 return bf_div(out, g, prod)
226}
227
228// tgamma(raw f64) -> raw f64, C99 semantics.
229func bf_gamma_f64(x: i64) -> i64 {
230 let ax: i64 = x & 0x7FFFFFFFFFFFFFFF
231 let ef: i64 = (x >> 52) & 0x7FF
232 let sgn: i64 = (x >> 63) & 1
233 if ef == K_MAGIC_2047 {
234 if (ax & 0x000FFFFFFFFFFFFF) != 0 { return 0x7FF8000000000000 }
235 if sgn == 1 { return 0x7FF8000000000000 } // tgamma(-inf) = NaN
236 return 0x7FF0000000000000 // tgamma(+inf) = +inf
237 }
238 if ax == 0 {
239 if sgn == 1 { return 0xFFF0000000000000 } // tgamma(-0) = -inf
240 return 0x7FF0000000000000 // tgamma(+0) = +inf
241 }
242 if sgn == 0 {
243 // overflow cut: Gamma(172) > maxdouble (true boundary 171.624);
244 // positive raws order as integers, so one raw compare decides
245 if ax >= 0x4065800000000000 { return 0x7FF0000000000000 } // >= 172
246 let xb: *i64 = bf_new()
247 bf_set_f64(xb, ax)
248 let g: *i64 = bf_new()
249 bfg_gamma_pos(g, xb)
250 return bf_to_f64(g, 0, 0)
251 }
252 // negative x: poles at integers; reflection elsewhere
253 let ycls: i64 = pw_int_class(ax)
254 if ycls != 0 { return 0x7FF8000000000000 } // negative integer
255 // |Gamma(x)| < min subnormal for x < -184: signed zero by interval parity
256 let e2: i64 = ef - 1023
257 var hardz: i64 = 0
258 if ax >= 0x4067000000000000 { hardz = 1 } // |x| >= 184
259 // interval parity: x in (-(n+1), -n): sign negative iff n even
260 // n = floor(|x|): for f64 with e2 >= 0, extract integer part
261 var ipart: i64 = 0
262 if e2 >= 0 {
263 if e2 <= 52 { ipart = ((ax & 0x000FFFFFFFFFFFFF) | 0x0010000000000000) >> (52 - e2) }
264 }
265 var rsign: i64 = 0
266 if (ipart & 1) == 0 { rsign = 1 } // n even -> negative
267 if hardz == 1 { return rsign << 63 }
268 // Gamma(x) = pi / (sin(pi x) * Gamma(1 - x)); x = -|x|:
269 // sin(pi x) = -sin(pi |x|): overall sign = 1 ^ sign(sin(pi |x|))
270 let axb: *i64 = bf_new()
271 bf_set_f64(axb, ax)
272 let pi: *i64 = bf_new()
273 bf_pi(pi)
274 let px: *i64 = bf_new()
275 bf_mul(px, pi, axb)
276 let sv: *i64 = bf_new()
277 let ssig: i64 = bfg_sin_big(sv, px)
278 let one: *i64 = bf_new()
279 bf_set_int(one, 1)
280 let omx: *i64 = bf_new()
281 bf_add(omx, one, axb) // 1 - x = 1 + |x|
282 let g1: *i64 = bf_new()
283 bfg_gamma_pos(g1, omx)
284 let den: *i64 = bf_new()
285 bf_mul(den, sv, g1)
286 let res: *i64 = bf_new()
287 bf_div(res, pi, den)
288 return bf_to_f64(res, 1 - ssig, 0)
289}