1/* Implementation of gamma function according to ISO C.
2 Copyright (C) 1997-2022 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4
5 The GNU C Library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Lesser General Public
7 License as published by the Free Software Foundation; either
8 version 2.1 of the License, or (at your option) any later version.
9
10 The GNU C Library is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public
16 License along with the GNU C Library; if not, see
17 <https://www.gnu.org/licenses/>. */
18
19#include <math.h>
20#include <math-narrow-eval.h>
21#include <math_private.h>
22#include <fenv_private.h>
23#include <math-underflow.h>
24#include <float.h>
25#include <libm-alias-finite.h>
26
27/* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
28 approximation to gamma function. */
29
30static const float gamma_coeff[] =
31 {
32 0x1.555556p-4f,
33 -0xb.60b61p-12f,
34 0x3.403404p-12f,
35 };
36
37#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0]))
38
39/* Return gamma (X), for positive X less than 42, in the form R *
40 2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to
41 avoid overflow or underflow in intermediate calculations. */
42
43static float
44gammaf_positive (float x, int *exp2_adj)
45{
46 int local_signgam;
47 if (x < 0.5f)
48 {
49 *exp2_adj = 0;
50 return __ieee754_expf (__ieee754_lgammaf_r (x + 1, &local_signgam)) / x;
51 }
52 else if (x <= 1.5f)
53 {
54 *exp2_adj = 0;
55 return __ieee754_expf (__ieee754_lgammaf_r (x, &local_signgam));
56 }
57 else if (x < 2.5f)
58 {
59 *exp2_adj = 0;
60 float x_adj = x - 1;
61 return (__ieee754_expf (__ieee754_lgammaf_r (x_adj, &local_signgam))
62 * x_adj);
63 }
64 else
65 {
66 float eps = 0;
67 float x_eps = 0;
68 float x_adj = x;
69 float prod = 1;
70 if (x < 4.0f)
71 {
72 /* Adjust into the range for applying Stirling's
73 approximation. */
74 float n = ceilf (4.0f - x);
75 x_adj = math_narrow_eval (x + n);
76 x_eps = (x - (x_adj - n));
77 prod = __gamma_productf (x_adj - n, x_eps, n, &eps);
78 }
79 /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)).
80 Compute gamma (X_ADJ + X_EPS) using Stirling's approximation,
81 starting by computing pow (X_ADJ, X_ADJ) with a power of 2
82 factored out. */
83 float exp_adj = -eps;
84 float x_adj_int = roundf (x_adj);
85 float x_adj_frac = x_adj - x_adj_int;
86 int x_adj_log2;
87 float x_adj_mant = __frexpf (x_adj, &x_adj_log2);
88 if (x_adj_mant < M_SQRT1_2f)
89 {
90 x_adj_log2--;
91 x_adj_mant *= 2.0f;
92 }
93 *exp2_adj = x_adj_log2 * (int) x_adj_int;
94 float ret = (__ieee754_powf (x_adj_mant, x_adj)
95 * __ieee754_exp2f (x_adj_log2 * x_adj_frac)
96 * __ieee754_expf (-x_adj)
97 * sqrtf (2 * M_PIf / x_adj)
98 / prod);
99 exp_adj += x_eps * __ieee754_logf (x_adj);
100 float bsum = gamma_coeff[NCOEFF - 1];
101 float x_adj2 = x_adj * x_adj;
102 for (size_t i = 1; i <= NCOEFF - 1; i++)
103 bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i];
104 exp_adj += bsum / x_adj;
105 return ret + ret * __expm1f (exp_adj);
106 }
107}
108
109float
110__ieee754_gammaf_r (float x, int *signgamp)
111{
112 int32_t hx;
113 float ret;
114
115 GET_FLOAT_WORD (hx, x);
116
117 if (__glibc_unlikely ((hx & 0x7fffffff) == 0))
118 {
119 /* Return value for x == 0 is Inf with divide by zero exception. */
120 *signgamp = 0;
121 return 1.0 / x;
122 }
123 if (__builtin_expect (hx < 0, 0)
124 && (uint32_t) hx < 0xff800000 && rintf (x) == x)
125 {
126 /* Return value for integer x < 0 is NaN with invalid exception. */
127 *signgamp = 0;
128 return (x - x) / (x - x);
129 }
130 if (__glibc_unlikely (hx == 0xff800000))
131 {
132 /* x == -Inf. According to ISO this is NaN. */
133 *signgamp = 0;
134 return x - x;
135 }
136 if (__glibc_unlikely ((hx & 0x7f800000) == 0x7f800000))
137 {
138 /* Positive infinity (return positive infinity) or NaN (return
139 NaN). */
140 *signgamp = 0;
141 return x + x;
142 }
143
144 if (x >= 36.0f)
145 {
146 /* Overflow. */
147 *signgamp = 0;
148 ret = math_narrow_eval (FLT_MAX * FLT_MAX);
149 return ret;
150 }
151 else
152 {
153 SET_RESTORE_ROUNDF (FE_TONEAREST);
154 if (x > 0.0f)
155 {
156 *signgamp = 0;
157 int exp2_adj;
158 float tret = gammaf_positive (x, &exp2_adj);
159 ret = __scalbnf (tret, exp2_adj);
160 }
161 else if (x >= -FLT_EPSILON / 4.0f)
162 {
163 *signgamp = 0;
164 ret = 1.0f / x;
165 }
166 else
167 {
168 float tx = truncf (x);
169 *signgamp = (tx == 2.0f * truncf (tx / 2.0f)) ? -1 : 1;
170 if (x <= -42.0f)
171 /* Underflow. */
172 ret = FLT_MIN * FLT_MIN;
173 else
174 {
175 float frac = tx - x;
176 if (frac > 0.5f)
177 frac = 1.0f - frac;
178 float sinpix = (frac <= 0.25f
179 ? __sinf (M_PIf * frac)
180 : __cosf (M_PIf * (0.5f - frac)));
181 int exp2_adj;
182 float tret = M_PIf / (-x * sinpix
183 * gammaf_positive (-x, &exp2_adj));
184 ret = __scalbnf (tret, -exp2_adj);
185 math_check_force_underflow_nonneg (ret);
186 }
187 }
188 ret = math_narrow_eval (ret);
189 }
190 if (isinf (ret) && x != 0)
191 {
192 if (*signgamp < 0)
193 {
194 ret = math_narrow_eval (-copysignf (FLT_MAX, ret) * FLT_MAX);
195 ret = -ret;
196 }
197 else
198 ret = math_narrow_eval (copysignf (FLT_MAX, ret) * FLT_MAX);
199 return ret;
200 }
201 else if (ret == 0)
202 {
203 if (*signgamp < 0)
204 {
205 ret = math_narrow_eval (-copysignf (FLT_MIN, ret) * FLT_MIN);
206 ret = -ret;
207 }
208 else
209 ret = math_narrow_eval (copysignf (FLT_MIN, ret) * FLT_MIN);
210 return ret;
211 }
212 else
213 return ret;
214}
215libm_alias_finite (__ieee754_gammaf_r, __gammaf_r)
216