sysdeps/ieee754/ldbl-128ibm/e_gammal_r.c

   1 /* Implementation of gamma function according to ISO C.
   2    Copyright (C) 1997-2020 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    <https://www.gnu.org/licenses/>.  */
  20
  21 #include <math.h>
  22 #include <math_private.h>
  23 #include <fenv_private.h>
  24 #include <math-underflow.h>
  25 #include <float.h>
  26
  27 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
  28    approximation to gamma function.  */
  29
  30 static const long double gamma_coeff[] =
  31   {
  32     0x1.555555555555555555555555558p-4L,
  33     -0xb.60b60b60b60b60b60b60b60b6p-12L,
  34     0x3.4034034034034034034034034p-12L,
  35     -0x2.7027027027027027027027027p-12L,
  36     0x3.72a3c5631fe46ae1d4e700dca9p-12L,
  37     -0x7.daac36664f1f207daac36664f2p-12L,
  38     0x1.a41a41a41a41a41a41a41a41a4p-8L,
  39     -0x7.90a1b2c3d4e5f708192a3b4c5ep-8L,
  40     0x2.dfd2c703c0cfff430edfd2c704p-4L,
  41     -0x1.6476701181f39edbdb9ce625988p+0L,
  42     0xd.672219167002d3a7a9c886459cp+0L,
  43     -0x9.cd9292e6660d55b3f712eb9e08p+4L,
  44     0x8.911a740da740da740da740da74p+8L,
  45   };
  46
  47 #define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0]))
  48
  49 /* Return gamma (X), for positive X less than 191, in the form R *
  50    2^(*EXP2_ADJ), where R is the return value and *EXP2_ADJ is set to
  51    avoid overflow or underflow in intermediate calculations.  */
  52
  53 static long double
  54 gammal_positive (long double x, int *exp2_adj)
  55 {
  56   int local_signgam;
  57   if (x < 0.5L)
  58     {
  59       *exp2_adj = 0;
  60       return __ieee754_expl (__ieee754_lgammal_r (x + 1, &local_signgam)) / x;
  61     }
  62   else if (x <= 1.5L)
  63     {
  64       *exp2_adj = 0;
  65       return __ieee754_expl (__ieee754_lgammal_r (x, &local_signgam));
  66     }
  67   else if (x < 11.5L)
  68     {
  69       /* Adjust into the range for using exp (lgamma).  */
  70       *exp2_adj = 0;
  71       long double n = ceill (x - 1.5L);
  72       long double x_adj = x - n;
  73       long double eps;
  74       long double prod = __gamma_productl (x_adj, 0, n, &eps);
  75       return (__ieee754_expl (__ieee754_lgammal_r (x_adj, &local_signgam))
  76               * prod * (1.0L + eps));
  77     }
  78   else
  79     {
  80       long double eps = 0;
  81       long double x_eps = 0;
  82       long double x_adj = x;
  83       long double prod = 1;
  84       if (x < 23.0L)
  85         {
  86           /* Adjust into the range for applying Stirling's
  87              approximation.  */
  88           long double n = ceill (23.0L - x);
  89           x_adj = x + n;
  90           x_eps = (x - (x_adj - n));
  91           prod = __gamma_productl (x_adj - n, x_eps, n, &eps);
  92         }
  93       /* The result is now gamma (X_ADJ + X_EPS) / (PROD * (1 + EPS)).
  94          Compute gamma (X_ADJ + X_EPS) using Stirling's approximation,
  95          starting by computing pow (X_ADJ, X_ADJ) with a power of 2
  96          factored out.  */
  97       long double exp_adj = -eps;
  98       long double x_adj_int = roundl (x_adj);
  99       long double x_adj_frac = x_adj - x_adj_int;
 100       int x_adj_log2;
 101       long double x_adj_mant = __frexpl (x_adj, &x_adj_log2);
 102       if (x_adj_mant < M_SQRT1_2l)
 103         {
 104           x_adj_log2--;
 105           x_adj_mant *= 2.0L;
 106         }
 107       *exp2_adj = x_adj_log2 * (int) x_adj_int;
 108       long double ret = (__ieee754_powl (x_adj_mant, x_adj)
 109                          * __ieee754_exp2l (x_adj_log2 * x_adj_frac)
 110                          * __ieee754_expl (-x_adj)
 111                          * sqrtl (2 * M_PIl / x_adj)
 112                          / prod);
 113       exp_adj += x_eps * __ieee754_logl (x_adj);
 114       long double bsum = gamma_coeff[NCOEFF - 1];
 115       long double x_adj2 = x_adj * x_adj;
 116       for (size_t i = 1; i <= NCOEFF - 1; i++)
 117         bsum = bsum / x_adj2 + gamma_coeff[NCOEFF - 1 - i];
 118       exp_adj += bsum / x_adj;
 119       return ret + ret * __expm1l (exp_adj);
 120     }
 121 }
 122
 123 long double
 124 __ieee754_gammal_r (long double x, int *signgamp)
 125 {
 126   int64_t hx;
 127   double xhi;
 128   long double ret;
 129
 130   xhi = ldbl_high (x);
 131   EXTRACT_WORDS64 (hx, xhi);
 132
 133   if ((hx & 0x7fffffffffffffffLL) == 0)
 134     {
 135       /* Return value for x == 0 is Inf with divide by zero exception.  */
 136       *signgamp = 0;
 137       return 1.0 / x;
 138     }
 139   if (hx < 0 && (uint64_t) hx < 0xfff0000000000000ULL && rintl (x) == x)
 140     {
 141       /* Return value for integer x < 0 is NaN with invalid exception.  */
 142       *signgamp = 0;
 143       return (x - x) / (x - x);
 144     }
 145   if (hx == 0xfff0000000000000ULL)
 146     {
 147       /* x == -Inf.  According to ISO this is NaN.  */
 148       *signgamp = 0;
 149       return x - x;
 150     }
 151   if ((hx & 0x7ff0000000000000ULL) == 0x7ff0000000000000ULL)
 152     {
 153       /* Positive infinity (return positive infinity) or NaN (return
 154          NaN).  */
 155       *signgamp = 0;
 156       return x + x;
 157     }
 158
 159   if (x >= 172.0L)
 160     {
 161       /* Overflow.  */
 162       *signgamp = 0;
 163       return LDBL_MAX * LDBL_MAX;
 164     }
 165   else
 166     {
 167       SET_RESTORE_ROUNDL (FE_TONEAREST);
 168       if (x > 0.0L)
 169         {
 170           *signgamp = 0;
 171           int exp2_adj;
 172           ret = gammal_positive (x, &exp2_adj);
 173           ret = __scalbnl (ret, exp2_adj);
 174         }
 175       else if (x >= -0x1p-110L)
 176         {
 177           *signgamp = 0;
 178           ret = 1.0L / x;
 179         }
 180       else
 181         {
 182           long double tx = truncl (x);
 183           *signgamp = (tx == 2.0L * truncl (tx / 2.0L)) ? -1 : 1;
 184           if (x <= -191.0L)
 185             /* Underflow.  */
 186             ret = LDBL_MIN * LDBL_MIN;
 187           else
 188             {
 189               long double frac = tx - x;
 190               if (frac > 0.5L)
 191                 frac = 1.0L - frac;
 192               long double sinpix = (frac <= 0.25L
 193                                     ? __sinl (M_PIl * frac)
 194                                     : __cosl (M_PIl * (0.5L - frac)));
 195               int exp2_adj;
 196               ret = M_PIl / (-x * sinpix
 197                              * gammal_positive (-x, &exp2_adj));
 198               ret = __scalbnl (ret, -exp2_adj);
 199               math_check_force_underflow_nonneg (ret);
 200             }
 201         }
 202     }
 203   if (isinf (ret) && x != 0)
 204     {
 205       if (*signgamp < 0)
 206         return -(-copysignl (LDBL_MAX, ret) * LDBL_MAX);
 207       else
 208         return copysignl (LDBL_MAX, ret) * LDBL_MAX;
 209     }
 210   else if (ret == 0)
 211     {
 212       if (*signgamp < 0)
 213         return -(-copysignl (LDBL_MIN, ret) * LDBL_MIN);
 214       else
 215         return copysignl (LDBL_MIN, ret) * LDBL_MIN;
 216     }
 217   else
 218     return ret;
 219 }
 220 strong_alias (__ieee754_gammal_r, __gammal_r_finite)