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, write to the Free Software
  31     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
  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  *      minimun 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
  97 #ifdef __STDC__
  98 static const long double
  99 #else
 100 static long double
 101 #endif
 102   half = 0.5L,
 103   one = 1.0L,
 104   pi = 3.14159265358979323846264L,
 105   two63 = 9.223372036854775808e18L,
 106
 107   /* lgam(1+x) = 0.5 x + x a(x)/b(x)
 108      -0.268402099609375 <= x <= 0
 109      peak relative error 6.6e-22 */
 110   a0 = -6.343246574721079391729402781192128239938E2L,
 111   a1 =  1.856560238672465796768677717168371401378E3L,
 112   a2 =  2.404733102163746263689288466865843408429E3L,
 113   a3 =  8.804188795790383497379532868917517596322E2L,
 114   a4 =  1.135361354097447729740103745999661157426E2L,
 115   a5 =  3.766956539107615557608581581190400021285E0L,
 116
 117   b0 =  8.214973713960928795704317259806842490498E3L,
 118   b1 =  1.026343508841367384879065363925870888012E4L,
 119   b2 =  4.553337477045763320522762343132210919277E3L,
 120   b3 =  8.506975785032585797446253359230031874803E2L,
 121   b4 =  6.042447899703295436820744186992189445813E1L,
 122   /* b5 =  1.000000000000000000000000000000000000000E0 */
 123
 124
 125   tc =  1.4616321449683623412626595423257213284682E0L,
 126   tf = -1.2148629053584961146050602565082954242826E-1,/* double precision */
 127 /* tt = (tail of tf), i.e. tf + tt has extended precision. */
 128   tt = 3.3649914684731379602768989080467587736363E-18L,
 129   /* lgam ( 1.4616321449683623412626595423257213284682E0 ) =
 130 -1.2148629053584960809551455717769158215135617312999903886372437313313530E-1 */
 131
 132   /* lgam (x + tc) = tf + tt + x g(x)/h(x)
 133      - 0.230003726999612341262659542325721328468 <= x
 134      <= 0.2699962730003876587373404576742786715318
 135      peak relative error 2.1e-21 */
 136   g0 = 3.645529916721223331888305293534095553827E-18L,
 137   g1 = 5.126654642791082497002594216163574795690E3L,
 138   g2 = 8.828603575854624811911631336122070070327E3L,
 139   g3 = 5.464186426932117031234820886525701595203E3L,
 140   g4 = 1.455427403530884193180776558102868592293E3L,
 141   g5 = 1.541735456969245924860307497029155838446E2L,
 142   g6 = 4.335498275274822298341872707453445815118E0L,
 143
 144   h0 = 1.059584930106085509696730443974495979641E4L,
 145   h1 =  2.147921653490043010629481226937850618860E4L,
 146   h2 = 1.643014770044524804175197151958100656728E4L,
 147   h3 =  5.869021995186925517228323497501767586078E3L,
 148   h4 =  9.764244777714344488787381271643502742293E2L,
 149   h5 =  6.442485441570592541741092969581997002349E1L,
 150   /* h6 = 1.000000000000000000000000000000000000000E0 */
 151
 152
 153   /* lgam (x+1) = -0.5 x + x u(x)/v(x)
 154      -0.100006103515625 <= x <= 0.231639862060546875
 155      peak relative error 1.3e-21 */
 156   u0 = -8.886217500092090678492242071879342025627E1L,
 157   u1 =  6.840109978129177639438792958320783599310E2L,
 158   u2 =  2.042626104514127267855588786511809932433E3L,
 159   u3 =  1.911723903442667422201651063009856064275E3L,
 160   u4 =  7.447065275665887457628865263491667767695E2L,
 161   u5 =  1.132256494121790736268471016493103952637E2L,
 162   u6 =  4.484398885516614191003094714505960972894E0L,
 163
 164   v0 =  1.150830924194461522996462401210374632929E3L,
 165   v1 =  3.399692260848747447377972081399737098610E3L,
 166   v2 =  3.786631705644460255229513563657226008015E3L,
 167   v3 =  1.966450123004478374557778781564114347876E3L,
 168   v4 =  4.741359068914069299837355438370682773122E2L,
 169   v5 =  4.508989649747184050907206782117647852364E1L,
 170   /* v6 =  1.000000000000000000000000000000000000000E0 */
 171
 172
 173   /* lgam (x+2) = .5 x + x s(x)/r(x)
 174      0 <= x <= 1
 175      peak relative error 7.2e-22 */
 176   s0 =  1.454726263410661942989109455292824853344E6L,
 177   s1 = -3.901428390086348447890408306153378922752E6L,
 178   s2 = -6.573568698209374121847873064292963089438E6L,
 179   s3 = -3.319055881485044417245964508099095984643E6L,
 180   s4 = -7.094891568758439227560184618114707107977E5L,
 181   s5 = -6.263426646464505837422314539808112478303E4L,
 182   s6 = -1.684926520999477529949915657519454051529E3L,
 183
 184   r0 = -1.883978160734303518163008696712983134698E7L,
 185   r1 = -2.815206082812062064902202753264922306830E7L,
 186   r2 = -1.600245495251915899081846093343626358398E7L,
 187   r3 = -4.310526301881305003489257052083370058799E6L,
 188   r4 = -5.563807682263923279438235987186184968542E5L,
 189   r5 = -3.027734654434169996032905158145259713083E4L,
 190   r6 = -4.501995652861105629217250715790764371267E2L,
 191   /* r6 =  1.000000000000000000000000000000000000000E0 */
 192
 193
 194 /* lgam(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x w(1/x^2)
 195    x >= 8
 196    Peak relative error 1.51e-21
 197    w0 = LS2PI - 0.5 */
 198   w0 =  4.189385332046727417803e-1L,
 199   w1 =  8.333333333333331447505E-2L,
 200   w2 = -2.777777777750349603440E-3L,
 201   w3 =  7.936507795855070755671E-4L,
 202   w4 = -5.952345851765688514613E-4L,
 203   w5 =  8.412723297322498080632E-4L,
 204   w6 = -1.880801938119376907179E-3L,
 205   w7 =  4.885026142432270781165E-3L;
 206
 207 #ifdef __STDC__
 208 static const long double zero = 0.0L;
 209 #else
 210 static long double zero = 0.0L;
 211 #endif
 212
 213 #ifdef __STDC__
 214 static long double
 215 sin_pi (long double x)
 216 #else
 217 static long double
 218 sin_pi (x)
 219      long double x;
 220 #endif
 221 {
 222   long double y, z;
 223   int n, ix;
 224   u_int32_t se, i0, i1;
 225
 226   GET_LDOUBLE_WORDS (se, i0, i1, x);
 227   ix = se & 0x7fff;
 228   ix = (ix << 16) | (i0 >> 16);
 229   if (ix < 0x3ffd8000) /* 0.25 */
 230     return __sinl (pi * x);
 231   y = -x;                       /* x is assume negative */
 232
 233   /*
 234    * argument reduction, make sure inexact flag not raised if input
 235    * is an integer
 236    */
 237   z = __floorl (y);
 238   if (z != y)
 239     {                           /* inexact anyway */
 240       y  *= 0.5;
 241       y = 2.0*(y - __floorl(y));                /* y = |x| mod 2.0 */
 242       n = (int) (y*4.0);
 243     }
 244   else
 245     {
 246       if (ix >= 0x403f8000)  /* 2^64 */
 247         {
 248           y = zero; n = 0;                 /* y must be even */
 249         }
 250       else
 251         {
 252         if (ix < 0x403e8000)  /* 2^63 */
 253           z = y + two63;        /* exact */
 254         GET_LDOUBLE_WORDS (se, i0, i1, z);
 255         n = i1 & 1;
 256         y  = n;
 257         n <<= 2;
 258       }
 259     }
 260
 261   switch (n)
 262     {
 263     case 0:
 264       y = __sinl (pi * y);
 265       break;
 266     case 1:
 267     case 2:
 268       y = __cosl (pi * (half - y));
 269       break;
 270     case 3:
 271     case 4:
 272       y = __sinl (pi * (one - y));
 273       break;
 274     case 5:
 275     case 6:
 276       y = -__cosl (pi * (y - 1.5));
 277       break;
 278     default:
 279       y = __sinl (pi * (y - 2.0));
 280       break;
 281     }
 282   return -y;
 283 }
 284
 285
 286 #ifdef __STDC__
 287 long double
 288 __ieee754_lgammal_r (long double x, int *signgamp)
 289 #else
 290 long double
 291 __ieee754_lgammal_r (x, signgamp)
 292      long double x;
 293      int *signgamp;
 294 #endif
 295 {
 296   long double t, y, z, nadj, p, p1, p2, q, r, w;
 297   int i, ix;
 298   u_int32_t se, i0, i1;
 299
 300   *signgamp = 1;
 301   GET_LDOUBLE_WORDS (se, i0, i1, x);
 302   ix = se & 0x7fff;
 303
 304   if ((ix | i0 | i1) == 0)
 305     {
 306       if (se & 0x8000)
 307         *signgamp = -1;
 308       return one / fabsl (x);
 309     }
 310
 311   ix = (ix << 16) | (i0 >> 16);
 312
 313   /* purge off +-inf, NaN, +-0, and negative arguments */
 314   if (ix >= 0x7fff0000)
 315     return x * x;
 316
 317   if (ix < 0x3fc08000) /* 2^-63 */
 318     {                           /* |x|<2**-63, return -log(|x|) */
 319       if (se & 0x8000)
 320         {
 321           *signgamp = -1;
 322           return -__ieee754_logl (-x);
 323         }
 324       else
 325         return -__ieee754_logl (x);
 326     }
 327   if (se & 0x8000)
 328     {
 329       t = sin_pi (x);
 330       if (t == zero)
 331         return one / fabsl (t); /* -integer */
 332       nadj = __ieee754_logl (pi / fabsl (t * x));
 333       if (t < zero)
 334         *signgamp = -1;
 335       x = -x;
 336     }
 337
 338   /* purge off 1 and 2 */
 339   if ((((ix - 0x3fff8000) | i0 | i1) == 0)
 340       || (((ix - 0x40008000) | i0 | i1) == 0))
 341     r = 0;
 342   else if (ix < 0x40008000) /* 2.0 */
 343     {
 344       /* x < 2.0 */
 345       if (ix <= 0x3ffee666) /* 8.99993896484375e-1 */
 346         {
 347           /* lgamma(x) = lgamma(x+1) - log(x) */
 348           r = -__ieee754_logl (x);
 349           if (ix >= 0x3ffebb4a) /* 7.31597900390625e-1 */
 350             {
 351               y = x - one;
 352               i = 0;
 353             }
 354           else if (ix >= 0x3ffced33)/* 2.31639862060546875e-1 */
 355             {
 356               y = x - (tc - one);
 357               i = 1;
 358             }
 359           else
 360             {
 361               /* x < 0.23 */
 362               y = x;
 363               i = 2;
 364             }
 365         }
 366       else
 367         {
 368           r = zero;
 369           if (ix >= 0x3fffdda6) /* 1.73162841796875 */
 370             {
 371               /* [1.7316,2] */
 372               y = x - 2.0;
 373               i = 0;
 374             }
 375           else if (ix >= 0x3fff9da6)/* 1.23162841796875 */
 376             {
 377               /* [1.23,1.73] */
 378               y = x - tc;
 379               i = 1;
 380             }
 381           else
 382             {
 383               /* [0.9, 1.23] */
 384               y = x - one;
 385               i = 2;
 386             }
 387         }
 388       switch (i)
 389         {
 390         case 0:
 391           p1 = a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * a5))));
 392           p2 = b0 + y * (b1 + y * (b2 + y * (b3 + y * (b4 + y))));
 393           r += half * y + y * p1/p2;
 394           break;
 395         case 1:
 396     p1 = g0 + y * (g1 + y * (g2 + y * (g3 + y * (g4 + y * (g5 + y * g6)))));
 397     p2 = h0 + y * (h1 + y * (h2 + y * (h3 + y * (h4 + y * (h5 + y)))));
 398     p = tt + y * p1/p2;
 399           r += (tf + p);
 400           break;
 401         case 2:
 402  p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * (u5 + y * u6))))));
 403       p2 = v0 + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * (v5 + y)))));
 404           r += (-half * y + p1 / p2);
 405         }
 406     }
 407   else if (ix < 0x40028000) /* 8.0 */
 408     {
 409       /* x < 8.0 */
 410       i = (int) x;
 411       t = zero;
 412       y = x - (double) i;
 413   p = y *
 414      (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
 415   q = r0 + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * (r6 + y))))));
 416       r = half * y + p / q;
 417       z = one;                  /* lgamma(1+s) = log(s) + lgamma(s) */
 418       switch (i)
 419         {
 420         case 7:
 421           z *= (y + 6.0);       /* FALLTHRU */
 422         case 6:
 423           z *= (y + 5.0);       /* FALLTHRU */
 424         case 5:
 425           z *= (y + 4.0);       /* FALLTHRU */
 426         case 4:
 427           z *= (y + 3.0);       /* FALLTHRU */
 428         case 3:
 429           z *= (y + 2.0);       /* FALLTHRU */
 430           r += __ieee754_logl (z);
 431           break;
 432         }
 433     }
 434   else if (ix < 0x40418000) /* 2^66 */
 435     {
 436       /* 8.0 <= x < 2**66 */
 437       t = __ieee754_logl (x);
 438       z = one / x;
 439       y = z * z;
 440       w = w0 + z * (w1
 441           + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * (w6 + y * w7))))));
 442       r = (x - half) * (t - one) + w;
 443     }
 444   else
 445     /* 2**66 <= x <= inf */
 446     r = x * (__ieee754_logl (x) - one);
 447   if (se & 0x8000)
 448     r = nadj - r;
 449   return r;
 450 }