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. 6 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. 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 15 Lesser General Public License for more details. 16 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/>. */ 20 21 #include "quadmath-imp.h" 22 __float128 23 tgammaq (__float128 x) 24 { 25 int sign; 26 __float128 ret; 27 ret = __quadmath_gammaq_r (x, &sign); 28 return sign < 0 ? -ret : ret; 29 } 30 31 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's 32 approximation to gamma function. */ 33 34 static const __float128 gamma_coeff[] = 35 { 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, 50 }; 51 52 #define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0])) 53 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. */ 57 58 static __float128 59 gammal_positive (__float128 x, int *exp2_adj) 60 { 61 int local_signgam; 62 if (x < 0.5Q) 63 { 64 *exp2_adj = 0; 65 return expq (__quadmath_lgammaq_r (x + 1, &local_signgam)) / x; 66 } 67 else if (x <= 1.5Q) 68 { 69 *exp2_adj = 0; 70 return expq (__quadmath_lgammaq_r (x, &local_signgam)); 71 } 72 else if (x < 12.5Q) 73 { 74 /* Adjust into the range for using exp (lgamma). */ 75 *exp2_adj = 0; 76 __float128 n = ceilq (x - 1.5Q); 77 __float128 x_adj = x - n; 78 __float128 eps; 79 __float128 prod = __quadmath_gamma_productq (x_adj, 0, n, &eps); 80 return (expq (__quadmath_lgammaq_r (x_adj, &local_signgam)) 81 * prod * (1 + eps)); 82 } 83 else 84 { 85 __float128 eps = 0; 86 __float128 x_eps = 0; 87 __float128 x_adj = x; 88 __float128 prod = 1; 89 if (x < 24) 90 { 91 /* Adjust into the range for applying Stirling's 92 approximation. */ 93 __float128 n = ceilq (24 - x); 94 x_adj = x + n; 95 x_eps = (x - (x_adj - n)); 96 prod = __quadmath_gamma_productq (x_adj - n, x_eps, n, &eps); 97 } 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 101 factored out. */ 102 __float128 exp_adj = -eps; 103 __float128 x_adj_int = roundq (x_adj); 104 __float128 x_adj_frac = x_adj - x_adj_int; 105 int x_adj_log2; 106 __float128 x_adj_mant = frexpq (x_adj, &x_adj_log2); 107 if (x_adj_mant < M_SQRT1_2q) 108 { 109 x_adj_log2--; 110 x_adj_mant *= 2; 111 } 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) 115 * expq (-x_adj) 116 * sqrtq (2 * M_PIq / x_adj) 117 / prod); 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); 125 } 126 } 127 128 __float128 129 __quadmath_gammaq_r (__float128 x, int *signgamp) 130 { 131 int64_t hx; 132 uint64_t lx; 133 __float128 ret; 134 135 GET_FLT128_WORDS64 (hx, lx, x); 136 137 if (((hx & 0x7fffffffffffffffLL) | lx) == 0) 138 { 139 /* Return value for x == 0 is Inf with divide by zero exception. */ 140 *signgamp = 0; 141 return 1.0 / x; 142 } 143 if (hx < 0 && (uint64_t) hx < 0xffff000000000000ULL && rintq (x) == x) 144 { 145 /* Return value for integer x < 0 is NaN with invalid exception. */ 146 *signgamp = 0; 147 return (x - x) / (x - x); 148 } 149 if (hx == 0xffff000000000000ULL && lx == 0) 150 { 151 /* x == -Inf. According to ISO this is NaN. */ 152 *signgamp = 0; 153 return x - x; 154 } 155 if ((hx & 0x7fff000000000000ULL) == 0x7fff000000000000ULL) 156 { 157 /* Positive infinity (return positive infinity) or NaN (return 158 NaN). */ 159 *signgamp = 0; 160 return x + x; 161 } 162 163 if (x >= 1756) 164 { 165 /* Overflow. */ 166 *signgamp = 0; 167 return FLT128_MAX * FLT128_MAX; 168 } 169 else 170 { 171 SET_RESTORE_ROUNDF128 (FE_TONEAREST); 172 if (x > 0) 173 { 174 *signgamp = 0; 175 int exp2_adj; 176 ret = gammal_positive (x, &exp2_adj); 177 ret = scalbnq (ret, exp2_adj); 178 } 179 else if (x >= -FLT128_EPSILON / 4) 180 { 181 *signgamp = 0; 182 ret = 1 / x; 183 } 184 else 185 { 186 __float128 tx = truncq (x); 187 *signgamp = (tx == 2 * truncq (tx / 2)) ? -1 : 1; 188 if (x <= -1775) 189 /* Underflow. */ 190 ret = FLT128_MIN * FLT128_MIN; 191 else 192 { 193 __float128 frac = tx - x; 194 if (frac > 0.5Q) 195 frac = 1 - frac; 196 __float128 sinpix = (frac <= 0.25Q 197 ? sinq (M_PIq * frac) 198 : cosq (M_PIq * (0.5Q - frac))); 199 int exp2_adj; 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); 204 } 205 } 206 } 207 if (isinfq (ret) && x != 0) 208 { 209 if (*signgamp < 0) 210 return -(-copysignq (FLT128_MAX, ret) * FLT128_MAX); 211 else 212 return copysignq (FLT128_MAX, ret) * FLT128_MAX; 213 } 214 else if (ret == 0) 215 { 216 if (*signgamp < 0) 217 return -(-copysignq (FLT128_MIN, ret) * FLT128_MIN); 218 else 219 return copysignq (FLT128_MIN, ret) * FLT128_MIN; 220 } 221 else 222 return ret; 223 } 224