]>
Commit | Line | Data |
---|---|---|
d705269e | 1 | /* Implementation of gamma function according to ISO C. |
04277e02 | 2 | Copyright (C) 1997-2019 Free Software Foundation, Inc. |
c131718c | 3 | This file is part of the GNU C Library. |
abfbdde1 | 4 | Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and |
0ac5ae23 | 5 | Jakub Jelinek <jj@ultra.linux.cz, 1999. |
c131718c UD |
6 | |
7 | The GNU C Library is free software; you can redistribute it and/or | |
41bdb6e2 AJ |
8 | modify it under the terms of the GNU Lesser General Public |
9 | License as published by the Free Software Foundation; either | |
10 | version 2.1 of the License, or (at your option) any later version. | |
c131718c UD |
11 | |
12 | The GNU C Library is distributed in the hope that it will be useful, | |
13 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
41bdb6e2 | 15 | Lesser General Public License for more details. |
c131718c | 16 | |
41bdb6e2 | 17 | You should have received a copy of the GNU Lesser General Public |
59ba27a6 | 18 | License along with the GNU C Library; if not, see |
5a82c748 | 19 | <https://www.gnu.org/licenses/>. */ |
c131718c | 20 | |
d705269e UD |
21 | #include <math.h> |
22 | #include <math_private.h> | |
70e2ba33 | 23 | #include <fenv_private.h> |
8f5b00d3 | 24 | #include <math-underflow.h> |
d8cd06db | 25 | #include <float.h> |
d705269e | 26 | |
d8cd06db JM |
27 | /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's |
28 | approximation to gamma function. */ | |
29 | ||
15089e04 | 30 | static const _Float128 gamma_coeff[] = |
d8cd06db | 31 | { |
02bbfb41 PM |
32 | L(0x1.5555555555555555555555555555p-4), |
33 | L(-0xb.60b60b60b60b60b60b60b60b60b8p-12), | |
34 | L(0x3.4034034034034034034034034034p-12), | |
35 | L(-0x2.7027027027027027027027027028p-12), | |
36 | L(0x3.72a3c5631fe46ae1d4e700dca8f2p-12), | |
37 | L(-0x7.daac36664f1f207daac36664f1f4p-12), | |
38 | L(0x1.a41a41a41a41a41a41a41a41a41ap-8), | |
39 | L(-0x7.90a1b2c3d4e5f708192a3b4c5d7p-8), | |
40 | L(0x2.dfd2c703c0cfff430edfd2c703cp-4), | |
41 | L(-0x1.6476701181f39edbdb9ce625987dp+0), | |
42 | L(0xd.672219167002d3a7a9c886459cp+0), | |
43 | L(-0x9.cd9292e6660d55b3f712eb9e07c8p+4), | |
44 | L(0x8.911a740da740da740da740da741p+8), | |
45 | L(-0x8.d0cc570e255bf59ff6eec24b49p+12), | |
d8cd06db JM |
46 | }; |
47 | ||
48 | #define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0])) | |
49 | ||
50 | /* Return gamma (X), for positive X less than 1775, in the form R * | |
51 | 2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to | |
52 | avoid overflow or underflow in intermediate calculations. */ | |
53 | ||
15089e04 PM |
54 | static _Float128 |
55 | gammal_positive (_Float128 x, int *exp2_adj) | |
d8cd06db JM |
56 | { |
57 | int local_signgam; | |
02bbfb41 | 58 | if (x < L(0.5)) |
d8cd06db JM |
59 | { |
60 | *exp2_adj = 0; | |
61 | return __ieee754_expl (__ieee754_lgammal_r (x + 1, &local_signgam)) / x; | |
62 | } | |
02bbfb41 | 63 | else if (x <= L(1.5)) |
d8cd06db JM |
64 | { |
65 | *exp2_adj = 0; | |
66 | return __ieee754_expl (__ieee754_lgammal_r (x, &local_signgam)); | |
67 | } | |
02bbfb41 | 68 | else if (x < L(12.5)) |
d8cd06db JM |
69 | { |
70 | /* Adjust into the range for using exp (lgamma). */ | |
71 | *exp2_adj = 0; | |
71223ef9 | 72 | _Float128 n = ceill (x - L(1.5)); |
15089e04 PM |
73 | _Float128 x_adj = x - n; |
74 | _Float128 eps; | |
75 | _Float128 prod = __gamma_productl (x_adj, 0, n, &eps); | |
d8cd06db | 76 | return (__ieee754_expl (__ieee754_lgammal_r (x_adj, &local_signgam)) |
02bbfb41 | 77 | * prod * (1 + eps)); |
d8cd06db JM |
78 | } |
79 | else | |
80 | { | |
15089e04 PM |
81 | _Float128 eps = 0; |
82 | _Float128 x_eps = 0; | |
83 | _Float128 x_adj = x; | |
84 | _Float128 prod = 1; | |
02bbfb41 | 85 | if (x < 24) |
d8cd06db JM |
86 | { |
87 | /* Adjust into the range for applying Stirling's | |
88 | approximation. */ | |
71223ef9 | 89 | _Float128 n = ceill (24 - x); |
d8cd06db JM |
90 | x_adj = x + n; |
91 | x_eps = (x - (x_adj - n)); | |
92 | prod = __gamma_productl (x_adj - n, x_eps, n, &eps); | |
93 | } | |
94 | /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)). | |
95 | Compute gamma (X_ADJ + X_EPS) using Stirling's approximation, | |
96 | starting by computing pow (X_ADJ, X_ADJ) with a power of 2 | |
97 | factored out. */ | |
15089e04 | 98 | _Float128 exp_adj = -eps; |
9755bc46 | 99 | _Float128 x_adj_int = roundl (x_adj); |
15089e04 | 100 | _Float128 x_adj_frac = x_adj - x_adj_int; |
d8cd06db | 101 | int x_adj_log2; |
15089e04 | 102 | _Float128 x_adj_mant = __frexpl (x_adj, &x_adj_log2); |
d8cd06db JM |
103 | if (x_adj_mant < M_SQRT1_2l) |
104 | { | |
105 | x_adj_log2--; | |
02bbfb41 | 106 | x_adj_mant *= 2; |
d8cd06db JM |
107 | } |
108 | *exp2_adj = x_adj_log2 * (int) x_adj_int; | |
15089e04 | 109 | _Float128 ret = (__ieee754_powl (x_adj_mant, x_adj) |
de6b6d14 PM |
110 | * __ieee754_exp2l (x_adj_log2 * x_adj_frac) |
111 | * __ieee754_expl (-x_adj) | |
f67a8147 | 112 | * sqrtl (2 * M_PIl / x_adj) |
de6b6d14 | 113 | / prod); |
e02920bc | 114 | exp_adj += x_eps * __ieee754_logl (x_adj); |
15089e04 PM |
115 | _Float128 bsum = gamma_coeff[NCOEFF - 1]; |
116 | _Float128 x_adj2 = x_adj * x_adj; | |
d8cd06db JM |
117 | for (size_t i = 1; i <= NCOEFF - 1; i++) |
118 | bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i]; | |
119 | exp_adj += bsum / x_adj; | |
120 | return ret + ret * __expm1l (exp_adj); | |
121 | } | |
122 | } | |
d705269e | 123 | |
15089e04 PM |
124 | _Float128 |
125 | __ieee754_gammal_r (_Float128 x, int *signgamp) | |
d705269e | 126 | { |
abfbdde1 | 127 | int64_t hx; |
24ab7723 | 128 | uint64_t lx; |
15089e04 | 129 | _Float128 ret; |
d705269e | 130 | |
abfbdde1 | 131 | GET_LDOUBLE_WORDS64 (hx, lx, x); |
d705269e | 132 | |
abfbdde1 | 133 | if (((hx & 0x7fffffffffffffffLL) | lx) == 0) |
b3fc5f84 | 134 | { |
52495f29 | 135 | /* Return value for x == 0 is Inf with divide by zero exception. */ |
b3fc5f84 | 136 | *signgamp = 0; |
52495f29 | 137 | return 1.0 / x; |
b3fc5f84 | 138 | } |
f29b6f17 | 139 | if (hx < 0 && (uint64_t) hx < 0xffff000000000000ULL && rintl (x) == x) |
d705269e UD |
140 | { |
141 | /* Return value for integer x < 0 is NaN with invalid exception. */ | |
b3fc5f84 | 142 | *signgamp = 0; |
d705269e UD |
143 | return (x - x) / (x - x); |
144 | } | |
52e1b618 UD |
145 | if (hx == 0xffff000000000000ULL && lx == 0) |
146 | { | |
147 | /* x == -Inf. According to ISO this is NaN. */ | |
148 | *signgamp = 0; | |
149 | return x - x; | |
150 | } | |
d8cd06db JM |
151 | if ((hx & 0x7fff000000000000ULL) == 0x7fff000000000000ULL) |
152 | { | |
153 | /* Positive infinity (return positive infinity) or NaN (return | |
154 | NaN). */ | |
155 | *signgamp = 0; | |
156 | return x + x; | |
157 | } | |
d705269e | 158 | |
02bbfb41 | 159 | if (x >= 1756) |
d8cd06db JM |
160 | { |
161 | /* Overflow. */ | |
162 | *signgamp = 0; | |
163 | return LDBL_MAX * LDBL_MAX; | |
164 | } | |
e02920bc | 165 | else |
d8cd06db | 166 | { |
e02920bc | 167 | SET_RESTORE_ROUNDL (FE_TONEAREST); |
02bbfb41 | 168 | if (x > 0) |
e02920bc JM |
169 | { |
170 | *signgamp = 0; | |
171 | int exp2_adj; | |
172 | ret = gammal_positive (x, &exp2_adj); | |
173 | ret = __scalbnl (ret, exp2_adj); | |
174 | } | |
02bbfb41 | 175 | else if (x >= -LDBL_EPSILON / 4) |
e02920bc JM |
176 | { |
177 | *signgamp = 0; | |
02bbfb41 | 178 | ret = 1 / x; |
e02920bc JM |
179 | } |
180 | else | |
181 | { | |
7abf97be JM |
182 | _Float128 tx = truncl (x); |
183 | *signgamp = (tx == 2 * truncl (tx / 2)) ? -1 : 1; | |
02bbfb41 | 184 | if (x <= -1775) |
e02920bc JM |
185 | /* Underflow. */ |
186 | ret = LDBL_MIN * LDBL_MIN; | |
187 | else | |
188 | { | |
15089e04 | 189 | _Float128 frac = tx - x; |
02bbfb41 PM |
190 | if (frac > L(0.5)) |
191 | frac = 1 - frac; | |
192 | _Float128 sinpix = (frac <= L(0.25) | |
de6b6d14 | 193 | ? __sinl (M_PIl * frac) |
02bbfb41 | 194 | : __cosl (M_PIl * (L(0.5) - frac))); |
e02920bc JM |
195 | int exp2_adj; |
196 | ret = M_PIl / (-x * sinpix | |
197 | * gammal_positive (-x, &exp2_adj)); | |
198 | ret = __scalbnl (ret, -exp2_adj); | |
d96164c3 | 199 | math_check_force_underflow_nonneg (ret); |
e02920bc JM |
200 | } |
201 | } | |
d8cd06db | 202 | } |
e02920bc | 203 | if (isinf (ret) && x != 0) |
d8cd06db | 204 | { |
e02920bc | 205 | if (*signgamp < 0) |
81dca813 | 206 | return -(-copysignl (LDBL_MAX, ret) * LDBL_MAX); |
e02920bc | 207 | else |
81dca813 | 208 | return copysignl (LDBL_MAX, ret) * LDBL_MAX; |
d8cd06db | 209 | } |
e02920bc | 210 | else if (ret == 0) |
d8cd06db | 211 | { |
e02920bc | 212 | if (*signgamp < 0) |
81dca813 | 213 | return -(-copysignl (LDBL_MIN, ret) * LDBL_MIN); |
e02920bc | 214 | else |
81dca813 | 215 | return copysignl (LDBL_MIN, ret) * LDBL_MIN; |
d8cd06db | 216 | } |
e02920bc JM |
217 | else |
218 | return ret; | |
d705269e | 219 | } |
0ac5ae23 | 220 | strong_alias (__ieee754_gammal_r, __gammal_r_finite) |