sysdeps/ieee754/ldbl-128/e_gammal_r.c

   1 /* Implementation of gamma function according to ISO C.
   2    Copyright (C) 1997-2024 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4
   5    The GNU C Library is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU Lesser General Public
   7    License as published by the Free Software Foundation; either
   8    version 2.1 of the License, or (at your option) any later version.
   9
  10    The GNU C Library is distributed in the hope that it will be useful,
  11    but WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    Lesser General Public License for more details.
  14
  15    You should have received a copy of the GNU Lesser General Public
  16    License along with the GNU C Library; if not, see
  17    <https://www.gnu.org/licenses/>.  */
  18
  19 #include <math.h>
  20 #include <math_private.h>
  21 #include <fenv_private.h>
  22 #include <math-underflow.h>
  23 #include <float.h>
  24 #include <libm-alias-finite.h>
  25
  26 /* Coefficients B_2k / 2k(2k-1) of x^-(2k-1) inside exp in Stirling's
  27    approximation to gamma function.  */
  28
  29 static const _Float128 gamma_coeff[] =
  30   {
  31     L(0x1.5555555555555555555555555555p-4),
  32     L(-0xb.60b60b60b60b60b60b60b60b60b8p-12),
  33     L(0x3.4034034034034034034034034034p-12),
  34     L(-0x2.7027027027027027027027027028p-12),
  35     L(0x3.72a3c5631fe46ae1d4e700dca8f2p-12),
  36     L(-0x7.daac36664f1f207daac36664f1f4p-12),
  37     L(0x1.a41a41a41a41a41a41a41a41a41ap-8),
  38     L(-0x7.90a1b2c3d4e5f708192a3b4c5d7p-8),
  39     L(0x2.dfd2c703c0cfff430edfd2c703cp-4),
  40     L(-0x1.6476701181f39edbdb9ce625987dp+0),
  41     L(0xd.672219167002d3a7a9c886459cp+0),
  42     L(-0x9.cd9292e6660d55b3f712eb9e07c8p+4),
  43     L(0x8.911a740da740da740da740da741p+8),
  44     L(-0x8.d0cc570e255bf59ff6eec24b49p+12),
  45   };
  46
  47 #define NCOEFF (sizeof (gamma_coeff) / sizeof (gamma_coeff[0]))
  48
  49 /* Return gamma (X), for positive X less than 1775, 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 _Float128
  54 gammal_positive (_Float128 x, int *exp2_adj)
  55 {
  56   int local_signgam;
  57   if (x < L(0.5))
  58     {
  59       *exp2_adj = 0;
  60       return __ieee754_expl (__ieee754_lgammal_r (x + 1, &local_signgam)) / x;
  61     }
  62   else if (x <= L(1.5))
  63     {
  64       *exp2_adj = 0;
  65       return __ieee754_expl (__ieee754_lgammal_r (x, &local_signgam));
  66     }
  67   else if (x < L(12.5))
  68     {
  69       /* Adjust into the range for using exp (lgamma).  */
  70       *exp2_adj = 0;
  71       _Float128 n = ceill (x - L(1.5));
  72       _Float128 x_adj = x - n;
  73       _Float128 eps;
  74       _Float128 prod = __gamma_productl (x_adj, 0, n, &eps);
  75       return (__ieee754_expl (__ieee754_lgammal_r (x_adj, &local_signgam))
  76               * prod * (1 + eps));
  77     }
  78   else
  79     {
  80       _Float128 eps = 0;
  81       _Float128 x_eps = 0;
  82       _Float128 x_adj = x;
  83       _Float128 prod = 1;
  84       if (x < 24)
  85         {
  86           /* Adjust into the range for applying Stirling's
  87              approximation.  */
  88           _Float128 n = ceill (24 - 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       _Float128 exp_adj = -eps;
  98       _Float128 x_adj_int = roundl (x_adj);
  99       _Float128 x_adj_frac = x_adj - x_adj_int;
 100       int x_adj_log2;
 101       _Float128 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;
 106         }
 107       *exp2_adj = x_adj_log2 * (int) x_adj_int;
 108       _Float128 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       _Float128 bsum = gamma_coeff[NCOEFF - 1];
 115       _Float128 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 _Float128
 124 __ieee754_gammal_r (_Float128 x, int *signgamp)
 125 {
 126   int64_t hx;
 127   uint64_t lx;
 128   _Float128 ret;
 129
 130   GET_LDOUBLE_WORDS64 (hx, lx, x);
 131
 132   if (((hx & 0x7fffffffffffffffLL) | lx) == 0)
 133     {
 134       /* Return value for x == 0 is Inf with divide by zero exception.  */
 135       *signgamp = 0;
 136       return 1.0 / x;
 137     }
 138   if (hx < 0 && (uint64_t) hx < 0xffff000000000000ULL && rintl (x) == x)
 139     {
 140       /* Return value for integer x < 0 is NaN with invalid exception.  */
 141       *signgamp = 0;
 142       return (x - x) / (x - x);
 143     }
 144   if (hx == 0xffff000000000000ULL && lx == 0)
 145     {
 146       /* x == -Inf.  According to ISO this is NaN.  */
 147       *signgamp = 0;
 148       return x - x;
 149     }
 150   if ((hx & 0x7fff000000000000ULL) == 0x7fff000000000000ULL)
 151     {
 152       /* Positive infinity (return positive infinity) or NaN (return
 153          NaN).  */
 154       *signgamp = 0;
 155       return x + x;
 156     }
 157
 158   if (x >= 1756)
 159     {
 160       /* Overflow.  */
 161       *signgamp = 0;
 162       return LDBL_MAX * LDBL_MAX;
 163     }
 164   else
 165     {
 166       SET_RESTORE_ROUNDL (FE_TONEAREST);
 167       if (x > 0)
 168         {
 169           *signgamp = 0;
 170           int exp2_adj;
 171           ret = gammal_positive (x, &exp2_adj);
 172           ret = __scalbnl (ret, exp2_adj);
 173         }
 174       else if (x >= -LDBL_EPSILON / 4)
 175         {
 176           *signgamp = 0;
 177           ret = 1 / x;
 178         }
 179       else
 180         {
 181           _Float128 tx = truncl (x);
 182           *signgamp = (tx == 2 * truncl (tx / 2)) ? -1 : 1;
 183           if (x <= -1775)
 184             /* Underflow.  */
 185             ret = LDBL_MIN * LDBL_MIN;
 186           else
 187             {
 188               _Float128 frac = tx - x;
 189               if (frac > L(0.5))
 190                 frac = 1 - frac;
 191               _Float128 sinpix = (frac <= L(0.25)
 192                                   ? __sinl (M_PIl * frac)
 193                                   : __cosl (M_PIl * (L(0.5) - frac)));
 194               int exp2_adj;
 195               ret = M_PIl / (-x * sinpix
 196                              * gammal_positive (-x, &exp2_adj));
 197               ret = __scalbnl (ret, -exp2_adj);
 198               math_check_force_underflow_nonneg (ret);
 199             }
 200         }
 201     }
 202   if (isinf (ret) && x != 0)
 203     {
 204       if (*signgamp < 0)
 205         return -(-copysignl (LDBL_MAX, ret) * LDBL_MAX);
 206       else
 207         return copysignl (LDBL_MAX, ret) * LDBL_MAX;
 208     }
 209   else if (ret == 0)
 210     {
 211       if (*signgamp < 0)
 212         return -(-copysignl (LDBL_MIN, ret) * LDBL_MIN);
 213       else
 214         return copysignl (LDBL_MIN, ret) * LDBL_MIN;
 215     }
 216   else
 217     return ret;
 218 }
 219 libm_alias_finite (__ieee754_gammal_r, __gammal_r)