sysdeps/ieee754/dbl-64/e_gamma_r.c

   1 /* Implementation of gamma function according to ISO C.
   2    Copyright (C) 1997-2014 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 double gamma_coeff[] =
  28   {
  29     0x1.5555555555555p-4,
  30     -0xb.60b60b60b60b8p-12,
  31     0x3.4034034034034p-12,
  32     -0x2.7027027027028p-12,
  33     0x3.72a3c5631fe46p-12,
  34     -0x7.daac36664f1f4p-12,
  35   };
  36
  37 #define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0]))
  38
  39 /* Return gamma (X), for positive X less than 184, in the form R *
  40    2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to
  41    avoid overflow or underflow in intermediate calculations.  */
  42
  43 static double
  44 gamma_positive (double x, int *exp2_adj)
  45 {
  46   int local_signgam;
  47   if (x < 0.5)
  48     {
  49       *exp2_adj = 0;
  50       return __ieee754_exp (__ieee754_lgamma_r (x + 1, &local_signgam)) / x;
  51     }
  52   else if (x <= 1.5)
  53     {
  54       *exp2_adj = 0;
  55       return __ieee754_exp (__ieee754_lgamma_r (x, &local_signgam));
  56     }
  57   else if (x < 6.5)
  58     {
  59       /* Adjust into the range for using exp (lgamma).  */
  60       *exp2_adj = 0;
  61       double n = __ceil (x - 1.5);
  62       double x_adj = x - n;
  63       double eps;
  64       double prod = __gamma_product (x_adj, 0, n, &eps);
  65       return (__ieee754_exp (__ieee754_lgamma_r (x_adj, &local_signgam))
  66               * prod * (1.0 + eps));
  67     }
  68   else
  69     {
  70       double eps = 0;
  71       double x_eps = 0;
  72       double x_adj = x;
  73       double prod = 1;
  74       if (x < 12.0)
  75         {
  76           /* Adjust into the range for applying Stirling's
  77              approximation.  */
  78           double n = __ceil (12.0 - x);
  79 #if FLT_EVAL_METHOD != 0
  80           volatile
  81 #endif
  82           double x_tmp = x + n;
  83           x_adj = x_tmp;
  84           x_eps = (x - (x_adj - n));
  85           prod = __gamma_product (x_adj - n, x_eps, n, &eps);
  86         }
  87       /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)).
  88          Compute gamma (X_ADJ + X_EPS) using Stirling's approximation,
  89          starting by computing pow (X_ADJ, X_ADJ) with a power of 2
  90          factored out.  */
  91       double exp_adj = -eps;
  92       double x_adj_int = __round (x_adj);
  93       double x_adj_frac = x_adj - x_adj_int;
  94       int x_adj_log2;
  95       double x_adj_mant = __frexp (x_adj, &x_adj_log2);
  96       if (x_adj_mant < M_SQRT1_2)
  97         {
  98           x_adj_log2--;
  99           x_adj_mant *= 2.0;
 100         }
 101       *exp2_adj = x_adj_log2 * (int) x_adj_int;
 102       double ret = (__ieee754_pow (x_adj_mant, x_adj)
 103                     * __ieee754_exp2 (x_adj_log2 * x_adj_frac)
 104                     * __ieee754_exp (-x_adj)
 105                     * __ieee754_sqrt (2 * M_PI / x_adj)
 106                     / prod);
 107       exp_adj += x_eps * __ieee754_log (x);
 108       double bsum = gamma_coeff[NCOEFF - 1];
 109       double x_adj2 = x_adj * x_adj;
 110       for (size_t i = 1; i <= NCOEFF - 1; i++)
 111         bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i];
 112       exp_adj += bsum / x_adj;
 113       return ret + ret * __expm1 (exp_adj);
 114     }
 115 }
 116
 117 double
 118 __ieee754_gamma_r (double x, int *signgamp)
 119 {
 120   int32_t hx;
 121   u_int32_t lx;
 122
 123   EXTRACT_WORDS (hx, lx, x);
 124
 125   if (__builtin_expect (((hx & 0x7fffffff) | lx) == 0, 0))
 126     {
 127       /* Return value for x == 0 is Inf with divide by zero exception.  */
 128       *signgamp = 0;
 129       return 1.0 / x;
 130     }
 131   if (__builtin_expect (hx < 0, 0)
 132       && (u_int32_t) hx < 0xfff00000 && __rint (x) == x)
 133     {
 134       /* Return value for integer x < 0 is NaN with invalid exception.  */
 135       *signgamp = 0;
 136       return (x - x) / (x - x);
 137     }
 138   if (__builtin_expect ((unsigned int) hx == 0xfff00000 && lx == 0, 0))
 139     {
 140       /* x == -Inf.  According to ISO this is NaN.  */
 141       *signgamp = 0;
 142       return x - x;
 143     }
 144   if (__builtin_expect ((hx & 0x7ff00000) == 0x7ff00000, 0))
 145     {
 146       /* Positive infinity (return positive infinity) or NaN (return
 147          NaN).  */
 148       *signgamp = 0;
 149       return x + x;
 150     }
 151
 152   if (x >= 172.0)
 153     {
 154       /* Overflow.  */
 155       *signgamp = 0;
 156       return DBL_MAX * DBL_MAX;
 157     }
 158   else if (x > 0.0)
 159     {
 160       *signgamp = 0;
 161       int exp2_adj;
 162       double ret = gamma_positive (x, &exp2_adj);
 163       return __scalbn (ret, exp2_adj);
 164     }
 165   else if (x >= -DBL_EPSILON / 4.0)
 166     {
 167       *signgamp = 0;
 168       return 1.0 / x;
 169     }
 170   else
 171     {
 172       double tx = __trunc (x);
 173       *signgamp = (tx == 2.0 * __trunc (tx / 2.0)) ? -1 : 1;
 174       if (x <= -184.0)
 175         /* Underflow.  */
 176         return DBL_MIN * DBL_MIN;
 177       double frac = tx - x;
 178       if (frac > 0.5)
 179         frac = 1.0 - frac;
 180       double sinpix = (frac <= 0.25
 181                        ? __sin (M_PI * frac)
 182                        : __cos (M_PI * (0.5 - frac)));
 183       int exp2_adj;
 184       double ret = M_PI / (-x * sinpix * gamma_positive (-x, &exp2_adj));
 185       return __scalbn (ret, -exp2_adj);
 186     }
 187 }
 188 strong_alias (__ieee754_gamma_r, __gamma_r_finite)