libquadmath/math/tgammaq.c

   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 }