sysdeps/ieee754/ldbl-96/e_gammal_r.c

   1 /* Implementation of gamma function according to ISO C.
   2    Copyright (C) 1997-2015 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_private.h>
  22 #include <float.h>
  23
  24 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
  25    approximation to gamma function.  */
  26
  27 static const long double gamma_coeff[] =
  28   {
  29     0x1.5555555555555556p-4L,
  30     -0xb.60b60b60b60b60bp-12L,
  31     0x3.4034034034034034p-12L,
  32     -0x2.7027027027027028p-12L,
  33     0x3.72a3c5631fe46aep-12L,
  34     -0x7.daac36664f1f208p-12L,
  35     0x1.a41a41a41a41a41ap-8L,
  36     -0x7.90a1b2c3d4e5f708p-8L,
  37   };
  38
  39 #define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0]))
  40
  41 /* Return gamma (X), for positive X less than 1766, in the form R *
  42    2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to
  43    avoid overflow or underflow in intermediate calculations.  */
  44
  45 static long double
  46 gammal_positive (long double x, int *exp2_adj)
  47 {
  48   int local_signgam;
  49   if (x < 0.5L)
  50     {
  51       *exp2_adj = 0;
  52       return __ieee754_expl (__ieee754_lgammal_r (x + 1, &local_signgam)) / x;
  53     }
  54   else if (x <= 1.5L)
  55     {
  56       *exp2_adj = 0;
  57       return __ieee754_expl (__ieee754_lgammal_r (x, &local_signgam));
  58     }
  59   else if (x < 7.5L)
  60     {
  61       /* Adjust into the range for using exp (lgamma).  */
  62       *exp2_adj = 0;
  63       long double n = __ceill (x - 1.5L);
  64       long double x_adj = x - n;
  65       long double eps;
  66       long double prod = __gamma_productl (x_adj, 0, n, &eps);
  67       return (__ieee754_expl (__ieee754_lgammal_r (x_adj, &local_signgam))
  68               * prod * (1.0L + eps));
  69     }
  70   else
  71     {
  72       long double eps = 0;
  73       long double x_eps = 0;
  74       long double x_adj = x;
  75       long double prod = 1;
  76       if (x < 13.0L)
  77         {
  78           /* Adjust into the range for applying Stirling's
  79              approximation.  */
  80           long double n = __ceill (13.0L - x);
  81           x_adj = x + n;
  82           x_eps = (x - (x_adj - n));
  83           prod = __gamma_productl (x_adj - n, x_eps, n, &eps);
  84         }
  85       /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)).
  86          Compute gamma (X_ADJ + X_EPS) using Stirling's approximation,
  87          starting by computing pow (X_ADJ, X_ADJ) with a power of 2
  88          factored out.  */
  89       long double exp_adj = -eps;
  90       long double x_adj_int = __roundl (x_adj);
  91       long double x_adj_frac = x_adj - x_adj_int;
  92       int x_adj_log2;
  93       long double x_adj_mant = __frexpl (x_adj, &x_adj_log2);
  94       if (x_adj_mant < M_SQRT1_2l)
  95         {
  96           x_adj_log2--;
  97           x_adj_mant *= 2.0L;
  98         }
  99       *exp2_adj = x_adj_log2 * (int) x_adj_int;
 100       long double ret = (__ieee754_powl (x_adj_mant, x_adj)
 101                          * __ieee754_exp2l (x_adj_log2 * x_adj_frac)
 102                          * __ieee754_expl (-x_adj)
 103                          * __ieee754_sqrtl (2 * M_PIl / x_adj)
 104                          / prod);
 105       exp_adj += x_eps * __ieee754_logl (x);
 106       long double bsum = gamma_coeff[NCOEFF - 1];
 107       long double x_adj2 = x_adj * x_adj;
 108       for (size_t i = 1; i <= NCOEFF - 1; i++)
 109         bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i];
 110       exp_adj += bsum / x_adj;
 111       return ret + ret * __expm1l (exp_adj);
 112     }
 113 }
 114
 115 long double
 116 __ieee754_gammal_r (long double x, int *signgamp)
 117 {
 118   u_int32_t es, hx, lx;
 119
 120   GET_LDOUBLE_WORDS (es, hx, lx, x);
 121
 122   if (__glibc_unlikely (((es & 0x7fff) | hx | lx) == 0))
 123     {
 124       /* Return value for x == 0 is Inf with divide by zero exception.  */
 125       *signgamp = 0;
 126       return 1.0 / x;
 127     }
 128   if (__glibc_unlikely (es == 0xffffffff && ((hx & 0x7fffffff) | lx) == 0))
 129     {
 130       /* x == -Inf.  According to ISO this is NaN.  */
 131       *signgamp = 0;
 132       return x - x;
 133     }
 134   if (__glibc_unlikely ((es & 0x7fff) == 0x7fff))
 135     {
 136       /* Positive infinity (return positive infinity) or NaN (return
 137          NaN).  */
 138       *signgamp = 0;
 139       return x + x;
 140     }
 141   if (__builtin_expect ((es & 0x8000) != 0, 0) && __rintl (x) == x)
 142     {
 143       /* Return value for integer x < 0 is NaN with invalid exception.  */
 144       *signgamp = 0;
 145       return (x - x) / (x - x);
 146     }
 147
 148   if (x >= 1756.0L)
 149     {
 150       /* Overflow.  */
 151       *signgamp = 0;
 152       return LDBL_MAX * LDBL_MAX;
 153     }
 154   else if (x > 0.0L)
 155     {
 156       *signgamp = 0;
 157       int exp2_adj;
 158       long double ret = gammal_positive (x, &exp2_adj);
 159       return __scalbnl (ret, exp2_adj);
 160     }
 161   else if (x >= -LDBL_EPSILON / 4.0L)
 162     {
 163       *signgamp = 0;
 164       return 1.0f / x;
 165     }
 166   else
 167     {
 168       long double tx = __truncl (x);
 169       *signgamp = (tx == 2.0L * __truncl (tx / 2.0L)) ? -1 : 1;
 170       if (x <= -1766.0L)
 171         /* Underflow.  */
 172         return LDBL_MIN * LDBL_MIN;
 173       long double frac = tx - x;
 174       if (frac > 0.5L)
 175         frac = 1.0L - frac;
 176       long double sinpix = (frac <= 0.25L
 177                             ? __sinl (M_PIl * frac)
 178                             : __cosl (M_PIl * (0.5L - frac)));
 179       int exp2_adj;
 180       long double ret = M_PIl / (-x * sinpix
 181                                  * gammal_positive (-x, &exp2_adj));
 182       return __scalbnl (ret, -exp2_adj);
 183     }
 184 }
 185 strong_alias (__ieee754_gammal_r, __gammal_r_finite)