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