Back to index

glibc  2.9
e_lgammal_r.c
Go to the documentation of this file.
00001 /*
00002  * ====================================================
00003  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00004  *
00005  * Developed at SunPro, a Sun Microsystems, Inc. business.
00006  * Permission to use, copy, modify, and distribute this
00007  * software is freely granted, provided that this notice
00008  * is preserved.
00009  * ====================================================
00010  */
00011 
00012 /* Long double expansions are
00013   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
00014   and are incorporated herein by permission of the author.  The author
00015   reserves the right to distribute this material elsewhere under different
00016   copying permissions.  These modifications are distributed here under
00017   the following terms:
00018 
00019     This library is free software; you can redistribute it and/or
00020     modify it under the terms of the GNU Lesser General Public
00021     License as published by the Free Software Foundation; either
00022     version 2.1 of the License, or (at your option) any later version.
00023 
00024     This library is distributed in the hope that it will be useful,
00025     but WITHOUT ANY WARRANTY; without even the implied warranty of
00026     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00027     Lesser General Public License for more details.
00028 
00029     You should have received a copy of the GNU Lesser General Public
00030     License along with this library; if not, write to the Free Software
00031     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00032 
00033 /* __ieee754_lgammal_r(x, signgamp)
00034  * Reentrant version of the logarithm of the Gamma function
00035  * with user provide pointer for the sign of Gamma(x).
00036  *
00037  * Method:
00038  *   1. Argument Reduction for 0 < x <= 8
00039  *     Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
00040  *     reduce x to a number in [1.5,2.5] by
00041  *            lgamma(1+s) = log(s) + lgamma(s)
00042  *     for example,
00043  *            lgamma(7.3) = log(6.3) + lgamma(6.3)
00044  *                       = log(6.3*5.3) + lgamma(5.3)
00045  *                       = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
00046  *   2. Polynomial approximation of lgamma around its
00047  *     minimun ymin=1.461632144968362245 to maintain monotonicity.
00048  *     On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
00049  *            Let z = x-ymin;
00050  *            lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
00051  *   2. Rational approximation in the primary interval [2,3]
00052  *     We use the following approximation:
00053  *            s = x-2.0;
00054  *            lgamma(x) = 0.5*s + s*P(s)/Q(s)
00055  *     Our algorithms are based on the following observation
00056  *
00057  *                             zeta(2)-1    2    zeta(3)-1    3
00058  * lgamma(2+s) = s*(1-Euler) + --------- * s  -  --------- * s  + ...
00059  *                                 2                 3
00060  *
00061  *     where Euler = 0.5771... is the Euler constant, which is very
00062  *     close to 0.5.
00063  *
00064  *   3. For x>=8, we have
00065  *     lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
00066  *     (better formula:
00067  *        lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
00068  *     Let z = 1/x, then we approximation
00069  *            f(z) = lgamma(x) - (x-0.5)(log(x)-1)
00070  *     by
00071  *                              3       5             11
00072  *            w = w0 + w1*z + w2*z  + w3*z  + ... + w6*z
00073  *
00074  *   4. For negative x, since (G is gamma function)
00075  *            -x*G(-x)*G(x) = pi/sin(pi*x),
00076  *     we have
00077  *            G(x) = pi/(sin(pi*x)*(-x)*G(-x))
00078  *     since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
00079  *     Hence, for x<0, signgam = sign(sin(pi*x)) and
00080  *            lgamma(x) = log(|Gamma(x)|)
00081  *                     = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
00082  *     Note: one should avoid compute pi*(-x) directly in the
00083  *           computation of sin(pi*(-x)).
00084  *
00085  *   5. Special Cases
00086  *            lgamma(2+s) ~ s*(1-Euler) for tiny s
00087  *            lgamma(1)=lgamma(2)=0
00088  *            lgamma(x) ~ -log(x) for tiny x
00089  *            lgamma(0) = lgamma(inf) = inf
00090  *            lgamma(-integer) = +-inf
00091  *
00092  */
00093 
00094 #include "math.h"
00095 #include "math_private.h"
00096 
00097 #ifdef __STDC__
00098 static const long double
00099 #else
00100 static long double
00101 #endif
00102   half = 0.5L,
00103   one = 1.0L,
00104   pi = 3.14159265358979323846264L,
00105   two63 = 9.223372036854775808e18L,
00106 
00107   /* lgam(1+x) = 0.5 x + x a(x)/b(x)
00108      -0.268402099609375 <= x <= 0
00109      peak relative error 6.6e-22 */
00110   a0 = -6.343246574721079391729402781192128239938E2L,
00111   a1 =  1.856560238672465796768677717168371401378E3L,
00112   a2 =  2.404733102163746263689288466865843408429E3L,
00113   a3 =  8.804188795790383497379532868917517596322E2L,
00114   a4 =  1.135361354097447729740103745999661157426E2L,
00115   a5 =  3.766956539107615557608581581190400021285E0L,
00116 
00117   b0 =  8.214973713960928795704317259806842490498E3L,
00118   b1 =  1.026343508841367384879065363925870888012E4L,
00119   b2 =  4.553337477045763320522762343132210919277E3L,
00120   b3 =  8.506975785032585797446253359230031874803E2L,
00121   b4 =  6.042447899703295436820744186992189445813E1L,
00122   /* b5 =  1.000000000000000000000000000000000000000E0 */
00123 
00124 
00125   tc =  1.4616321449683623412626595423257213284682E0L,
00126   tf = -1.2148629053584961146050602565082954242826E-1,/* double precision */
00127 /* tt = (tail of tf), i.e. tf + tt has extended precision. */
00128   tt = 3.3649914684731379602768989080467587736363E-18L,
00129   /* lgam ( 1.4616321449683623412626595423257213284682E0 ) =
00130 -1.2148629053584960809551455717769158215135617312999903886372437313313530E-1 */
00131 
00132   /* lgam (x + tc) = tf + tt + x g(x)/h(x)
00133      - 0.230003726999612341262659542325721328468 <= x
00134      <= 0.2699962730003876587373404576742786715318
00135      peak relative error 2.1e-21 */
00136   g0 = 3.645529916721223331888305293534095553827E-18L,
00137   g1 = 5.126654642791082497002594216163574795690E3L,
00138   g2 = 8.828603575854624811911631336122070070327E3L,
00139   g3 = 5.464186426932117031234820886525701595203E3L,
00140   g4 = 1.455427403530884193180776558102868592293E3L,
00141   g5 = 1.541735456969245924860307497029155838446E2L,
00142   g6 = 4.335498275274822298341872707453445815118E0L,
00143 
00144   h0 = 1.059584930106085509696730443974495979641E4L,
00145   h1 =  2.147921653490043010629481226937850618860E4L,
00146   h2 = 1.643014770044524804175197151958100656728E4L,
00147   h3 =  5.869021995186925517228323497501767586078E3L,
00148   h4 =  9.764244777714344488787381271643502742293E2L,
00149   h5 =  6.442485441570592541741092969581997002349E1L,
00150   /* h6 = 1.000000000000000000000000000000000000000E0 */
00151 
00152 
00153   /* lgam (x+1) = -0.5 x + x u(x)/v(x)
00154      -0.100006103515625 <= x <= 0.231639862060546875
00155      peak relative error 1.3e-21 */
00156   u0 = -8.886217500092090678492242071879342025627E1L,
00157   u1 =  6.840109978129177639438792958320783599310E2L,
00158   u2 =  2.042626104514127267855588786511809932433E3L,
00159   u3 =  1.911723903442667422201651063009856064275E3L,
00160   u4 =  7.447065275665887457628865263491667767695E2L,
00161   u5 =  1.132256494121790736268471016493103952637E2L,
00162   u6 =  4.484398885516614191003094714505960972894E0L,
00163 
00164   v0 =  1.150830924194461522996462401210374632929E3L,
00165   v1 =  3.399692260848747447377972081399737098610E3L,
00166   v2 =  3.786631705644460255229513563657226008015E3L,
00167   v3 =  1.966450123004478374557778781564114347876E3L,
00168   v4 =  4.741359068914069299837355438370682773122E2L,
00169   v5 =  4.508989649747184050907206782117647852364E1L,
00170   /* v6 =  1.000000000000000000000000000000000000000E0 */
00171 
00172 
00173   /* lgam (x+2) = .5 x + x s(x)/r(x)
00174      0 <= x <= 1
00175      peak relative error 7.2e-22 */
00176   s0 =  1.454726263410661942989109455292824853344E6L,
00177   s1 = -3.901428390086348447890408306153378922752E6L,
00178   s2 = -6.573568698209374121847873064292963089438E6L,
00179   s3 = -3.319055881485044417245964508099095984643E6L,
00180   s4 = -7.094891568758439227560184618114707107977E5L,
00181   s5 = -6.263426646464505837422314539808112478303E4L,
00182   s6 = -1.684926520999477529949915657519454051529E3L,
00183 
00184   r0 = -1.883978160734303518163008696712983134698E7L,
00185   r1 = -2.815206082812062064902202753264922306830E7L,
00186   r2 = -1.600245495251915899081846093343626358398E7L,
00187   r3 = -4.310526301881305003489257052083370058799E6L,
00188   r4 = -5.563807682263923279438235987186184968542E5L,
00189   r5 = -3.027734654434169996032905158145259713083E4L,
00190   r6 = -4.501995652861105629217250715790764371267E2L,
00191   /* r6 =  1.000000000000000000000000000000000000000E0 */
00192 
00193 
00194 /* lgam(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x w(1/x^2)
00195    x >= 8
00196    Peak relative error 1.51e-21
00197    w0 = LS2PI - 0.5 */
00198   w0 =  4.189385332046727417803e-1L,
00199   w1 =  8.333333333333331447505E-2L,
00200   w2 = -2.777777777750349603440E-3L,
00201   w3 =  7.936507795855070755671E-4L,
00202   w4 = -5.952345851765688514613E-4L,
00203   w5 =  8.412723297322498080632E-4L,
00204   w6 = -1.880801938119376907179E-3L,
00205   w7 =  4.885026142432270781165E-3L;
00206 
00207 #ifdef __STDC__
00208 static const long double zero = 0.0L;
00209 #else
00210 static long double zero = 0.0L;
00211 #endif
00212 
00213 #ifdef __STDC__
00214 static long double
00215 sin_pi (long double x)
00216 #else
00217 static long double
00218 sin_pi (x)
00219      long double x;
00220 #endif
00221 {
00222   long double y, z;
00223   int n, ix;
00224   u_int32_t se, i0, i1;
00225 
00226   GET_LDOUBLE_WORDS (se, i0, i1, x);
00227   ix = se & 0x7fff;
00228   ix = (ix << 16) | (i0 >> 16);
00229   if (ix < 0x3ffd8000) /* 0.25 */
00230     return __sinl (pi * x);
00231   y = -x;                   /* x is assume negative */
00232 
00233   /*
00234    * argument reduction, make sure inexact flag not raised if input
00235    * is an integer
00236    */
00237   z = __floorl (y);
00238   if (z != y)
00239     {                       /* inexact anyway */
00240       y  *= 0.5;
00241       y = 2.0*(y - __floorl(y));          /* y = |x| mod 2.0 */
00242       n = (int) (y*4.0);
00243     }
00244   else
00245     {
00246       if (ix >= 0x403f8000)  /* 2^64 */
00247        {
00248          y = zero; n = 0;                 /* y must be even */
00249        }
00250       else
00251        {
00252        if (ix < 0x403e8000)  /* 2^63 */
00253          z = y + two63;     /* exact */
00254        GET_LDOUBLE_WORDS (se, i0, i1, z);
00255        n = i1 & 1;
00256        y  = n;
00257        n <<= 2;
00258       }
00259     }
00260 
00261   switch (n)
00262     {
00263     case 0:
00264       y = __sinl (pi * y);
00265       break;
00266     case 1:
00267     case 2:
00268       y = __cosl (pi * (half - y));
00269       break;
00270     case 3:
00271     case 4:
00272       y = __sinl (pi * (one - y));
00273       break;
00274     case 5:
00275     case 6:
00276       y = -__cosl (pi * (y - 1.5));
00277       break;
00278     default:
00279       y = __sinl (pi * (y - 2.0));
00280       break;
00281     }
00282   return -y;
00283 }
00284 
00285 
00286 #ifdef __STDC__
00287 long double
00288 __ieee754_lgammal_r (long double x, int *signgamp)
00289 #else
00290 long double
00291 __ieee754_lgammal_r (x, signgamp)
00292      long double x;
00293      int *signgamp;
00294 #endif
00295 {
00296   long double t, y, z, nadj, p, p1, p2, q, r, w;
00297   int i, ix;
00298   u_int32_t se, i0, i1;
00299 
00300   *signgamp = 1;
00301   GET_LDOUBLE_WORDS (se, i0, i1, x);
00302   ix = se & 0x7fff;
00303 
00304   if ((ix | i0 | i1) == 0)
00305     {
00306       if (se & 0x8000)
00307        *signgamp = -1;
00308       return one / fabsl (x);
00309     }
00310 
00311   ix = (ix << 16) | (i0 >> 16);
00312 
00313   /* purge off +-inf, NaN, +-0, and negative arguments */
00314   if (ix >= 0x7fff0000)
00315     return x * x;
00316 
00317   if (ix < 0x3fc08000) /* 2^-63 */
00318     {                       /* |x|<2**-63, return -log(|x|) */
00319       if (se & 0x8000)
00320        {
00321          *signgamp = -1;
00322          return -__ieee754_logl (-x);
00323        }
00324       else
00325        return -__ieee754_logl (x);
00326     }
00327   if (se & 0x8000)
00328     {
00329       t = sin_pi (x);
00330       if (t == zero)
00331        return one / fabsl (t);     /* -integer */
00332       nadj = __ieee754_logl (pi / fabsl (t * x));
00333       if (t < zero)
00334        *signgamp = -1;
00335       x = -x;
00336     }
00337 
00338   /* purge off 1 and 2 */
00339   if ((((ix - 0x3fff8000) | i0 | i1) == 0)
00340       || (((ix - 0x40008000) | i0 | i1) == 0))
00341     r = 0;
00342   else if (ix < 0x40008000) /* 2.0 */
00343     {
00344       /* x < 2.0 */
00345       if (ix <= 0x3ffee666) /* 8.99993896484375e-1 */
00346        {
00347          /* lgamma(x) = lgamma(x+1) - log(x) */
00348          r = -__ieee754_logl (x);
00349          if (ix >= 0x3ffebb4a) /* 7.31597900390625e-1 */
00350            {
00351              y = x - one;
00352              i = 0;
00353            }
00354          else if (ix >= 0x3ffced33)/* 2.31639862060546875e-1 */
00355            {
00356              y = x - (tc - one);
00357              i = 1;
00358            }
00359          else
00360            {
00361              /* x < 0.23 */
00362              y = x;
00363              i = 2;
00364            }
00365        }
00366       else
00367        {
00368          r = zero;
00369          if (ix >= 0x3fffdda6) /* 1.73162841796875 */
00370            {
00371              /* [1.7316,2] */
00372              y = x - 2.0;
00373              i = 0;
00374            }
00375          else if (ix >= 0x3fff9da6)/* 1.23162841796875 */
00376            {
00377              /* [1.23,1.73] */
00378              y = x - tc;
00379              i = 1;
00380            }
00381          else
00382            {
00383              /* [0.9, 1.23] */
00384              y = x - one;
00385              i = 2;
00386            }
00387        }
00388       switch (i)
00389        {
00390        case 0:
00391          p1 = a0 + y * (a1 + y * (a2 + y * (a3 + y * (a4 + y * a5))));
00392          p2 = b0 + y * (b1 + y * (b2 + y * (b3 + y * (b4 + y))));
00393          r += half * y + y * p1/p2;
00394          break;
00395        case 1:
00396     p1 = g0 + y * (g1 + y * (g2 + y * (g3 + y * (g4 + y * (g5 + y * g6)))));
00397     p2 = h0 + y * (h1 + y * (h2 + y * (h3 + y * (h4 + y * (h5 + y)))));
00398     p = tt + y * p1/p2;
00399          r += (tf + p);
00400          break;
00401        case 2:
00402  p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * (u5 + y * u6))))));
00403       p2 = v0 + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * (v5 + y)))));
00404          r += (-half * y + p1 / p2);
00405        }
00406     }
00407   else if (ix < 0x40028000) /* 8.0 */
00408     {
00409       /* x < 8.0 */
00410       i = (int) x;
00411       t = zero;
00412       y = x - (double) i;
00413   p = y *
00414      (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
00415   q = r0 + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * (r6 + y))))));
00416       r = half * y + p / q;
00417       z = one;                     /* lgamma(1+s) = log(s) + lgamma(s) */
00418       switch (i)
00419        {
00420        case 7:
00421          z *= (y + 6.0);    /* FALLTHRU */
00422        case 6:
00423          z *= (y + 5.0);    /* FALLTHRU */
00424        case 5:
00425          z *= (y + 4.0);    /* FALLTHRU */
00426        case 4:
00427          z *= (y + 3.0);    /* FALLTHRU */
00428        case 3:
00429          z *= (y + 2.0);    /* FALLTHRU */
00430          r += __ieee754_logl (z);
00431          break;
00432        }
00433     }
00434   else if (ix < 0x40418000) /* 2^66 */
00435     {
00436       /* 8.0 <= x < 2**66 */
00437       t = __ieee754_logl (x);
00438       z = one / x;
00439       y = z * z;
00440       w = w0 + z * (w1
00441           + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * (w6 + y * w7))))));
00442       r = (x - half) * (t - one) + w;
00443     }
00444   else
00445     /* 2**66 <= x <= inf */
00446     r = x * (__ieee754_logl (x) - one);
00447   if (se & 0x8000)
00448     r = nadj - r;
00449   return r;
00450 }