/rust/registry/src/index.crates.io-1949cf8c6b5b557f/pxfm-0.1.25/src/gamma/tgamma.rs

Source
/*
 * // Copyright (c) Radzivon Bartoshyk 8/2025. All rights reserved.
 * //
 * // Redistribution and use in source and binary forms, with or without modification,
 * // are permitted provided that the following conditions are met:
 * //
 * // 1.  Redistributions of source code must retain the above copyright notice, this
 * // list of conditions and the following disclaimer.
 * //
 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
 * // this list of conditions and the following disclaimer in the documentation
 * // and/or other materials provided with the distribution.
 * //
 * // 3.  Neither the name of the copyright holder nor the names of its
 * // contributors may be used to endorse or promote products derived from
 * // this software without specific prior written permission.
 * //
 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
use crate::common::{f_fmla, is_integer};
use crate::double_double::DoubleDouble;
use crate::logs::fast_log_dd;
use crate::polyeval::{f_polyeval6, f_polyeval7, f_polyeval8};
use crate::pow_exec::exp_dd_fast;
use crate::rounding::CpuFloor;
use crate::sincospi::f_fast_sinpi_dd;

/// Computes gamma(x)
///
/// ulp 1
pub fn f_tgamma(x: f64) -> f64 {
    let x_a = f64::from_bits(x.to_bits() & 0x7fff_ffff_ffff_ffff);

    if !x.is_normal() {
        if x == 0.0 {
            return f64::INFINITY;
        }

        if x.is_nan() {
            return x + x;
        }

        if x.is_infinite() {
            if x.is_sign_negative() {
                return f64::NAN;
            }
            return x;
        }
    }

    if x >= 171.624 {
        return f64::INFINITY;
    }

    if is_integer(x) {
        if x < 0. {
            return f64::NAN;
        }
        if x < 38. {
            let mut t = DoubleDouble::new(0., 1.);
            let k = x as i64;
            let mut x0 = 1.0;
            for _i in 1..k {
                t = DoubleDouble::quick_mult_f64(t, x0);
                t = DoubleDouble::from_exact_add(t.hi, t.lo);
                x0 += 1.0;
            }
            return t.to_f64();
        }
    }

    if x <= -184.0 {
        /* negative non-integer */
        /* For x <= -184, x non-integer, |gamma(x)| < 2^-1078.  */
        static SIGN: [f64; 2] = [
            f64::from_bits(0x0010000000000000),
            f64::from_bits(0x8010000000000000),
        ];
        let k = x as i64;
        return f64::from_bits(0x0010000000000000) * SIGN[((k & 1) != 0) as usize];
    }

    const EULER_DD: DoubleDouble =
        DoubleDouble::from_bit_pair((0xbc56cb90701fbfab, 0x3fe2788cfc6fb619));

    if x_a < 0.006 {
        if x_a.to_bits() < (0x71e0000000000000u64 >> 1) {
            // |x| < 0x1p-112
            return 1. / x;
        }
        if x_a < 2e-10 {
            // x is tiny then Gamma(x) = 1/x - euler
            let p = DoubleDouble::full_dd_sub(DoubleDouble::from_quick_recip(x), EULER_DD);
            return p.to_f64();
        } else if x_a < 2e-6 {
            // x is small then Gamma(x) = 1/x - euler + a2*x
            // a2 = 1/12 * (6 * euler^2 + pi^2)
            const A2: DoubleDouble =
                DoubleDouble::from_bit_pair((0x3c8dd92b465a8221, 0x3fefa658c23b1578));
            let rcp = DoubleDouble::from_quick_recip(x);
            let p = DoubleDouble::full_dd_add(DoubleDouble::mul_f64_add(A2, x, -EULER_DD), rcp);
            return p.to_f64();
        }

        // Laurent series of Gamma(x)
        const C: [(u64, u64); 8] = [
            (0x3c8dd92b465a8221, 0x3fefa658c23b1578),
            (0x3c53a4f483760950, 0xbfed0a118f324b63),
            (0x3c7fabe4f7369157, 0x3fef6a51055096b5),
            (0x3c8c9fc795fc6142, 0xbfef6c80ec38b67b),
            (0xbc5042339d62e721, 0x3fefc7e0a6eb310b),
            (0xbc86fd0d8bdc0c1e, 0xbfefdf3f157b7a39),
            (0xbc89b912df09395d, 0x3feff07b5a17ff6c),
            (0x3c4e626faf780ff9, 0xbfeff803d68a0bd4),
        ];
        let rcp = DoubleDouble::from_quick_recip(x);
        let mut p = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(C[7]),
            x,
            DoubleDouble::from_bit_pair(C[6]),
        );
        p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[5]));
        p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[4]));
        p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[3]));
        p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[2]));
        p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[1]));
        p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[0]));
        let z = DoubleDouble::mul_f64_add(p, x, DoubleDouble::full_dd_sub(rcp, EULER_DD));
        return z.to_f64();
    }

    let mut fact = DoubleDouble::new(0., 0.0f64);
    let mut parity = 1.0;
    const PI: DoubleDouble = DoubleDouble::from_bit_pair((0x3ca1a62633145c07, 0x400921fb54442d18));
    let mut dy = DoubleDouble::new(0., x);
    let mut result: DoubleDouble;

    // reflection
    if dy.hi < 0. {
        if dy.hi.cpu_floor() == dy.hi {
            return f64::NAN;
        }
        dy.hi = f64::from_bits(dy.hi.to_bits() & 0x7fff_ffff_ffff_ffff);
        let y1 = x_a.cpu_floor();
        let fraction = x_a - y1;
        if fraction != 0.0
        // is it an integer?
        {
            // is it odd or even?
            if y1 != (y1 * 0.5).cpu_floor() * 2.0 {
                parity = -1.0;
            }
            fact = DoubleDouble::div(-PI, f_fast_sinpi_dd(fraction));
            fact = DoubleDouble::from_exact_add(fact.hi, fact.lo);
            dy = DoubleDouble::from_full_exact_add(dy.hi, 1.0);
        }
    }

    if dy.hi < 12.0 {
        let y1 = dy;
        let z: DoubleDouble;
        let mut n = 0;
        // x is in (0.06, 1.0).
        if dy.hi < 1.0 {
            z = dy;
            dy = DoubleDouble::full_add_f64(dy, 1.0);
        } else
        // x is in [1.0, max].
        {
            n = dy.hi as i32 - 1;
            dy = DoubleDouble::full_add_f64(dy, -n as f64);
            z = DoubleDouble::full_add_f64(dy, -1.0);
        }

        // Gamma(x+1) on [1;2] generated by Wolfram Mathematica:
        // <<FunctionApproximations`
        // ClearAll["Global`*"]
        // f[x_]:=Gamma[x+1]
        // {err0, approx}=MiniMaxApproximation[f[z],{z,{0,1 },9,8},WorkingPrecision->90]
        // num=Numerator[approx][[1]];
        // den=Denominator[approx][[1]];
        // poly=den;
        // coeffs=CoefficientList[poly,z];
        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
        let ps_num = f_polyeval8(
            z.hi,
            f64::from_bits(0x3fb38b80568a42aa),
            f64::from_bits(0xbf8e7685b00d63a6),
            f64::from_bits(0xbf80629ed2c48f1a),
            f64::from_bits(0xbf6dfc4cdbcee96a),
            f64::from_bits(0xbf471816939dc42b),
            f64::from_bits(0xbf24bede7d8b3c20),
            f64::from_bits(0xbeef56936d891e42),
            f64::from_bits(0xbec3e2b405350813),
        );
        let mut p_num = DoubleDouble::mul_f64_add(
            z,
            ps_num,
            DoubleDouble::from_bit_pair((0xbc700f441aea2edb, 0x3fd218bdde7878b8)),
        );
        p_num = DoubleDouble::mul_add(
            z,
            p_num,
            DoubleDouble::from_bit_pair((0x3b7056a45a2fa50e, 0x3ff0000000000000)),
        );
        p_num = DoubleDouble::from_exact_add(p_num.hi, p_num.lo);

        let ps_den = f_polyeval7(
            z.hi,
            f64::from_bits(0xbfdaa4f09f0caab1),
            f64::from_bits(0xbfc960ba48423f9d),
            f64::from_bits(0x3fb6873b64e8ccd6),
            f64::from_bits(0x3f69ea1ca5b8a225),
            f64::from_bits(0xbf77b166f68a2e63),
            f64::from_bits(0x3f4fd1eff9193728),
            f64::from_bits(0xbf0c1a43f4985c97),
        );
        let mut p_den = DoubleDouble::mul_f64_add(
            z,
            ps_den,
            DoubleDouble::from_bit_pair((0xbc759594c51ad8b7, 0x3feb84ebebabf275)),
        );
        p_den = DoubleDouble::mul_add_f64(z, p_den, f64::from_bits(0x3ff0000000000000));
        p_den = DoubleDouble::from_exact_add(p_den.hi, p_den.lo);
        result = DoubleDouble::div(p_num, p_den);
        if y1.hi < dy.hi {
            result = DoubleDouble::div(result, y1);
        } else if y1.hi > dy.hi {
            for _ in 0..n {
                result = DoubleDouble::mult(result, dy);
                dy = DoubleDouble::full_add_f64(dy, 1.0);
            }
        }
    } else {
        if x > 171.624e+0 {
            return f64::INFINITY;
        }
        // Stirling's approximation of Log(Gamma) and then Exp[Log[Gamma]]
        let y_recip = dy.recip();
        let y_sqr = DoubleDouble::mult(y_recip, y_recip);
        // Bernoulli coefficients generated by SageMath:
        // var('x')
        // def bernoulli_terms(x, N):
        //     S = 0
        //     for k in range(1, N+1):
        //         B = bernoulli(2*k)
        //         term = (B / (2*k*(2*k-1))) * x**((2*k-1))
        //         S += term
        //     return S
        //
        // terms = bernoulli_terms(x, 7)
        let bernoulli_poly_s = f_polyeval6(
            y_sqr.hi,
            f64::from_bits(0xbf66c16c16c16c17),
            f64::from_bits(0x3f4a01a01a01a01a),
            f64::from_bits(0xbf43813813813814),
            f64::from_bits(0x3f4b951e2b18ff23),
            f64::from_bits(0xbf5f6ab0d9993c7d),
            f64::from_bits(0x3f7a41a41a41a41a),
        );
        let bernoulli_poly = DoubleDouble::mul_f64_add(
            y_sqr,
            bernoulli_poly_s,
            DoubleDouble::from_bit_pair((0x3c55555555555555, 0x3fb5555555555555)),
        );
        // Log[Gamma(x)] = x*log(x) - x + 1/2*Log(2*PI/x) + bernoulli_terms
        const LOG2_PI_OVER_2: DoubleDouble =
            DoubleDouble::from_bit_pair((0xbc865b5a1b7ff5df, 0x3fed67f1c864beb5));
        let mut log_gamma = DoubleDouble::add(
            DoubleDouble::mul_add(bernoulli_poly, y_recip, -dy),
            LOG2_PI_OVER_2,
        );
        let dy_log = fast_log_dd(dy);
        log_gamma = DoubleDouble::mul_add(
            DoubleDouble::from_exact_add(dy_log.hi, dy_log.lo),
            DoubleDouble::add_f64(dy, -0.5),
            log_gamma,
        );
        let log_prod = log_gamma.to_f64();
        if log_prod >= 690. {
            // underflow/overflow case
            log_gamma = DoubleDouble::quick_mult_f64(log_gamma, 0.5);
            result = exp_dd_fast(log_gamma);
            let exp_result = result;
            result.hi *= parity;
            result.lo *= parity;
            if fact.lo != 0. && fact.hi != 0. {
                // y / x = y / (z*z) = y / z * 1/z
                result = DoubleDouble::from_exact_add(result.hi, result.lo);
                result = DoubleDouble::div(fact, result);
                result = DoubleDouble::div(result, exp_result);
            } else {
                result = DoubleDouble::quick_mult(result, exp_result);
            }

            let err = f_fmla(
                result.hi.abs(),
                f64::from_bits(0x3c20000000000000), // 2^-61
                f64::from_bits(0x3bd0000000000000), // 2^-65
            );
            let ub = result.hi + (result.lo + err);
            let lb = result.hi + (result.lo - err);
            if ub == lb {
                return result.to_f64();
            }
            return result.to_f64();
        }
        result = exp_dd_fast(log_gamma);
    }

    if fact.lo != 0. && fact.hi != 0. {
        // y / x = y / (z*z) = y / z * 1/z
        result = DoubleDouble::from_exact_add(result.hi, result.lo);
        result = DoubleDouble::div(fact, result);
    }
    result.to_f64() * parity
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_tgamma() {
        // assert_eq!(f_tgamma(6.757812502211891), 459.54924419209556);
        assert_eq!(f_tgamma(-1.70000000042915), 2.513923520668069);
        assert_eq!(f_tgamma(5.), 24.);
        assert_eq!(f_tgamma(24.), 25852016738884980000000.);
        assert_eq!(f_tgamma(6.4324324), 255.1369211339094);
        assert_eq!(f_tgamma(f64::INFINITY), f64::INFINITY);
        assert_eq!(f_tgamma(0.), f64::INFINITY);
        assert_eq!(f_tgamma(-0.), f64::INFINITY);
        assert!(f_tgamma(f64::NAN).is_nan());
    }
}

