pxfm/gamma/
lgamma_rf.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 9/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_integerf, is_odd_integerf};
30use crate::floor::FloorFinite;
31use crate::logs::simple_fast_log;
32use crate::polyeval::{f_estrin_polyeval4, f_polyeval5, f_polyeval7, f_polyeval8};
33use crate::sin_cosf::fast_sinpif;
34
35#[inline]
36pub(crate) fn lgamma_coref(x: f32) -> (f64, i32) {
37    let ax = f32::from_bits(x.to_bits() & 0x7fff_ffff);
38    let dx = ax as f64;
39
40    let is_positive = x.is_sign_positive();
41    let mut sum_parity = 1f64;
42
43    let mut f_res = 0f64;
44
45    let mut signgam: i32 = 1i32;
46
47    // For negative x, since (G is gamma function)
48    // -x*G(-x)*G(x) = pi/sin(pi*x),
49    // we have
50    // G(x) = pi/(sin(pi*x)*(-x)*G(-x))
51    // since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
52    // Hence, for x<0, signgam = sign(sin(pi*x)) and
53    // lgamma(x) = log(|Gamma(x)|)
54    //  = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
55    if !is_positive {
56        let y1 = ax.floor_finite();
57        let fraction = ax - y1; // excess over the boundary
58
59        let a = fast_sinpif(fraction);
60
61        sum_parity = -1.;
62        const LOG_PI: f64 = f64::from_bits(0x3ff250d048e7a1bd);
63        f_res = LOG_PI - simple_fast_log(a * dx);
64        // gamma(x) is negative in (-2n-1,-2n), thus when fx is odd
65        let is_odd = (!is_odd_integerf(y1)) as i32;
66        signgam = 1 - (is_odd << 1);
67    }
68
69    if ax <= 0.007 {
70        // Generated by Wolfram Mathematica:
71        // <<FunctionApproximations`
72        // ClearAll["Global`*"]
73        // f[x_]:=LogGamma[x+1]
74        // {err0,approx}=MiniMaxApproximation[f[x],{x,{2^-32,0.07},3,3},WorkingPrecision->90]
75        // num=Numerator[approx][[1]];
76        // den=Denominator[approx][[1]];
77        // poly=num;
78        // coeffs=CoefficientList[poly,x];
79        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
80        // poly=den;
81        // coeffs=CoefficientList[poly,x];
82        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
83        let num = f_estrin_polyeval4(
84            dx,
85            f64::from_bits(0xbb49bbc47e595350),
86            f64::from_bits(0xbfe2788cfc6fb2d3),
87            f64::from_bits(0x3fc829299b97c525),
88            f64::from_bits(0x3fd8e2e8aab1247c),
89        );
90        let den = f_estrin_polyeval4(
91            dx,
92            f64::from_bits(0x3ff0000000000000),
93            f64::from_bits(0x3ff190e5bf5a8162),
94            f64::from_bits(0x3fc927632366c502),
95            f64::from_bits(0xbf8b4ea7bd1b808e),
96        );
97        f_res = f_fmla(num / den - simple_fast_log(dx), sum_parity, f_res);
98    } else if ax < 1. {
99        // Poly for loggamma(x + 1) - log(x) generated by Wolfram Mathematica
100        // <<FunctionApproximations`
101        // ClearAll["Global`*"]
102        // f[x_]:=LogGamma[x+1]
103        // {err0, approx}=MiniMaxApproximation[f[z],{z,{0.06,0.9999999},6,6},WorkingPrecision->90]
104        // num=Numerator[approx][[1]];
105        // den=Denominator[approx][[1]];
106        // poly=den;
107        // coeffs=CoefficientList[poly,z];
108        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
109        let num = f_polyeval7(
110            dx,
111            f64::from_bits(0xbd24e4cf78c8818a),
112            f64::from_bits(0xbfe2788cfc6f64e6),
113            f64::from_bits(0xbfe1ea632fe853b2),
114            f64::from_bits(0x3fd988d2daad3806),
115            f64::from_bits(0x3fe1f4870eaafdf4),
116            f64::from_bits(0x3fc51e12d5b97330),
117            f64::from_bits(0x3f889ebeb9f0c8ff),
118        );
119        let den = f_polyeval7(
120            dx,
121            f64::from_bits(0x3ff0000000000000),
122            f64::from_bits(0x4003289872066f33),
123            f64::from_bits(0x4000373d88c32fe5),
124            f64::from_bits(0x3fe71e95bc874164),
125            f64::from_bits(0x3fb9936a0e008be7),
126            f64::from_bits(0x3f6cba20a1988a60),
127            f64::from_bits(0xbef3ca884ba4237a),
128        );
129        f_res = f_fmla(num / den - simple_fast_log(dx), sum_parity, f_res);
130    } else if ax <= 4. {
131        // Poly for loggamma(x + 1) - log(x) generated by Wolfram Mathematica
132        // <<FunctionApproximations`
133        // ClearAll["Global`*"]
134        // f[x_]:=LogGamma[x+1]
135        // {err0, approx}=MiniMaxApproximation[f[z],{z,{1.0000001,4},7,7},WorkingPrecision->90]
136        // num=Numerator[approx][[1]];
137        // den=Denominator[approx][[1]];
138        // poly=den;
139        // coeffs=CoefficientList[poly,z];
140        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
141        let num = f_polyeval8(
142            dx,
143            f64::from_bits(0xbe630847302a205f),
144            f64::from_bits(0xbfe2788c24a7deff),
145            f64::from_bits(0xbfdd0f61b8171907),
146            f64::from_bits(0x3fdab4ec27801e42),
147            f64::from_bits(0x3fde11fd49427a3c),
148            f64::from_bits(0x3fc0dbce39c660da),
149            f64::from_bits(0x3f88e0b4cd3e97f7),
150            f64::from_bits(0x3f328fcd9e9ee94e),
151        );
152        let den = f_polyeval8(
153            dx,
154            f64::from_bits(0x3ff0000000000000),
155            f64::from_bits(0x4001b135a31bffae),
156            f64::from_bits(0x3ffbbecb5e39f4ea),
157            f64::from_bits(0x3fe2e4da095d10bd),
158            f64::from_bits(0x3fb63cb7548d0d30),
159            f64::from_bits(0x3f73a67613163399),
160            f64::from_bits(0x3f10edd4e2b80e6f),
161            f64::from_bits(0xbe77f4e4068f9d32),
162        );
163        f_res = f_fmla(num / den - simple_fast_log(dx), sum_parity, f_res);
164    } else if ax <= 12. {
165        // Poly for loggamma(x) generated by Wolfram Mathematica
166        // <<FunctionApproximations`
167        // ClearAll["Global`*"]
168        // f[x_]:=LogGamma[x]
169        // {err0, approx}=MiniMaxApproximation[f[z],{z,{4,12},7,7},WorkingPrecision->90]
170        // num=Numerator[approx][[1]];
171        // den=Denominator[approx][[1]];
172        // poly=den;
173        // coeffs=CoefficientList[poly,z];
174        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
175        let num = f_polyeval8(
176            dx,
177            f64::from_bits(0x40086069fcc21690),
178            f64::from_bits(0x4008c8ef8602e78d),
179            f64::from_bits(0xc0186323149757ae),
180            f64::from_bits(0xbff6a9ded42d3f31),
181            f64::from_bits(0x3ff1adb5566cd807),
182            f64::from_bits(0x3fcff023390dead6),
183            f64::from_bits(0x3f8977a5454ec393),
184            f64::from_bits(0x3f21361088f694dc),
185        );
186        let den = f_polyeval8(
187            dx,
188            f64::from_bits(0x3ff0000000000000),
189            f64::from_bits(0x4016096b397dbe73),
190            f64::from_bits(0x401478926593a395),
191            f64::from_bits(0x3ff6a931f9e31511),
192            f64::from_bits(0x3fc0fbe57bf263a4),
193            f64::from_bits(0x3f7023f51ba39a98),
194            f64::from_bits(0x3efabd75ebfbde40),
195            f64::from_bits(0xbe4b2c9e2e4e7daa),
196        );
197        f_res = f_fmla(num / den, sum_parity, f_res);
198    } else {
199        // Stirling's approximation of Log(Gamma) and then Exp[Log[Gamma]]
200        let y_recip = 1. / dx;
201        let y_sqr = y_recip * y_recip;
202        // Bernoulli coefficients generated by SageMath:
203        // var('x')
204        // def bernoulli_terms(x, N):
205        //     S = 0
206        //     for k in range(1, N+1):
207        //         B = bernoulli(2*k)
208        //         term = (B / (2*k*(2*k-1))) * x**((2*k-1))
209        //         S += term
210        //     return S
211        //
212        // terms = bernoulli_terms(x, 5)
213        let bernoulli_poly = f_polyeval5(
214            y_sqr,
215            f64::from_bits(0x3fb5555555555555),
216            f64::from_bits(0xbf66c16c16c16c17),
217            f64::from_bits(0x3f4a01a01a01a01a),
218            f64::from_bits(0xbf43813813813814),
219            f64::from_bits(0x3f4b951e2b18ff23),
220        );
221        // Log[Gamma(x)] = x*log(x) - x + 1/2*Log(2*PI/x) + bernoulli_terms
222        const LOG2_PI_OVER_2: f64 = f64::from_bits(0x3fed67f1c864beb5);
223        let mut log_gamma = f_fmla(bernoulli_poly, y_recip, -dx) + LOG2_PI_OVER_2;
224        log_gamma = f_fmla(simple_fast_log(dx), dx - 0.5, log_gamma);
225        f_res = f_fmla(log_gamma, sum_parity, f_res);
226    }
227
228    (f_res, signgam)
229}
230
231/// Computes log(gamma(x))
232///
233/// Returns gamma value + sign.
234///
235/// ulp 0.5
236pub fn f_lgamma_rf(x: f32) -> (f32, i32) {
237    // filter out exceptional cases
238    let xb = x.to_bits().wrapping_shl(1);
239    if xb >= 0xffu32 << 24 || xb == 0 {
240        if x.is_infinite() {
241            return (f32::INFINITY, 1);
242        }
243        if x.is_nan() {
244            return (f32::NAN, 1);
245        }
246        if xb == 0 {
247            return (f32::INFINITY, 1 - 2 * (x.to_bits() >> 31) as i32);
248        }
249    }
250
251    if is_integerf(x) {
252        if x == 2. || x == 1. {
253            return (0., 1);
254        }
255        if x.is_sign_negative() {
256            let is_odd = (!is_odd_integerf(x)) as i32;
257            return (f32::INFINITY, 1 - (is_odd << 1));
258        }
259    }
260    let c_gamma = lgamma_coref(x);
261    let fx = c_gamma.0 as f32;
262    (fx, c_gamma.1)
263}
264
265#[cfg(test)]
266mod tests {
267    use crate::f_lgamma_rf;
268
269    #[test]
270    fn test_lgamma_rf() {
271        assert_eq!(f_lgamma_rf(f32::NEG_INFINITY), (f32::INFINITY, 1));
272        assert_eq!(f_lgamma_rf(f32::INFINITY), (f32::INFINITY, 1));
273        assert_eq!(f_lgamma_rf(0.), (f32::INFINITY, 1));
274        assert_eq!(f_lgamma_rf(-0.), (f32::INFINITY, -1));
275        assert_eq!(f_lgamma_rf(5.), (3.1780539, 1));
276        assert_eq!(f_lgamma_rf(-4.5), (-2.8130841, -1));
277        assert_eq!(f_lgamma_rf(-2.0015738), (5.7596655, -1));
278    }
279}
pxfm/gamma/lgamma_rf.rs

pxfm/gamma/
lgamma_rf.rs
