1 /* Implementation of gamma function according to ISO C.
2 Copyright (C) 1997-2018 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
4 Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
5 Jakub Jelinek <jj@ultra.linux.cz, 1999.
7 The GNU C Library is free software; you can redistribute it and/or
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.
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
15 Lesser General Public License for more details.
17 You should have received a copy of the GNU Lesser General Public
18 License along with the GNU C Library; if not, see
19 <http://www.gnu.org/licenses/>. */
21 #include "quadmath-imp.h"
23 tgammaq (__float128 x
)
27 ret
= __quadmath_gammaq_r (x
, &sign
);
28 return sign
< 0 ? -ret
: ret
;
31 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
32 approximation to gamma function. */
34 static const __float128 gamma_coeff
[] =
36 0x1.5555555555555555555555555555p
-4Q
,
37 -0xb.60b60b60b60b60b60b60b60b60b8p
-12Q
,
38 0x3.4034034034034034034034034034p
-12Q
,
39 -0x2.7027027027027027027027027028p
-12Q
,
40 0x3.72a3c5631fe46ae1d4e700dca8f2p
-12Q
,
41 -0x7.daac36664f1f207daac36664f1f4p
-12Q
,
42 0x1.a41a41a41a41a41a41a41a41a41ap
-8Q
,
43 -0x7.90a1b2c3d4e5f708192a3b4c5d7p
-8Q
,
44 0x2.dfd2c703c0cfff430edfd2c703cp
-4Q
,
45 -0x1.6476701181f39edbdb9ce625987dp
+0Q
,
46 0xd.672219167002d3a7a9c886459cp
+0Q
,
47 -0x9.cd9292e6660d55b3f712eb9e07c8p
+4Q
,
48 0x8.911a740da740da740da740da741p
+8Q
,
49 -0x8.d0cc570e255bf59ff6eec24b49p
+12Q
,
52 #define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0]))
54 /* Return gamma (X), for positive X less than 1775, in the form R *
55 2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to
56 avoid overflow or underflow in intermediate calculations. */
59 gammal_positive (__float128 x
, int *exp2_adj
)
65 return expq (__quadmath_lgammaq_r (x
+ 1, &local_signgam
)) / x
;
70 return expq (__quadmath_lgammaq_r (x
, &local_signgam
));
74 /* Adjust into the range for using exp (lgamma). */
76 __float128 n
= ceilq (x
- 1.5Q
);
77 __float128 x_adj
= x
- n
;
79 __float128 prod
= __quadmath_gamma_productq (x_adj
, 0, n
, &eps
);
80 return (expq (__quadmath_lgammaq_r (x_adj
, &local_signgam
))
91 /* Adjust into the range for applying Stirling's
93 __float128 n
= ceilq (24 - x
);
95 x_eps
= (x
- (x_adj
- n
));
96 prod
= __quadmath_gamma_productq (x_adj
- n
, x_eps
, n
, &eps
);
98 /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)).
99 Compute gamma (X_ADJ + X_EPS) using Stirling's approximation,
100 starting by computing pow (X_ADJ, X_ADJ) with a power of 2
102 __float128 exp_adj
= -eps
;
103 __float128 x_adj_int
= roundq (x_adj
);
104 __float128 x_adj_frac
= x_adj
- x_adj_int
;
106 __float128 x_adj_mant
= frexpq (x_adj
, &x_adj_log2
);
107 if (x_adj_mant
< M_SQRT1_2q
)
112 *exp2_adj
= x_adj_log2
* (int) x_adj_int
;
113 __float128 ret
= (powq (x_adj_mant
, x_adj
)
114 * exp2q (x_adj_log2
* x_adj_frac
)
116 * sqrtq (2 * M_PIq
/ x_adj
)
118 exp_adj
+= x_eps
* logq (x_adj
);
119 __float128 bsum
= gamma_coeff
[NCOEFF
- 1];
120 __float128 x_adj2
= x_adj
* x_adj
;
121 for (size_t i
= 1; i
<= NCOEFF
- 1; i
++)
122 bsum
= bsum
/ x_adj2
+ gamma_coeff
[NCOEFF
- 1 - i
];
123 exp_adj
+= bsum
/ x_adj
;
124 return ret
+ ret
* expm1q (exp_adj
);
129 __quadmath_gammaq_r (__float128 x
, int *signgamp
)
135 GET_FLT128_WORDS64 (hx
, lx
, x
);
137 if (((hx
& 0x7fffffffffffffffLL
) | lx
) == 0)
139 /* Return value for x == 0 is Inf with divide by zero exception. */
143 if (hx
< 0 && (uint64_t) hx
< 0xffff000000000000ULL
&& rintq (x
) == x
)
145 /* Return value for integer x < 0 is NaN with invalid exception. */
147 return (x
- x
) / (x
- x
);
149 if (hx
== 0xffff000000000000ULL
&& lx
== 0)
151 /* x == -Inf. According to ISO this is NaN. */
155 if ((hx
& 0x7fff000000000000ULL
) == 0x7fff000000000000ULL
)
157 /* Positive infinity (return positive infinity) or NaN (return
167 return FLT128_MAX
* FLT128_MAX
;
171 SET_RESTORE_ROUNDF128 (FE_TONEAREST
);
176 ret
= gammal_positive (x
, &exp2_adj
);
177 ret
= scalbnq (ret
, exp2_adj
);
179 else if (x
>= -FLT128_EPSILON
/ 4)
186 __float128 tx
= truncq (x
);
187 *signgamp
= (tx
== 2 * truncq (tx
/ 2)) ? -1 : 1;
190 ret
= FLT128_MIN
* FLT128_MIN
;
193 __float128 frac
= tx
- x
;
196 __float128 sinpix
= (frac
<= 0.25Q
197 ? sinq (M_PIq
* frac
)
198 : cosq (M_PIq
* (0.5Q
- frac
)));
200 ret
= M_PIq
/ (-x
* sinpix
201 * gammal_positive (-x
, &exp2_adj
));
202 ret
= scalbnq (ret
, -exp2_adj
);
203 math_check_force_underflow_nonneg (ret
);
207 if (isinfq (ret
) && x
!= 0)
210 return -(-copysignq (FLT128_MAX
, ret
) * FLT128_MAX
);
212 return copysignq (FLT128_MAX
, ret
) * FLT128_MAX
;
217 return -(-copysignq (FLT128_MIN
, ret
) * FLT128_MIN
);
219 return copysignq (FLT128_MIN
, ret
) * FLT128_MIN
;