pxfm/gamma/tgammaf.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::f_exp;
31use crate::floor::FloorFinite;
32use crate::logs::simple_fast_log;
33use crate::polyeval::{f_polyeval5, f_polyeval18};
34use crate::sin_cosf::fast_sinpif;
35
36/// True gamma function
37///
38/// ulp 0.5
39pub fn f_tgammaf(x: f32) -> f32 {
40 let xb = x.to_bits();
41 // filter out exceptional cases
42 if xb >= 0xffu32 << 23 || xb == 0 {
43 if xb.wrapping_shl(1) == 0 {
44 return f32::INFINITY;
45 }
46
47 if x.is_nan() {
48 return x + x;
49 }
50
51 if x.is_infinite() {
52 if x.is_sign_negative() {
53 return f32::NAN;
54 }
55 return x;
56 }
57 }
58
59 let x_a = f32::from_bits(x.to_bits() & 0x7fff_ffff);
60 if x.is_sign_positive() && x_a.to_bits() >= 0x420c2910u32 {
61 // x >= 35.0401
62 return f32::INFINITY;
63 } else if x.is_sign_negative() && x_a.to_bits() >= 0x421a67dau32 {
64 // x <= -38.601418
65 if x == x.floor_finite() {
66 // integer where x < 0
67 return f32::NAN;
68 }
69 return -0.0;
70 }
71
72 if x_a < 2e-8 {
73 // x is tiny then Gamma(x) = 1/x - euler
74 const EULER: f64 = f64::from_bits(0x3fe2788cfc6fb619);
75 let dx = x as f64;
76 let p = 1. / dx - EULER;
77 return p as f32;
78 } else if x_a < 2e-6 {
79 // x is small then Gamma(x) = 1/x - euler + a2*x
80 // a2 = 1/12 * (6 * euler^2 + pi^2)
81 const EULER: f64 = f64::from_bits(0x3fe2788cfc6fb619);
82 const A2: f64 = f64::from_bits(0x3fefa658c23b1578);
83 let dx = x as f64;
84 let p = 1. / dx + f_fmla(dx, A2, -EULER);
85 return p as f32;
86 }
87
88 let mut fact = 1.0f64;
89 let mut parity = 1.0;
90 const PI: f64 = f64::from_bits(0x400921fb54442d18);
91 let mut dy = x as f64;
92 let mut result;
93
94 // reflection
95 if dy < 0. {
96 if is_integer(dy) {
97 return f32::NAN;
98 }
99 dy = f64::from_bits(dy.to_bits() & 0x7fff_ffff_ffff_ffff);
100 // x < ulp(x)
101 if x_a < f32::EPSILON {
102 return (1. / x as f64) as f32;
103 }
104 let y1 = x_a.floor_finite();
105 let fraction = x_a - y1;
106 if fraction != 0.0
107 // is it an integer?
108 {
109 // is it odd or even?
110 if y1 != (y1 * 0.5).floor_finite() * 2.0 {
111 parity = -1.0;
112 }
113 fact = -PI / fast_sinpif(fraction);
114 dy += 1.0;
115 }
116 }
117
118 if dy < 12.0 {
119 let y1 = dy;
120 let z;
121 let mut n = 0;
122 // x is in (eps, 1.0).
123 if dy < 1.0 {
124 z = dy;
125 dy += 1.0;
126 } else
127 // x is in [1.0, max].
128 {
129 n = dy as i32 - 1;
130 dy -= n as f64;
131 z = dy - 1.0;
132 }
133
134 // Gamma(x+1) on [1;2] generated by Wolfram Mathematica:
135 // <<FunctionApproximations`
136 // ClearAll["Global`*"]
137 // f[x_]:=Gamma[x+1]
138 // {err0, approx}=MiniMaxApproximation[f[z],{z,{0,1},17,0},WorkingPrecision->90]
139 // num=Numerator[approx][[1]];
140 // den=Denominator[approx][[1]];
141 // poly=num;
142 // coeffs=CoefficientList[poly,z];
143 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]], {50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
144 result = f_polyeval18(
145 z,
146 f64::from_bits(0x3fefffffffffff19),
147 f64::from_bits(0xbfe2788cfc6d59a6),
148 f64::from_bits(0x3fefa658c133abd2),
149 f64::from_bits(0xbfed0a116184e1d0),
150 f64::from_bits(0x3fef6a4cd2befb54),
151 f64::from_bits(0xbfef6c4479fe9aca),
152 f64::from_bits(0x3fefc59bfa97e9bf),
153 f64::from_bits(0xbfefcfde03bfb0d9),
154 f64::from_bits(0x3fefa3ee5eab6681),
155 f64::from_bits(0xbfeed8130a67dd46),
156 f64::from_bits(0x3fecb6e603fdb5ed),
157 f64::from_bits(0xbfe8801e901b4c10),
158 f64::from_bits(0x3fe2356219082d3e),
159 f64::from_bits(0xbfd64ad556b6062a),
160 f64::from_bits(0x3fc5219e235288f4),
161 f64::from_bits(0xbfacad444cb0f4ec),
162 f64::from_bits(0x3f888ea3bfa9155a),
163 f64::from_bits(0xbf53d1776abc9eaa),
164 );
165 if y1 < dy {
166 result /= y1;
167 } else if y1 > dy {
168 for _ in 0..n {
169 result *= dy;
170 dy += 1.0;
171 }
172 }
173 } else {
174 // Stirling's approximation of Log(Gamma) and then Exp[Log[Gamma]]
175 let y_recip = 1. / dy;
176 let y_sqr = y_recip * y_recip;
177 // Bernoulli coefficients generated by SageMath:
178 // var('x')
179 // def bernoulli_terms(x, N):
180 // S = 0
181 // for k in range(1, N+1):
182 // B = bernoulli(2*k)
183 // term = (B / (2*k*(2*k-1))) * x**((2*k-1))
184 // S += term
185 // return S
186 //
187 // terms = bernoulli_terms(x, 5)
188 let bernoulli_poly = f_polyeval5(
189 y_sqr,
190 f64::from_bits(0x3fb5555555555555),
191 f64::from_bits(0xbf66c16c16c16c17),
192 f64::from_bits(0x3f4a01a01a01a01a),
193 f64::from_bits(0xbf43813813813814),
194 f64::from_bits(0x3f4b951e2b18ff23),
195 );
196 // Log[Gamma(x)] = x*log(x) - x + 1/2*Log(2*PI/x) + bernoulli_terms
197 const LOG2_PI_OVER_2: f64 = f64::from_bits(0x3fed67f1c864beb5);
198 let mut log_gamma = f_fmla(bernoulli_poly, y_recip, -dy) + LOG2_PI_OVER_2;
199 log_gamma = f_fmla(simple_fast_log(dy), dy - 0.5, log_gamma);
200 result = f_exp(log_gamma);
201 }
202
203 result *= parity;
204 if fact != 1.0 {
205 result = fact / result;
206 }
207
208 result as f32
209}
210
211#[cfg(test)]
212mod tests {
213 use crate::f_tgammaf;
214
215 #[test]
216 fn test_gammaf() {
217 assert!(f_tgammaf(-3.).is_nan());
218 assert_eq!(f_tgammaf(12.58038769), 167149000.0);
219 assert_eq!(f_tgammaf(4.566854e-7), 2189690.8);
220 assert_eq!(f_tgammaf(4.9587783e-8), 20166258.0);
221 assert_eq!(f_tgammaf(-46.1837), -0.0);
222 assert_eq!(f_tgammaf(1.1502532e-19), 8.6937384e18);
223 assert_eq!(f_tgammaf(0.04038769), 24.221333);
224 assert_eq!(f_tgammaf(2.58038769), 1.4088479);
225 assert_eq!(f_tgammaf(-1.3559477), 2.892815);
226 assert_eq!(f_tgammaf(-2.3559477), -1.2278775);
227 assert_eq!(f_tgammaf(0.0), f32::INFINITY);
228 assert_eq!(f_tgammaf(-0.0), f32::INFINITY);
229 assert_eq!(f_tgammaf(f32::INFINITY), f32::INFINITY);
230 assert!(f_tgammaf(f32::NEG_INFINITY).is_nan());
231 assert!(f_tgammaf(f32::NAN).is_nan());
232 }
233}