sysdeps/ieee754/ldbl-96/e_lgammal_r.c

   1 /*
   2  * ====================================================
   3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   4  *
   5  * Developed at SunPro, a Sun Microsystems, Inc. business.
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  10  */
  11
  12 /* Long double expansions are
  13   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
  14   and are incorporated herein by permission of the author.  The author
  15   reserves the right to distribute this material elsewhere under different
  16   copying permissions.  These modifications are distributed here under
  17   the following terms:
  18
  19     This library is free software; you can redistribute it and/or
  20     modify it under the terms of the GNU Lesser General Public
  21     License as published by the Free Software Foundation; either
  22     version 2.1 of the License, or (at your option) any later version.
  23
  24     This library is distributed in the hope that it will be useful,
  25     but WITHOUT ANY WARRANTY; without even the implied warranty of
  26     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  27     Lesser General Public License for more details.
  28
  29     You should have received a copy of the GNU Lesser General Public
  30     License along with this library; if not, see
  31     <https://www.gnu.org/licenses/>.  */
  32
  33 /* __ieee754_lgammal_r(x, signgamp)
  34  * Reentrant version of the logarithm of the Gamma function
  35  * with user provide pointer for the sign of Gamma(x).
  36  *
  37  * Method:
  38  *   1. Argument Reduction for 0 < x <= 8
  39  *      Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
  40  *      reduce x to a number in [1.5,2.5] by
  41  *              lgamma(1+s) = log(s) + lgamma(s)
  42  *      for example,
  43  *              lgamma(7.3) = log(6.3) + lgamma(6.3)
  44  *                          = log(6.3*5.3) + lgamma(5.3)
  45  *                          = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
  46  *   2. Polynomial approximation of lgamma around its
  47  *      minimum ymin=1.461632144968362245 to maintain monotonicity.
  48  *      On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
  49  *              Let z = x-ymin;
  50  *              lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
  51  *   2. Rational approximation in the primary interval [2,3]
  52  *      We use the following approximation:
  53  *              s = x-2.0;
  54  *              lgamma(x) = 0.5*s + s*P(s)/Q(s)
  55  *      Our algorithms are based on the following observation
  56  *
  57  *                             zeta(2)-1    2    zeta(3)-1    3
  58  * lgamma(2+s) = s*(1-Euler) + --------- * s  -  --------- * s  + ...
  59  *                                 2                 3
  60  *
  61  *      where Euler = 0.5771... is the Euler constant, which is very
  62  *      close to 0.5.
  63  *
  64  *   3. For x>=8, we have
  65  *      lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
  66  *      (better formula:
  67  *         lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
  68  *      Let z = 1/x, then we approximation
  69  *              f(z) = lgamma(x) - (x-0.5)(log(x)-1)
  70  *      by
  71  *                                  3       5             11
  72  *              w = w0 + w1*z + w2*z  + w3*z  + ... + w6*z
  73  *
  74  *   4. For negative x, since (G is gamma function)
  75  *              -x*G(-x)*G(x) = pi/sin(pi*x),
  76  *      we have
  77  *              G(x) = pi/(sin(pi*x)*(-x)*G(-x))
  78  *      since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
  79  *      Hence, for x<0, signgam = sign(sin(pi*x)) and
  80  *              lgamma(x) = log(|Gamma(x)|)
  81  *                        = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
  82  *      Note: one should avoid compute pi*(-x) directly in the
  83  *            computation of sin(pi*(-x)).
  84  *
  85  *   5. Special Cases
  86  *              lgamma(2+s) ~ s*(1-Euler) for tiny s
  87  *              lgamma(1)=lgamma(2)=0
  88  *              lgamma(x) ~ -log(x) for tiny x
  89  *              lgamma(0) = lgamma(inf) = inf
  90  *              lgamma(-integer) = +-inf
  91  *
  92  */
  93
  94 #include <math.h>
  95 #include <math_private.h>
  96 #include <libc-diag.h>
  97 #include <libm-alias-finite.h>
  98
  99 static const long double
 100   half = 0.5L,
 101   one = 1.0L,
 102   pi = 3.14159265358979323846264L,
 103   two63 = 9.223372036854775808e18L,
 104
 105   /* lgam(1+x) = 0.5 x + x a(x)/b(x)
 106      -0.268402099609375 <= x <= 0
 107      peak relative error 6.6e-22 */
 108   a0 = -6.343246574721079391729402781192128239938E2L,
 109   a1 =  1.856560238672465796768677717168371401378E3L,
 110   a2 =  2.404733102163746263689288466865843408429E3L,
 111   a3 =  8.804188795790383497379532868917517596322E2L,
 112   a4 =  1.135361354097447729740103745999661157426E2L,
 113   a5 =  3.766956539107615557608581581190400021285E0L,
 114
 115   b0 =  8.214973713960928795704317259806842490498E3L,
 116   b1 =  1.026343508841367384879065363925870888012E4L,
 117   b2 =  4.553337477045763320522762343132210919277E3L,
 118   b3 =  8.506975785032585797446253359230031874803E2L,
 119   b4 =  6.042447899703295436820744186992189445813E1L,
 120   /* b5 =  1.000000000000000000000000000000000000000E0 */
 121
 122
 123   tc =  1.4616321449683623412626595423257213284682E0L,
 124   tf = -1.2148629053584961146050602565082954242826E-1,/* double precision */
 125 /* tt = (tail of tf), i.e. tf + tt has extended precision. */
 126   tt = 3.3649914684731379602768989080467587736363E-18L,
 127   /* lgam ( 1.4616321449683623412626595423257213284682E0 ) =
 128 -1.2148629053584960809551455717769158215135617312999903886372437313313530E-1 */
 129
 130   /* lgam (x + tc) = tf + tt + x g(x)/h(x)
 131      - 0.230003726999612341262659542325721328468 <= x
 132      <= 0.2699962730003876587373404576742786715318
 133      peak relative error 2.1e-21 */
 134   g0 = 3.645529916721223331888305293534095553827E-18L,
 135   g1 = 5.126654642791082497002594216163574795690E3L,
 136   g2 = 8.828603575854624811911631336122070070327E3L,
 137   g3 = 5.464186426932117031234820886525701595203E3L,
 138   g4 = 1.455427403530884193180776558102868592293E3L,
 139   g5 = 1.541735456969245924860307497029155838446E2L,
 140   g6 = 4.335498275274822298341872707453445815118E0L,
 141
 142   h0 = 1.059584930106085509696730443974495979641E4L,
 143   h1 =  2.147921653490043010629481226937850618860E4L,
 144   h2 = 1.643014770044524804175197151958100656728E4L,
 145   h3 =  5.869021995186925517228323497501767586078E3L,
 146   h4 =  9.764244777714344488787381271643502742293E2L,
 147   h5 =  6.442485441570592541741092969581997002349E1L,
 148   /* h6 = 1.000000000000000000000000000000000000000E0 */
 149
 150
 151   /* lgam (x+1) = -0.5 x + x u(x)/v(x)
 152      -0.100006103515625 <= x <= 0.231639862060546875
 153      peak relative error 1.3e-21 */
 154   u0 = -8.886217500092090678492242071879342025627E1L,
 155   u1 =  6.840109978129177639438792958320783599310E2L,
 156   u2 =  2.042626104514127267855588786511809932433E3L,
 157   u3 =  1.911723903442667422201651063009856064275E3L,
 158   u4 =  7.447065275665887457628865263491667767695E2L,
 159   u5 =  1.132256494121790736268471016493103952637E2L,
 160   u6 =  4.484398885516614191003094714505960972894E0L,
 161
 162   v0 =  1.150830924194461522996462401210374632929E3L,
 163   v1 =  3.399692260848747447377972081399737098610E3L,
 164   v2 =  3.786631705644460255229513563657226008015E3L,
 165   v3 =  1.966450123004478374557778781564114347876E3L,
 166   v4 =  4.741359068914069299837355438370682773122E2L,
 167   v5 =  4.508989649747184050907206782117647852364E1L,
 168   /* v6 =  1.000000000000000000000000000000000000000E0 */
 169
 170
 171   /* lgam (x+2) = .5 x + x s(x)/r(x)
 172      0 <= x <= 1
 173      peak relative error 7.2e-22 */
 174   s0 =  1.454726263410661942989109455292824853344E6L,
 175   s1 = -3.901428390086348447890408306153378922752E6L,
 176   s2 = -6.573568698209374121847873064292963089438E6L,
 177   s3 = -3.319055881485044417245964508099095984643E6L,
 178   s4 = -7.094891568758439227560184618114707107977E5L,
 179   s5 = -6.263426646464505837422314539808112478303E4L,
 180   s6 = -1.684926520999477529949915657519454051529E3L,
 181
 182   r0 = -1.883978160734303518163008696712983134698E7L,
 183   r1 = -2.815206082812062064902202753264922306830E7L,
 184   r2 = -1.600245495251915899081846093343626358398E7L,
 185   r3 = -4.310526301881305003489257052083370058799E6L,
 186   r4 = -5.563807682263923279438235987186184968542E5L,
 187   r5 = -3.027734654434169996032905158145259713083E4L,
 188   r6 = -4.501995652861105629217250715790764371267E2L,
 189   /* r6 =  1.000000000000000000000000000000000000000E0 */
 190
 191
 192 /* lgam(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x w(1/x^2)
 193    x >= 8
 194    Peak relative error 1.51e-21
 195    w0 = LS2PI - 0.5 */
 196   w0 =  4.189385332046727417803e-1L,
 197   w1 =  8.333333333333331447505E-2L,
 198   w2 = -2.777777777750349603440E-3L,
 199   w3 =  7.936507795855070755671E-4L,
 200   w4 = -5.952345851765688514613E-4L,
 201   w5 =  8.412723297322498080632E-4L,
 202   w6 = -1.880801938119376907179E-3L,
 203   w7 =  4.885026142432270781165E-3L;
 204
 205 static const long double zero = 0.0L;
 206
 207 static long double
 208 sin_pi (long double x)
 209 {
 210   long double y, z;
 211   int n, ix;
 212   uint32_t se, i0, i1;
 213
 214   GET_LDOUBLE_WORDS (se, i0, i1, x);
 215   ix = se & 0x7fff;
 216   ix = (ix << 16) | (i0 >> 16);
 217   if (ix < 0x3ffd8000) /* 0.25 */
 218     return __sinl (pi * x);
 219   y = -x;                       /* x is assume negative */
 220
 221   /*
 222    * argument reduction, make sure inexact flag not raised if input
 223    * is an integer
 224    */
 225   z = floorl (y);
 226   if (z != y)
 227     {                           /* inexact anyway */
 228       y  *= 0.5;
 229       y = 2.0*(y - floorl(y));          /* y = |x| mod 2.0 */
 230       n = (int) (y*4.0);
 231     }
 232   else
 233     {
 234       if (ix >= 0x403f8000)  /* 2^64 */
 235         {
 236           y = zero; n = 0;                 /* y must be even */
 237         }
 238       else
 239         {
 240         if (ix < 0x403e8000)  /* 2^63 */
 241           z = y + two63;        /* exact */
 242         GET_LDOUBLE_WORDS (se, i0, i1, z);
 243         n = i1 & 1;
 244         y  = n;
 245         n <<= 2;
 246       }
 247     }
 248
 249   switch (n)
 250     {
 251     case 0:
 252       y = __sinl (pi * y);
 253       break;
 254     case 1:
 255     case 2:
 256       y = __cosl (pi * (half - y));
 257       break;
 258     case 3:
 259     case 4:
 260       y = __sinl (pi * (one - y));
 261       break;
 262     case 5:
 263     case 6:
 264       y = -__cosl (pi * (y - 1.5));
 265       break;
 266     default:
 267       y = __sinl (pi * (y - 2.0));
 268       break;
 269     }
 270   return -y;
 271 }
 272
 273
 274 long double
 275 __ieee754_lgammal_r (long double x, int *signgamp)
 276 {
 277   long double t, y, z, nadj, p, p1, p2, q, r, w;
 278   int i, ix;
 279   uint32_t se, i0, i1;
 280
 281   *signgamp = 1;
 282   GET_LDOUBLE_WORDS (se, i0, i1, x);
 283   ix = se & 0x7fff;
 284
 285   if (__builtin_expect((ix | i0 | i1) == 0, 0))
 286     {
 287       if (se & 0x8000)
 288         *signgamp = -1;
 289       return one / fabsl (x);
 290     }
 291
 292   ix = (ix << 16) | (i0 >> 16);
 293
 294   /* purge off +-inf, NaN, +-0, and negative arguments */
 295   if (__builtin_expect(ix >= 0x7fff0000, 0))
 296     return x * x;
 297
 298   if (__builtin_expect(ix < 0x3fc08000, 0)) /* 2^-63 */
 299     {                           /* |x|<2**-63, return -log(|x|) */
 300       if (se & 0x8000)
 301         {
 302           *signgamp = -1;
 303           return -__ieee754_logl (-x);
 304         }
 305       else
 306         return -__ieee754_logl (x);
 307     }
 308   if (se & 0x8000)
 309     {
 310       if (x < -2.0L && x > -33.0L)
 311         return __lgamma_negl (x, signgamp);
 312       t = sin_pi (x);
 313       if (t == zero)
 314         return one / fabsl (t); /* -integer */
 315       nadj = __ieee754_logl (pi / fabsl (t * x));
 316       if (t < zero)
 317         *signgamp = -1;
 318       x = -x;
 319     }
 320
 321   /* purge off 1 and 2 */
 322   if ((((ix - 0x3fff8000) | i0 | i1) == 0)
 323       || (((ix - 0x40008000) | i0 | i1) == 0))
 324     r = 0;
 325   else if (ix < 0x40008000) /* 2.0 */
 326     {
 327       /* x < 2.0 */
 328       if (ix <= 0x3ffee666) /* 8.99993896484375e-1 */
 329         {
 330           /* lgamma(x) = lgamma(x+1) - log(x) */
 331           r = -__ieee754_logl (x);
 332           if (ix >= 0x3ffebb4a) /* 7.31597900390625e-1 */
 333             {
 334               y = x - one;
 335               i = 0;
 336             }
 337           else if (ix >= 0x3ffced33)/* 2.31639862060546875e-1 */
 338             {
 339               y = x - (tc - one);
 340               i = 1;
 341             }
 342           else
 343             {
 344               /* x < 0.23 */
 345               y = x;
 346               i = 2;
 347             }
 348         }
 349       else
 350         {
 351           r = zero;
 352           if (ix >= 0x3fffdda6) /* 1.73162841796875 */
 353             {
 354               /* [1.7316,2] */
 355               y = x - 2.0;
 356               i = 0;
 357             }
 358           else if (ix >= 0x3fff9da6)/* 1.23162841796875 */
 359             {
 360               /* [1.23,1.73] */
 361               y = x - tc;
 362               i = 1;
 363             }
 364           else
 365             {
 366               /* [0.9, 1.23] */
 367               y = x - one;
 368               i = 2;
 369             }
 370         }
 371       switch (i)
 372         {
 373         case 0:
 374           p1 = a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * a5))));
 375           p2 = b0 + y * (b1 + y * (b2 + y * (b3 + y * (b4 + y))));
 376           r += half * y + y * p1/p2;
 377           break;
 378         case 1:
 379     p1 = g0 + y * (g1 + y * (g2 + y * (g3 + y * (g4 + y * (g5 + y * g6)))));
 380     p2 = h0 + y * (h1 + y * (h2 + y * (h3 + y * (h4 + y * (h5 + y)))));
 381     p = tt + y * p1/p2;
 382           r += (tf + p);
 383           break;
 384         case 2:
 385  p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * (u5 + y * u6))))));
 386       p2 = v0 + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * (v5 + y)))));
 387           r += (-half * y + p1 / p2);
 388         }
 389     }
 390   else if (ix < 0x40028000) /* 8.0 */
 391     {
 392       /* x < 8.0 */
 393       i = (int) x;
 394       t = zero;
 395       y = x - (double) i;
 396   p = y *
 397      (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
 398   q = r0 + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * (r6 + y))))));
 399       r = half * y + p / q;
 400       z = one;                  /* lgamma(1+s) = log(s) + lgamma(s) */
 401       switch (i)
 402         {
 403         case 7:
 404           z *= (y + 6.0);       /* FALLTHRU */
 405         case 6:
 406           z *= (y + 5.0);       /* FALLTHRU */
 407         case 5:
 408           z *= (y + 4.0);       /* FALLTHRU */
 409         case 4:
 410           z *= (y + 3.0);       /* FALLTHRU */
 411         case 3:
 412           z *= (y + 2.0);       /* FALLTHRU */
 413           r += __ieee754_logl (z);
 414           break;
 415         }
 416     }
 417   else if (ix < 0x40418000) /* 2^66 */
 418     {
 419       /* 8.0 <= x < 2**66 */
 420       t = __ieee754_logl (x);
 421       z = one / x;
 422       y = z * z;
 423       w = w0 + z * (w1
 424           + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * (w6 + y * w7))))));
 425       r = (x - half) * (t - one) + w;
 426     }
 427   else
 428     /* 2**66 <= x <= inf */
 429     r = x * (__ieee754_logl (x) - one);
 430   /* NADJ is set for negative arguments but not otherwise, resulting
 431      in warnings that it may be used uninitialized although in the
 432      cases where it is used it has always been set.  */
 433   DIAG_PUSH_NEEDS_COMMENT;
 434   DIAG_IGNORE_NEEDS_COMMENT (4.9, "-Wmaybe-uninitialized");
 435   if (se & 0x8000)
 436     r = nadj - r;
 437   DIAG_POP_NEEDS_COMMENT;
 438   return r;
 439 }
 440 libm_alias_finite (__ieee754_lgammal_r, __lgammal_r)