1/* Implementation of gamma function according to ISO C.
2 Copyright (C) 1997-2019 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
5
6 The GNU C Library is free software; you can redistribute it and/or
7 modify it under the terms of the GNU Lesser General Public
8 License as published by the Free Software Foundation; either
9 version 2.1 of the License, or (at your option) any later version.
10
11 The GNU C Library is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 Lesser General Public License for more details.
15
16 You should have received a copy of the GNU Lesser General Public
17 License along with the GNU C Library; if not, see
18 <http://www.gnu.org/licenses/>. */
19
20#include <math.h>
21#include <math-narrow-eval.h>
22#include <math_private.h>
23#include <fenv_private.h>
24#include <math-underflow.h>
25#include <float.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 double gamma_coeff[] =
31 {
32 0x1.5555555555555p-4,
33 -0xb.60b60b60b60b8p-12,
34 0x3.4034034034034p-12,
35 -0x2.7027027027028p-12,
36 0x3.72a3c5631fe46p-12,
37 -0x7.daac36664f1f4p-12,
38 };
39
40#define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0]))
41
42/* Return gamma (X), for positive X less than 184, in the form R *
43 2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to
44 avoid overflow or underflow in intermediate calculations. */
45
46static double
47gamma_positive (double x, int *exp2_adj)
48{
49 int local_signgam;
50 if (x < 0.5)
51 {
52 *exp2_adj = 0;
53 return __ieee754_exp (__ieee754_lgamma_r (x + 1, &local_signgam)) / x;
54 }
55 else if (x <= 1.5)
56 {
57 *exp2_adj = 0;
58 return __ieee754_exp (__ieee754_lgamma_r (x, &local_signgam));
59 }
60 else if (x < 6.5)
61 {
62 /* Adjust into the range for using exp (lgamma). */
63 *exp2_adj = 0;
64 double n = ceil (x - 1.5);
65 double x_adj = x - n;
66 double eps;
67 double prod = __gamma_product (x_adj, 0, n, &eps);
68 return (__ieee754_exp (__ieee754_lgamma_r (x_adj, &local_signgam))
69 * prod * (1.0 + eps));
70 }
71 else
72 {
73 double eps = 0;
74 double x_eps = 0;
75 double x_adj = x;
76 double prod = 1;
77 if (x < 12.0)
78 {
79 /* Adjust into the range for applying Stirling's
80 approximation. */
81 double n = ceil (12.0 - x);
82 x_adj = math_narrow_eval (x + n);
83 x_eps = (x - (x_adj - n));
84 prod = __gamma_product (x_adj - n, x_eps, n, &eps);
85 }
86 /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)).
87 Compute gamma (X_ADJ + X_EPS) using Stirling's approximation,
88 starting by computing pow (X_ADJ, X_ADJ) with a power of 2
89 factored out. */
90 double exp_adj = -eps;
91 double x_adj_int = round (x_adj);
92 double x_adj_frac = x_adj - x_adj_int;
93 int x_adj_log2;
94 double x_adj_mant = __frexp (x_adj, &x_adj_log2);
95 if (x_adj_mant < M_SQRT1_2)
96 {
97 x_adj_log2--;
98 x_adj_mant *= 2.0;
99 }
100 *exp2_adj = x_adj_log2 * (int) x_adj_int;
101 double ret = (__ieee754_pow (x_adj_mant, x_adj)
102 * __ieee754_exp2 (x_adj_log2 * x_adj_frac)
103 * __ieee754_exp (-x_adj)
104 * sqrt (2 * M_PI / x_adj)
105 / prod);
106 exp_adj += x_eps * __ieee754_log (x_adj);
107 double bsum = gamma_coeff[NCOEFF - 1];
108 double x_adj2 = x_adj * x_adj;
109 for (size_t i = 1; i <= NCOEFF - 1; i++)
110 bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i];
111 exp_adj += bsum / x_adj;
112 return ret + ret * __expm1 (exp_adj);
113 }
114}
115
116double
117__ieee754_gamma_r (double x, int *signgamp)
118{
119 int32_t hx;
120 uint32_t lx;
121 double ret;
122
123 EXTRACT_WORDS (hx, lx, x);
124
125 if (__glibc_unlikely (((hx & 0x7fffffff) | lx) == 0))
126 {
127 /* Return value for x == 0 is Inf with divide by zero exception. */
128 *signgamp = 0;
129 return 1.0 / x;
130 }
131 if (__builtin_expect (hx < 0, 0)
132 && (uint32_t) hx < 0xfff00000 && rint (x) == x)
133 {
134 /* Return value for integer x < 0 is NaN with invalid exception. */
135 *signgamp = 0;
136 return (x - x) / (x - x);
137 }
138 if (__glibc_unlikely ((unsigned int) hx == 0xfff00000 && lx == 0))
139 {
140 /* x == -Inf. According to ISO this is NaN. */
141 *signgamp = 0;
142 return x - x;
143 }
144 if (__glibc_unlikely ((hx & 0x7ff00000) == 0x7ff00000))
145 {
146 /* Positive infinity (return positive infinity) or NaN (return
147 NaN). */
148 *signgamp = 0;
149 return x + x;
150 }
151
152 if (x >= 172.0)
153 {
154 /* Overflow. */
155 *signgamp = 0;
156 ret = math_narrow_eval (DBL_MAX * DBL_MAX);
157 return ret;
158 }
159 else
160 {
161 SET_RESTORE_ROUND (FE_TONEAREST);
162 if (x > 0.0)
163 {
164 *signgamp = 0;
165 int exp2_adj;
166 double tret = gamma_positive (x, &exp2_adj);
167 ret = __scalbn (tret, exp2_adj);
168 }
169 else if (x >= -DBL_EPSILON / 4.0)
170 {
171 *signgamp = 0;
172 ret = 1.0 / x;
173 }
174 else
175 {
176 double tx = trunc (x);
177 *signgamp = (tx == 2.0 * trunc (tx / 2.0)) ? -1 : 1;
178 if (x <= -184.0)
179 /* Underflow. */
180 ret = DBL_MIN * DBL_MIN;
181 else
182 {
183 double frac = tx - x;
184 if (frac > 0.5)
185 frac = 1.0 - frac;
186 double sinpix = (frac <= 0.25
187 ? __sin (M_PI * frac)
188 : __cos (M_PI * (0.5 - frac)));
189 int exp2_adj;
190 double tret = M_PI / (-x * sinpix
191 * gamma_positive (-x, &exp2_adj));
192 ret = __scalbn (tret, -exp2_adj);
193 math_check_force_underflow_nonneg (ret);
194 }
195 }
196 ret = math_narrow_eval (ret);
197 }
198 if (isinf (ret) && x != 0)
199 {
200 if (*signgamp < 0)
201 {
202 ret = math_narrow_eval (-copysign (DBL_MAX, ret) * DBL_MAX);
203 ret = -ret;
204 }
205 else
206 ret = math_narrow_eval (copysign (DBL_MAX, ret) * DBL_MAX);
207 return ret;
208 }
209 else if (ret == 0)
210 {
211 if (*signgamp < 0)
212 {
213 ret = math_narrow_eval (-copysign (DBL_MIN, ret) * DBL_MIN);
214 ret = -ret;
215 }
216 else
217 ret = math_narrow_eval (copysign (DBL_MIN, ret) * DBL_MIN);
218 return ret;
219 }
220 else
221 return ret;
222}
223strong_alias (__ieee754_gamma_r, __gamma_r_finite)
224