Coverage Report

Created: 2025-11-24 07:30

Line	Count	Source
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		*/
29		use crate::common::{f_fmla, is_integer};
30		use crate::double_double::DoubleDouble;
31		use crate::logs::fast_log_dd;
32		use crate::polyeval::{f_polyeval6, f_polyeval7, f_polyeval8};
33		use crate::pow_exec::exp_dd_fast;
34		use crate::rounding::CpuFloor;
35		use crate::sincospi::f_fast_sinpi_dd;
36
37		/// Computes gamma(x)
38		///
39		/// ulp 1
40	0	pub fn f_tgamma(x: f64) -> f64 {
41	0	let x_a = f64::from_bits(x.to_bits() & 0x7fff_ffff_ffff_ffff);
42
43	0	if !x.is_normal() {
44	0	if x == 0.0 {
45	0	return f64::INFINITY;
46	0	}
47
48	0	if x.is_nan() {
49	0	return x + x;
50	0	}
51
52	0	if x.is_infinite() {
53	0	if x.is_sign_negative() {
54	0	return f64::NAN;
55	0	}
56	0	return x;
57	0	}
58	0	}
59
60	0	if x >= 171.624 {
61	0	return f64::INFINITY;
62	0	}
63
64	0	if is_integer(x) {
65	0	if x < 0. {
66	0	return f64::NAN;
67	0	}
68	0	if x < 38. {
69	0	let mut t = DoubleDouble::new(0., 1.);
70	0	let k = x as i64;
71	0	let mut x0 = 1.0;
72	0	for _i in 1..k {
73	0	t = DoubleDouble::quick_mult_f64(t, x0);
74	0	t = DoubleDouble::from_exact_add(t.hi, t.lo);
75	0	x0 += 1.0;
76	0	}
77	0	return t.to_f64();
78	0	}
79	0	}
80
81	0	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	0	let k = x as i64;
89	0	return f64::from_bits(0x0010000000000000) * SIGN[((k & 1) != 0) as usize];
90	0	}
91
92		const EULER_DD: DoubleDouble =
93		DoubleDouble::from_bit_pair((0xbc56cb90701fbfab, 0x3fe2788cfc6fb619));
94
95	0	if x_a < 0.006 {
96	0	if x_a.to_bits() < (0x71e0000000000000u64 >> 1) {
97		// \|x\| < 0x1p-112
98	0	return 1. / x;
99	0	}
100	0	if x_a < 2e-10 {
101		// x is tiny then Gamma(x) = 1/x - euler
102	0	let p = DoubleDouble::full_dd_sub(DoubleDouble::from_quick_recip(x), EULER_DD);
103	0	return p.to_f64();
104	0	} 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	0	let rcp = DoubleDouble::from_quick_recip(x);
110	0	let p = DoubleDouble::full_dd_add(DoubleDouble::mul_f64_add(A2, x, -EULER_DD), rcp);
111	0	return p.to_f64();
112	0	}
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	0	let rcp = DoubleDouble::from_quick_recip(x);
126	0	let mut p = DoubleDouble::mul_f64_add(
127	0	DoubleDouble::from_bit_pair(C[7]),
128	0	x,
129	0	DoubleDouble::from_bit_pair(C[6]),
130		);
131	0	p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[5]));
132	0	p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[4]));
133	0	p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[3]));
134	0	p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[2]));
135	0	p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[1]));
136	0	p = DoubleDouble::mul_f64_add(p, x, DoubleDouble::from_bit_pair(C[0]));
137	0	let z = DoubleDouble::mul_f64_add(p, x, DoubleDouble::full_dd_sub(rcp, EULER_DD));
138	0	return z.to_f64();
139	0	}
140
141	0	let mut fact = DoubleDouble::new(0., 0.0f64);
142	0	let mut parity = 1.0;
143		const PI: DoubleDouble = DoubleDouble::from_bit_pair((0x3ca1a62633145c07, 0x400921fb54442d18));
144	0	let mut dy = DoubleDouble::new(0., x);
145		let mut result: DoubleDouble;
146
147		// reflection
148	0	if dy.hi < 0. {
149	0	if dy.hi.cpu_floor() == dy.hi {
150	0	return f64::NAN;
151	0	}
152	0	dy.hi = f64::from_bits(dy.hi.to_bits() & 0x7fff_ffff_ffff_ffff);
153	0	let y1 = x_a.cpu_floor();
154	0	let fraction = x_a - y1;
155	0	if fraction != 0.0
156		// is it an integer?
157		{
158		// is it odd or even?
159	0	if y1 != (y1 * 0.5).cpu_floor() * 2.0 {
160	0	parity = -1.0;
161	0	}
162	0	fact = DoubleDouble::div(-PI, f_fast_sinpi_dd(fraction));
163	0	fact = DoubleDouble::from_exact_add(fact.hi, fact.lo);
164	0	dy = DoubleDouble::from_full_exact_add(dy.hi, 1.0);
165	0	}
166	0	}
167
168	0	if dy.hi < 12.0 {
169	0	let y1 = dy;
170		let z: DoubleDouble;
171	0	let mut n = 0;
172		// x is in (0.06, 1.0).
173	0	if dy.hi < 1.0 {
174	0	z = dy;
175	0	dy = DoubleDouble::full_add_f64(dy, 1.0);
176	0	} else
177		// x is in [1.0, max].
178	0	{
179	0	n = dy.hi as i32 - 1;
180	0	dy = DoubleDouble::full_add_f64(dy, -n as f64);
181	0	z = DoubleDouble::full_add_f64(dy, -1.0);
182	0	}
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	0	let ps_num = f_polyeval8(
195	0	z.hi,
196	0	f64::from_bits(0x3fb38b80568a42aa),
197	0	f64::from_bits(0xbf8e7685b00d63a6),
198	0	f64::from_bits(0xbf80629ed2c48f1a),
199	0	f64::from_bits(0xbf6dfc4cdbcee96a),
200	0	f64::from_bits(0xbf471816939dc42b),
201	0	f64::from_bits(0xbf24bede7d8b3c20),
202	0	f64::from_bits(0xbeef56936d891e42),
203	0	f64::from_bits(0xbec3e2b405350813),
204		);
205	0	let mut p_num = DoubleDouble::mul_f64_add(
206	0	z,
207	0	ps_num,
208	0	DoubleDouble::from_bit_pair((0xbc700f441aea2edb, 0x3fd218bdde7878b8)),
209		);
210	0	p_num = DoubleDouble::mul_add(
211	0	z,
212	0	p_num,
213	0	DoubleDouble::from_bit_pair((0x3b7056a45a2fa50e, 0x3ff0000000000000)),
214	0	);
215	0	p_num = DoubleDouble::from_exact_add(p_num.hi, p_num.lo);
216
217	0	let ps_den = f_polyeval7(
218	0	z.hi,
219	0	f64::from_bits(0xbfdaa4f09f0caab1),
220	0	f64::from_bits(0xbfc960ba48423f9d),
221	0	f64::from_bits(0x3fb6873b64e8ccd6),
222	0	f64::from_bits(0x3f69ea1ca5b8a225),
223	0	f64::from_bits(0xbf77b166f68a2e63),
224	0	f64::from_bits(0x3f4fd1eff9193728),
225	0	f64::from_bits(0xbf0c1a43f4985c97),
226		);
227	0	let mut p_den = DoubleDouble::mul_f64_add(
228	0	z,
229	0	ps_den,
230	0	DoubleDouble::from_bit_pair((0xbc759594c51ad8b7, 0x3feb84ebebabf275)),
231		);
232	0	p_den = DoubleDouble::mul_add_f64(z, p_den, f64::from_bits(0x3ff0000000000000));
233	0	p_den = DoubleDouble::from_exact_add(p_den.hi, p_den.lo);
234	0	result = DoubleDouble::div(p_num, p_den);
235	0	if y1.hi < dy.hi {
236	0	result = DoubleDouble::div(result, y1);
237	0	} else if y1.hi > dy.hi {
238	0	for _ in 0..n {
239	0	result = DoubleDouble::mult(result, dy);
240	0	dy = DoubleDouble::full_add_f64(dy, 1.0);
241	0	}
242	0	}
243		} else {
244	0	if x > 171.624e+0 {
245	0	return f64::INFINITY;
246	0	}
247		// Stirling's approximation of Log(Gamma) and then Exp[Log[Gamma]]
248	0	let y_recip = dy.recip();
249	0	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 / (2k(2k-1))) x*((2k-1))
257		// S += term
258		// return S
259		//
260		// terms = bernoulli_terms(x, 7)
261	0	let bernoulli_poly_s = f_polyeval6(
262	0	y_sqr.hi,
263	0	f64::from_bits(0xbf66c16c16c16c17),
264	0	f64::from_bits(0x3f4a01a01a01a01a),
265	0	f64::from_bits(0xbf43813813813814),
266	0	f64::from_bits(0x3f4b951e2b18ff23),
267	0	f64::from_bits(0xbf5f6ab0d9993c7d),
268	0	f64::from_bits(0x3f7a41a41a41a41a),
269		);
270	0	let bernoulli_poly = DoubleDouble::mul_f64_add(
271	0	y_sqr,
272	0	bernoulli_poly_s,
273	0	DoubleDouble::from_bit_pair((0x3c55555555555555, 0x3fb5555555555555)),
274		);
275		// Log[Gamma(x)] = xlog(x) - x + 1/2Log(2*PI/x) + bernoulli_terms
276		const LOG2_PI_OVER_2: DoubleDouble =
277		DoubleDouble::from_bit_pair((0xbc865b5a1b7ff5df, 0x3fed67f1c864beb5));
278	0	let mut log_gamma = DoubleDouble::add(
279	0	DoubleDouble::mul_add(bernoulli_poly, y_recip, -dy),
280		LOG2_PI_OVER_2,
281		);
282	0	let dy_log = fast_log_dd(dy);
283	0	log_gamma = DoubleDouble::mul_add(
284	0	DoubleDouble::from_exact_add(dy_log.hi, dy_log.lo),
285	0	DoubleDouble::add_f64(dy, -0.5),
286	0	log_gamma,
287	0	);
288	0	let log_prod = log_gamma.to_f64();
289	0	if log_prod >= 690. {
290		// underflow/overflow case
291	0	log_gamma = DoubleDouble::quick_mult_f64(log_gamma, 0.5);
292	0	result = exp_dd_fast(log_gamma);
293	0	let exp_result = result;
294	0	result.hi *= parity;
295	0	result.lo *= parity;
296	0	if fact.lo != 0. && fact.hi != 0. {
297	0	// y / x = y / (zz) = y / z 1/z
298	0	result = DoubleDouble::from_exact_add(result.hi, result.lo);
299	0	result = DoubleDouble::div(fact, result);
300	0	result = DoubleDouble::div(result, exp_result);
301	0	} else {
302	0	result = DoubleDouble::quick_mult(result, exp_result);
303	0	}
304
305	0	let err = f_fmla(
306	0	result.hi.abs(),
307	0	f64::from_bits(0x3c20000000000000), // 2^-61
308	0	f64::from_bits(0x3bd0000000000000), // 2^-65
309		);
310	0	let ub = result.hi + (result.lo + err);
311	0	let lb = result.hi + (result.lo - err);
312	0	if ub == lb {
313	0	return result.to_f64();
314	0	}
315	0	return result.to_f64();
316	0	}
317	0	result = exp_dd_fast(log_gamma);
318		}
319
320	0	if fact.lo != 0. && fact.hi != 0. {
321	0	// y / x = y / (zz) = y / z 1/z
322	0	result = DoubleDouble::from_exact_add(result.hi, result.lo);
323	0	result = DoubleDouble::div(fact, result);
324	0	}
325	0	result.to_f64() * parity
326	0	}
327
328		#[cfg(test)]
329		mod 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		}