pxfm/gamma/
tgamma.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 8/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::common::{f_fmla, is_integer};
30use crate::double_double::DoubleDouble;
31use crate::floor::FloorFinite;
32use crate::logs::fast_log_dd;
33use crate::polyeval::{f_polyeval6, f_polyeval7, f_polyeval8};
34use crate::pow_exec::exp_dd_fast;
35use crate::sincospi::f_fast_sinpi_dd;
36
37/// Computes gamma(x)
38///
39/// ulp 1
40pub fn f_tgamma(x: f64) -> f64 {
41    let x_a = f64::from_bits(x.to_bits() & 0x7fff_ffff_ffff_ffff);
42
43    if !x.is_normal() {
44        if x == 0.0 {
45            return f64::INFINITY;
46        }
47
48        if x.is_nan() {
49            return x + x;
50        }
51
52        if x.is_infinite() {
53            if x.is_sign_negative() {
54                return f64::NAN;
55            }
56            return x;
57        }
58    }
59
60    if x >= 171.624 {
61        return f64::INFINITY;
62    }
63
64    if is_integer(x) {
65        if x < 0. {
66            return f64::NAN;
67        }
68        if x < 38. {
69            let mut t = DoubleDouble::new(0., 1.);
70            let k = x as i64;
71            let mut x0 = 1.0;
72            for _i in 1..k {
73                t = DoubleDouble::quick_mult_f64(t, x0);
74                t = DoubleDouble::from_exact_add(t.hi, t.lo);
75                x0 += 1.0;
76            }
77            return t.to_f64();
78        }
79    }
80
81    if x <= -184.0 {
82        /* negative non-integer */
83        /* For x <= -184, x non-integer, |gamma(x)| < 2^-1078.  */
84        static SIGN: [f64; 2] = [
85            f64::from_bits(0x0010000000000000),
86            f64::from_bits(0x8010000000000000),
87        ];
88        let k = x as i64;
89        return f64::from_bits(0x0010000000000000) * SIGN[((k & 1) != 0) as usize];
90    }
91
92    const EULER_DD: DoubleDouble =
93        DoubleDouble::from_bit_pair((0xbc56cb90701fbfab, 0x3fe2788cfc6fb619));
94
95    if x_a < 0.006 {
96        if x_a.to_bits() < (0x71e0000000000000u64 >> 1) {
97            // |x| < 0x1p-112
98            return 1. / x;
99        }
100        if x_a < 2e-10 {
101            // x is tiny then Gamma(x) = 1/x - euler
102            let p = DoubleDouble::full_dd_sub(DoubleDouble::from_quick_recip(x), EULER_DD);
103            return p.to_f64();
104        } else if x_a < 2e-6 {
105            // x is small then Gamma(x) = 1/x - euler + a2*x
106            // a2 = 1/12 * (6 * euler^2 + pi^2)
107            const A2: DoubleDouble =
108                DoubleDouble::from_bit_pair((0x3c8dd92b465a8221, 0x3fefa658c23b1578));
109            let rcp = DoubleDouble::from_quick_recip(x);
110            let p = DoubleDouble::full_dd_add(DoubleDouble::mul_f64_add(A2, x, -EULER_DD), rcp);
111            return p.to_f64();
112        }
113
114        // Laurent series of Gamma(x)
115        const C: [(u64, u64); 8] = [
116            (0x3c8dd92b465a8221, 0x3fefa658c23b1578),
117            (0x3c53a4f483760950, 0xbfed0a118f324b63),
118            (0x3c7fabe4f7369157, 0x3fef6a51055096b5),
119            (0x3c8c9fc795fc6142, 0xbfef6c80ec38b67b),
120            (0xbc5042339d62e721, 0x3fefc7e0a6eb310b),
121            (0xbc86fd0d8bdc0c1e, 0xbfefdf3f157b7a39),
122            (0xbc89b912df09395d, 0x3feff07b5a17ff6c),
123            (0x3c4e626faf780ff9, 0xbfeff803d68a0bd4),
124        ];
125        let rcp = DoubleDouble::from_quick_recip(x);
126        let mut p = DoubleDouble::mul_f64_add(
127            DoubleDouble::from_bit_pair(C[7]),
128            x,
129            DoubleDouble::from_bit_pair(C[6]),
130        );
131        p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[5]));
132        p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[4]));
133        p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[3]));
134        p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[2]));
135        p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[1]));
136        p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[0]));
137        let z = DoubleDouble::mul_f64_add(p, x, DoubleDouble::full_dd_sub(rcp, EULER_DD));
138        return z.to_f64();
139    }
140
141    let mut fact = DoubleDouble::new(0., 0.0f64);
142    let mut parity = 1.0;
143    const PI: DoubleDouble = DoubleDouble::from_bit_pair((0x3ca1a62633145c07, 0x400921fb54442d18));
144    let mut dy = DoubleDouble::new(0., x);
145    let mut result: DoubleDouble;
146
147    // reflection
148    if dy.hi < 0. {
149        if dy.hi.floor_finite() == dy.hi {
150            return f64::NAN;
151        }
152        dy.hi = f64::from_bits(dy.hi.to_bits() & 0x7fff_ffff_ffff_ffff);
153        let y1 = x_a.floor_finite();
154        let fraction = x_a - y1;
155        if fraction != 0.0
156        // is it an integer?
157        {
158            // is it odd or even?
159            if y1 != (y1 * 0.5).floor_finite() * 2.0 {
160                parity = -1.0;
161            }
162            fact = DoubleDouble::div(-PI, f_fast_sinpi_dd(fraction));
163            fact = DoubleDouble::from_exact_add(fact.hi, fact.lo);
164            dy = DoubleDouble::from_full_exact_add(dy.hi, 1.0);
165        }
166    }
167
168    if dy.hi < 12.0 {
169        let y1 = dy;
170        let z: DoubleDouble;
171        let mut n = 0;
172        // x is in (0.06, 1.0).
173        if dy.hi < 1.0 {
174            z = dy;
175            dy = DoubleDouble::full_add_f64(dy, 1.0);
176        } else
177        // x is in [1.0, max].
178        {
179            n = dy.hi as i32 - 1;
180            dy = DoubleDouble::full_add_f64(dy, -n as f64);
181            z = DoubleDouble::full_add_f64(dy, -1.0);
182        }
183
184        // Gamma(x+1) on [1;2] generated by Wolfram Mathematica:
185        // <<FunctionApproximations`
186        // ClearAll["Global`*"]
187        // f[x_]:=Gamma[x+1]
188        // {err0, approx}=MiniMaxApproximation[f[z],{z,{0,1 },9,8},WorkingPrecision->90]
189        // num=Numerator[approx][[1]];
190        // den=Denominator[approx][[1]];
191        // poly=den;
192        // coeffs=CoefficientList[poly,z];
193        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
194        let ps_num = f_polyeval8(
195            z.hi,
196            f64::from_bits(0x3fb38b80568a42aa),
197            f64::from_bits(0xbf8e7685b00d63a6),
198            f64::from_bits(0xbf80629ed2c48f1a),
199            f64::from_bits(0xbf6dfc4cdbcee96a),
200            f64::from_bits(0xbf471816939dc42b),
201            f64::from_bits(0xbf24bede7d8b3c20),
202            f64::from_bits(0xbeef56936d891e42),
203            f64::from_bits(0xbec3e2b405350813),
204        );
205        let mut p_num = DoubleDouble::mul_f64_add(
206            z,
207            ps_num,
208            DoubleDouble::from_bit_pair((0xbc700f441aea2edb, 0x3fd218bdde7878b8)),
209        );
210        p_num = DoubleDouble::mul_add(
211            z,
212            p_num,
213            DoubleDouble::from_bit_pair((0x3b7056a45a2fa50e, 0x3ff0000000000000)),
214        );
215        p_num = DoubleDouble::from_exact_add(p_num.hi, p_num.lo);
216
217        let ps_den = f_polyeval7(
218            z.hi,
219            f64::from_bits(0xbfdaa4f09f0caab1),
220            f64::from_bits(0xbfc960ba48423f9d),
221            f64::from_bits(0x3fb6873b64e8ccd6),
222            f64::from_bits(0x3f69ea1ca5b8a225),
223            f64::from_bits(0xbf77b166f68a2e63),
224            f64::from_bits(0x3f4fd1eff9193728),
225            f64::from_bits(0xbf0c1a43f4985c97),
226        );
227        let mut p_den = DoubleDouble::mul_f64_add(
228            z,
229            ps_den,
230            DoubleDouble::from_bit_pair((0xbc759594c51ad8b7, 0x3feb84ebebabf275)),
231        );
232        p_den = DoubleDouble::mul_add_f64(z, p_den, f64::from_bits(0x3ff0000000000000));
233        p_den = DoubleDouble::from_exact_add(p_den.hi, p_den.lo);
234        result = DoubleDouble::div(p_num, p_den);
235        if y1.hi < dy.hi {
236            result = DoubleDouble::div(result, y1);
237        } else if y1.hi > dy.hi {
238            for _ in 0..n {
239                result = DoubleDouble::mult(result, dy);
240                dy = DoubleDouble::full_add_f64(dy, 1.0);
241            }
242        }
243    } else {
244        if x > 171.624e+0 {
245            return f64::INFINITY;
246        }
247        // Stirling's approximation of Log(Gamma) and then Exp[Log[Gamma]]
248        let y_recip = dy.recip();
249        let y_sqr = DoubleDouble::mult(y_recip, y_recip);
250        // Bernoulli coefficients generated by SageMath:
251        // var('x')
252        // def bernoulli_terms(x, N):
253        //     S = 0
254        //     for k in range(1, N+1):
255        //         B = bernoulli(2*k)
256        //         term = (B / (2*k*(2*k-1))) * x**((2*k-1))
257        //         S += term
258        //     return S
259        //
260        // terms = bernoulli_terms(x, 7)
261        let bernoulli_poly_s = f_polyeval6(
262            y_sqr.hi,
263            f64::from_bits(0xbf66c16c16c16c17),
264            f64::from_bits(0x3f4a01a01a01a01a),
265            f64::from_bits(0xbf43813813813814),
266            f64::from_bits(0x3f4b951e2b18ff23),
267            f64::from_bits(0xbf5f6ab0d9993c7d),
268            f64::from_bits(0x3f7a41a41a41a41a),
269        );
270        let bernoulli_poly = DoubleDouble::mul_f64_add(
271            y_sqr,
272            bernoulli_poly_s,
273            DoubleDouble::from_bit_pair((0x3c55555555555555, 0x3fb5555555555555)),
274        );
275        // Log[Gamma(x)] = x*log(x) - x + 1/2*Log(2*PI/x) + bernoulli_terms
276        const LOG2_PI_OVER_2: DoubleDouble =
277            DoubleDouble::from_bit_pair((0xbc865b5a1b7ff5df, 0x3fed67f1c864beb5));
278        let mut log_gamma = DoubleDouble::add(
279            DoubleDouble::mul_add(bernoulli_poly, y_recip, -dy),
280            LOG2_PI_OVER_2,
281        );
282        let dy_log = fast_log_dd(dy);
283        log_gamma = DoubleDouble::mul_add(
284            DoubleDouble::from_exact_add(dy_log.hi, dy_log.lo),
285            DoubleDouble::add_f64(dy, -0.5),
286            log_gamma,
287        );
288        let log_prod = log_gamma.to_f64();
289        if log_prod >= 690. {
290            // underflow/overflow case
291            log_gamma = DoubleDouble::quick_mult_f64(log_gamma, 0.5);
292            result = exp_dd_fast(log_gamma);
293            let exp_result = result;
294            result.hi *= parity;
295            result.lo *= parity;
296            if fact.lo != 0. && fact.hi != 0. {
297                // y / x = y / (z*z) = y / z * 1/z
298                result = DoubleDouble::from_exact_add(result.hi, result.lo);
299                result = DoubleDouble::div(fact, result);
300                result = DoubleDouble::div(result, exp_result);
301            } else {
302                result = DoubleDouble::quick_mult(result, exp_result);
303            }
304
305            let err = f_fmla(
306                result.hi.abs(),
307                f64::from_bits(0x3c20000000000000), // 2^-61
308                f64::from_bits(0x3bd0000000000000), // 2^-65
309            );
310            let ub = result.hi + (result.lo + err);
311            let lb = result.hi + (result.lo - err);
312            if ub == lb {
313                return result.to_f64();
314            }
315            return result.to_f64();
316        }
317        result = exp_dd_fast(log_gamma);
318    }
319
320    if fact.lo != 0. && fact.hi != 0. {
321        // y / x = y / (z*z) = y / z * 1/z
322        result = DoubleDouble::from_exact_add(result.hi, result.lo);
323        result = DoubleDouble::div(fact, result);
324    }
325    result.to_f64() * parity
326}
327
328#[cfg(test)]
329mod tests {
330    use super::*;
331
332    #[test]
333    fn test_tgamma() {
334        // assert_eq!(f_tgamma(6.757812502211891), 459.54924419209556);
335        assert_eq!(f_tgamma(-1.70000000042915), 2.513923520668069);
336        assert_eq!(f_tgamma(5.), 24.);
337        assert_eq!(f_tgamma(24.), 25852016738884980000000.);
338        assert_eq!(f_tgamma(6.4324324), 255.1369211339094);
339        assert_eq!(f_tgamma(f64::INFINITY), f64::INFINITY);
340        assert_eq!(f_tgamma(0.), f64::INFINITY);
341        assert_eq!(f_tgamma(-0.), f64::INFINITY);
342        assert!(f_tgamma(f64::NAN).is_nan());
343    }
344}
pxfm/gamma/tgamma.rs

pxfm/gamma/
tgamma.rs
