Print this page
11210 libm should be cstyle(1ONBLD) clean


   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */

  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29 
  30 /*
  31  * long double __k_lgammal(long double x, int *signgamlp);
  32  * K.C. Ng, August, 1989.
  33  *
  34  * We choose [1.5,2.5] to be the primary interval. Our algorithms
  35  * are mainly derived from
  36  *
  37  *
  38  *                             zeta(2)-1    2    zeta(3)-1    3
  39  * lgamma(2+s) = s*(1-euler) + --------- * s  -  --------- * s  + ...
  40  *                                 2                 3
  41  *
  42  *
  43  * Note 1. Since gamma(1+s)=s*gamma(s), hence
  44  *              lgamma(1+s) = log(s) + lgamma(s), or
  45  *              lgamma(s) = lgamma(1+s) - log(s).
  46  *         When s is really tiny (like roundoff), lgamma(1+s) ~ s(1-enler)
  47  *         Hence lgamma(s) ~ -log(s) for tiny s
  48  *
  49  */
  50 
  51 #include "libm.h"
  52 #include "longdouble.h"
  53 
  54 static long double neg(long double, int *);
  55 static long double poly(long double, const long double *, int);
  56 static long double polytail(long double);
  57 static long double primary(long double);
  58 
  59 static const long double
  60 c0 =     0.0L,
  61 ch =     0.5L,
  62 c1 =     1.0L,
  63 c2 =     2.0L,
  64 c3 =     3.0L,
  65 c4 =     4.0L,
  66 c5 =     5.0L,
  67 c6 =     6.0L,
  68 pi =     3.1415926535897932384626433832795028841971L,
  69 tiny =   1.0e-40L;
  70 
  71 long double
  72 __k_lgammal(long double x, int *signgamlp) {

  73         long double t, y;
  74         int i;
  75 
  76         /* purge off +-inf, NaN and negative arguments */
  77         if (!finitel(x))
  78                 return (x*x);

  79         *signgamlp = 1;

  80         if (signbitl(x))
  81                 return (neg(x, signgamlp));
  82 
  83         /* for x < 8.0 */
  84         if (x < 8.0L) {
  85             y = anintl(x);
  86             i = (int) y;

  87             switch (i) {
  88             case 0:

  89                 if (x < 1.0e-40L)
  90                         return (-logl(x));
  91                 else
  92                         return (primary(x)-log1pl(x))-logl(x);

  93             case 1:
  94                 return (primary(x-y)-logl(x));
  95             case 2:
  96                 return (primary(x-y));
  97             case 3:
  98                 return (primary(x-y)+logl(x-c1));
  99             case 4:
 100                 return (primary(x-y)+logl((x-c1)*(x-c2)));
 101             case 5:
 102                 return (primary(x-y)+logl((x-c1)*(x-c2)*(x-c3)));

 103             case 6:
 104                 return (primary(x-y)+logl((x-c1)*(x-c2)*(x-c3)*(x-c4)));

 105             case 7:
 106                 return (primary(x-y)+logl((x-c1)*(x-c2)*(x-c3)*(x-c4)*(x-c5)));

 107             case 8:
 108                 return primary(x-y)+
 109                         logl((x-c1)*(x-c2)*(x-c3)*(x-c4)*(x-c5)*(x-c6));
 110             }
 111         }
 112 
 113         /* 8.0 <= x < 1.0e40 */
 114         if (x < 1.0e40L) {
 115             t = logl(x);
 116             return (x*(t-c1)-(ch*t-polytail(c1/x)));
 117         }
 118 
 119         /* 1.0e40 <= x <= inf */
 120         return (x*(logl(x)-c1));
 121 }
 122 
 123 static const long double an1[] = {              /* 20 terms */
 124         -0.0772156649015328606065120900824024309741L,
 125         3.224670334241132182362075833230130289059e-0001L,
 126         -6.735230105319809513324605383668929964120e-0002L,
 127         2.058080842778454787900092432928910226297e-0002L,
 128         -7.385551028673985266273054086081102125704e-0003L,
 129         2.890510330741523285758867304409628648727e-0003L,
 130         -1.192753911703260976581414338096267498555e-0003L,
 131         5.096695247430424562831956662855697824035e-0004L,
 132         -2.231547584535777978926798502084300123638e-0004L,
 133         9.945751278186384670278268034322157947635e-0005L,
 134         -4.492623673665547726647838474125147631082e-0005L,
 135         2.050721280617796810096993154281561168706e-0005L,
 136         -9.439487785617396552092393234044767313568e-0006L,
 137         4.374872903516051510689234173139793159340e-0006L,
 138         -2.039156676413643091040459825776029327487e-0006L,
 139         9.555777181318621470466563543806211523634e-0007L,
 140         -4.468344919709630637558538313482398989638e-0007L,


 285   -6.735230105302832007479431772160948499254e-0002L,
 286    2.058080842553481183648529360967441889912e-0002L,
 287   -7.385551007602909242024706804659879199244e-0003L,
 288    2.890510182473907253939821312248303471206e-0003L,
 289   -1.192753098427856770847894497586825614450e-0003L,
 290    5.096659636418811568063339214203693550804e-0004L,
 291   -2.231421144004355691166194259675004483639e-0004L,
 292    9.942073842343832132754332881883387625136e-0005L,
 293   -4.483809261973204531263252655050701205397e-0005L,
 294    2.033260142610284888319116654931994447173e-0005L,
 295   -9.153539544026646699870528191410440585796e-0006L,
 296    3.988460469925482725894144688699584997971e-0006L,
 297   -1.609692980087029172567957221850825977621e-0006L,
 298    5.634916377249975825399706694496688803488e-0007L,
 299   -1.560065465929518563549083208482591437696e-0007L,
 300    2.961350193868935325526962209019387821584e-0008L,
 301   -2.834602215195368130104649234505033159842e-0009L,
 302 };
 303 
 304 static long double
 305 primary(long double s) {        /* assume |s|<=0.5 */

 306         int i;
 307 
 308         i = (int) (8.0L * (s + 0.5L));

 309         switch (i) {
 310         case 0: return ch*s+s*poly(s, an4, 21);
 311         case 1: return ch*s+s*poly(s, an3, 20);
 312         case 2: return ch*s+s*poly(s, an2, 20);
 313         case 3: return ch*s+s*poly(s, an1, 20);
 314         case 4: return ch*s+s*poly(s, ap1, 19);
 315         case 5: return ch*s+s*poly(s, ap2, 19);
 316         case 6: return ch*s+s*poly(s, ap3, 19);
 317         case 7: return ch*s+s*poly(s, ap4, 19);








 318         }

 319         /* NOTREACHED */
 320     return (0.0L);
 321 }
 322 
 323 static long double
 324 poly(long double s, const long double *p, int n) {

 325         long double y;
 326         int i;
 327         y = p[n-1];
 328         for (i = n-2; i >= 0; i--) y = p[i]+s*y;




 329         return (y);
 330 }
 331 
 332 static const long double pt[] = {
 333    9.189385332046727417803297364056176804663e-0001L,
 334    8.333333333333333333333333333331286969123e-0002L,
 335   -2.777777777777777777777777553194796036402e-0003L,
 336    7.936507936507936507927283071433584248176e-0004L,
 337   -5.952380952380952362351042163192634108297e-0004L,
 338    8.417508417508395661774286645578379460131e-0004L,
 339   -1.917526917525263651186066417934685675649e-0003L,
 340    6.410256409395203164659292973142293199083e-0003L,
 341   -2.955065327248303301763594514012418438188e-0002L,
 342    1.796442830099067542945998615411893822886e-0001L,
 343   -1.392413465829723742489974310411118662919e+0000L,
 344    1.339984238037267658352656597960492029261e+0001L,
 345   -1.564707657605373662425785904278645727813e+0002L,
 346    2.156323807499211356127813962223067079300e+0003L,
 347   -3.330486427626223184647299834137041307569e+0004L,
 348    5.235535072011889213611369254140123518699e+0005L,
 349   -7.258160984602220710491988573430212593080e+0006L,
 350    7.316526934569686459641438882340322673357e+0007L,
 351   -3.806450279064900548836571789284896711473e+0008L,
 352 };
 353 
 354 static long double
 355 polytail(long double s) {

 356         long double t, z;
 357         int i;
 358         z = s*s;

 359         t = pt[18];
 360         for (i = 17; i >= 1; i--) t = pt[i]+z*t;
 361         return (pt[0]+s*t);



 362 }
 363 
 364 static long double
 365 neg(long double z, int *signgamlp) {

 366         long double t, p;
 367 

 368         /*
 369          * written by K.C. Ng,  Feb 2, 1989.
 370          *
 371          * Since
 372          *              -z*G(-z)*G(z) = pi/sin(pi*z),
 373          * we have
 374          *      G(-z) = -pi/(sin(pi*z)*G(z)*z)
 375          *                =  pi/(sin(pi*(-z))*G(z)*z)
 376          * Algorithm
 377          *              z = |z|
 378          *              t = sinpi(z); ...note that when z>2**112, z is an int
 379          *              and hence t=0.
 380          *
 381          *              if (t == 0.0) return 1.0/0.0;
 382          *              if (t< 0.0) *signgamlp = -1; else t= -t;
 383          *              if (z<1.0e-40)       ...tiny z
 384          *                  return -log(z);
 385          *              else
 386          *                  return log(pi/(t*z))-lgamma(z);
 387          *
 388          */

 389 
 390         t = sinpil(z);                  /* t := sin(pi*z) */

 391         if (t == c0)                    /* return   1.0/0.0 =  +INF */
 392             return (c1/c0);
 393 
 394         z = -z;

 395         if (z <= tiny)
 396             p = -logl(z);
 397         else
 398                 p = logl(pi/(fabsl(t)*z)) - __k_lgammal(z, signgamlp);
 399         if (t < c0) *signgamlp = -1;



 400         return (p);
 401 }


   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 
  26 /*
  27  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 /*
  32  * long double __k_lgammal(long double x, int *signgamlp);
  33  * K.C. Ng, August, 1989.
  34  *
  35  * We choose [1.5,2.5] to be the primary interval. Our algorithms
  36  * are mainly derived from
  37  *
  38  *
  39  *                             zeta(2)-1    2    zeta(3)-1    3
  40  * lgamma(2+s) = s*(1-euler) + --------- * s  -  --------- * s  + ...
  41  *                                 2                 3
  42  *
  43  *
  44  * Note 1. Since gamma(1+s)=s*gamma(s), hence
  45  *              lgamma(1+s) = log(s) + lgamma(s), or
  46  *              lgamma(s) = lgamma(1+s) - log(s).
  47  *         When s is really tiny (like roundoff), lgamma(1+s) ~ s(1-enler)
  48  *         Hence lgamma(s) ~ -log(s) for tiny s
  49  *
  50  */
  51 
  52 #include "libm.h"
  53 #include "longdouble.h"
  54 
  55 static long double neg(long double, int *);
  56 static long double poly(long double, const long double *, int);
  57 static long double polytail(long double);
  58 static long double primary(long double);
  59 
  60 static const long double c0 = 0.0L,
  61         ch = 0.5L,
  62         c1 = 1.0L,
  63         c2 = 2.0L,
  64         c3 = 3.0L,
  65         c4 = 4.0L,
  66         c5 = 5.0L,
  67         c6 = 6.0L,
  68         pi = 3.1415926535897932384626433832795028841971L,
  69         tiny = 1.0e-40L;

  70 
  71 long double
  72 __k_lgammal(long double x, int *signgamlp)
  73 {
  74         long double t, y;
  75         int i;
  76 
  77         /* purge off +-inf, NaN and negative arguments */
  78         if (!finitel(x))
  79                 return (x * x);
  80 
  81         *signgamlp = 1;
  82 
  83         if (signbitl(x))
  84                 return (neg(x, signgamlp));
  85 
  86         /* for x < 8.0 */
  87         if (x < 8.0L) {
  88                 y = anintl(x);
  89                 i = (int)y;
  90 
  91                 switch (i) {
  92                 case 0:
  93 
  94                         if (x < 1.0e-40L)
  95                                 return (-logl(x));
  96                         else
  97                                 return ((primary(x) - log1pl(x)) - logl(x));
  98 
  99                 case 1:
 100                         return (primary(x - y) - logl(x));
 101                 case 2:
 102                         return (primary(x - y));
 103                 case 3:
 104                         return (primary(x - y) + logl(x - c1));
 105                 case 4:
 106                         return (primary(x - y) + logl((x - c1) * (x - c2)));
 107                 case 5:
 108                         return (primary(x - y) + logl((x - c1) * (x - c2) * (x -
 109                             c3)));
 110                 case 6:
 111                         return (primary(x - y) + logl((x - c1) * (x - c2) * (x -
 112                             c3) * (x - c4)));
 113                 case 7:
 114                         return (primary(x - y) + logl((x - c1) * (x - c2) * (x -
 115                             c3) * (x - c4) * (x - c5)));
 116                 case 8:
 117                         return (primary(x - y) + logl((x - c1) * (x - c2) * (x -
 118                             c3) * (x - c4) * (x - c5) * (x - c6)));
 119                 }
 120         }
 121 
 122         /* 8.0 <= x < 1.0e40 */
 123         if (x < 1.0e40L) {
 124                 t = logl(x);
 125                 return (x * (t - c1) - (ch * t - polytail(c1 / x)));
 126         }
 127 
 128         /* 1.0e40 <= x <= inf */
 129         return (x * (logl(x) - c1));
 130 }
 131 
 132 static const long double an1[] = {      /* 20 terms */
 133         -0.0772156649015328606065120900824024309741L,
 134         3.224670334241132182362075833230130289059e-0001L,
 135         -6.735230105319809513324605383668929964120e-0002L,
 136         2.058080842778454787900092432928910226297e-0002L,
 137         -7.385551028673985266273054086081102125704e-0003L,
 138         2.890510330741523285758867304409628648727e-0003L,
 139         -1.192753911703260976581414338096267498555e-0003L,
 140         5.096695247430424562831956662855697824035e-0004L,
 141         -2.231547584535777978926798502084300123638e-0004L,
 142         9.945751278186384670278268034322157947635e-0005L,
 143         -4.492623673665547726647838474125147631082e-0005L,
 144         2.050721280617796810096993154281561168706e-0005L,
 145         -9.439487785617396552092393234044767313568e-0006L,
 146         4.374872903516051510689234173139793159340e-0006L,
 147         -2.039156676413643091040459825776029327487e-0006L,
 148         9.555777181318621470466563543806211523634e-0007L,
 149         -4.468344919709630637558538313482398989638e-0007L,


 294         -6.735230105302832007479431772160948499254e-0002L,
 295         2.058080842553481183648529360967441889912e-0002L,
 296         -7.385551007602909242024706804659879199244e-0003L,
 297         2.890510182473907253939821312248303471206e-0003L,
 298         -1.192753098427856770847894497586825614450e-0003L,
 299         5.096659636418811568063339214203693550804e-0004L,
 300         -2.231421144004355691166194259675004483639e-0004L,
 301         9.942073842343832132754332881883387625136e-0005L,
 302         -4.483809261973204531263252655050701205397e-0005L,
 303         2.033260142610284888319116654931994447173e-0005L,
 304         -9.153539544026646699870528191410440585796e-0006L,
 305         3.988460469925482725894144688699584997971e-0006L,
 306         -1.609692980087029172567957221850825977621e-0006L,
 307         5.634916377249975825399706694496688803488e-0007L,
 308         -1.560065465929518563549083208482591437696e-0007L,
 309         2.961350193868935325526962209019387821584e-0008L,
 310         -2.834602215195368130104649234505033159842e-0009L,
 311 };
 312 
 313 static long double
 314 primary(long double s)                  /* assume |s|<=0.5 */
 315 {
 316         int i;
 317 
 318         i = (int)(8.0L * (s + 0.5L));
 319 
 320         switch (i) {
 321         case 0:
 322                 return (ch * s + s * poly(s, an4, 21));
 323         case 1:
 324                 return (ch * s + s * poly(s, an3, 20));
 325         case 2:
 326                 return (ch * s + s * poly(s, an2, 20));
 327         case 3:
 328                 return (ch * s + s * poly(s, an1, 20));
 329         case 4:
 330                 return (ch * s + s * poly(s, ap1, 19));
 331         case 5:
 332                 return (ch * s + s * poly(s, ap2, 19));
 333         case 6:
 334                 return (ch * s + s * poly(s, ap3, 19));
 335         case 7:
 336                 return (ch * s + s * poly(s, ap4, 19));
 337         }
 338 
 339         /* NOTREACHED */
 340         return (0.0L);
 341 }
 342 
 343 static long double
 344 poly(long double s, const long double *p, int n)
 345 {
 346         long double y;
 347         int i;
 348 
 349         y = p[n - 1];
 350 
 351         for (i = n - 2; i >= 0; i--)
 352                 y = p[i] + s * y;
 353 
 354         return (y);
 355 }
 356 
 357 static const long double pt[] = {
 358         9.189385332046727417803297364056176804663e-0001L,
 359         8.333333333333333333333333333331286969123e-0002L,
 360         -2.777777777777777777777777553194796036402e-0003L,
 361         7.936507936507936507927283071433584248176e-0004L,
 362         -5.952380952380952362351042163192634108297e-0004L,
 363         8.417508417508395661774286645578379460131e-0004L,
 364         -1.917526917525263651186066417934685675649e-0003L,
 365         6.410256409395203164659292973142293199083e-0003L,
 366         -2.955065327248303301763594514012418438188e-0002L,
 367         1.796442830099067542945998615411893822886e-0001L,
 368         -1.392413465829723742489974310411118662919e+0000L,
 369         1.339984238037267658352656597960492029261e+0001L,
 370         -1.564707657605373662425785904278645727813e+0002L,
 371         2.156323807499211356127813962223067079300e+0003L,
 372         -3.330486427626223184647299834137041307569e+0004L,
 373         5.235535072011889213611369254140123518699e+0005L,
 374         -7.258160984602220710491988573430212593080e+0006L,
 375         7.316526934569686459641438882340322673357e+0007L,
 376         -3.806450279064900548836571789284896711473e+0008L,
 377 };
 378 
 379 static long double
 380 polytail(long double s)
 381 {
 382         long double t, z;
 383         int i;
 384 
 385         z = s * s;
 386         t = pt[18];
 387 
 388         for (i = 17; i >= 1; i--)
 389                 t = pt[i] + z * t;
 390 
 391         return (pt[0] + s * t);
 392 }
 393 
 394 static long double
 395 neg(long double z, int *signgamlp)
 396 {
 397         long double t, p;
 398 
 399         /* BEGIN CSTYLED */
 400         /*
 401          * written by K.C. Ng,  Feb 2, 1989.
 402          *
 403          * Since
 404          *              -z*G(-z)*G(z) = pi/sin(pi*z),
 405          * we have
 406          *      G(-z) = -pi/(sin(pi*z)*G(z)*z)
 407          *                =  pi/(sin(pi*(-z))*G(z)*z)
 408          * Algorithm
 409          *              z = |z|
 410          *              t = sinpi(z); ...note that when z>2**112, z is an int
 411          *              and hence t=0.
 412          *
 413          *              if (t == 0.0) return 1.0/0.0;
 414          *              if (t< 0.0) *signgamlp = -1; else t= -t;
 415          *              if (z<1.0e-40)       ...tiny z
 416          *                  return -log(z);
 417          *              else
 418          *                  return log(pi/(t*z))-lgamma(z);
 419          *
 420          */
 421         /* END CSTYLED */
 422 
 423         t = sinpil(z);                  /* t := sin(pi*z) */
 424 
 425         if (t == c0)                    /* return   1.0/0.0 =  +INF */
 426                 return (c1 / c0);
 427 
 428         z = -z;
 429 
 430         if (z <= tiny)
 431                 p = -logl(z);
 432         else
 433                 p = logl(pi / (fabsl(t) * z)) - __k_lgammal(z, signgamlp);
 434 
 435         if (t < c0)
 436                 *signgamlp = -1;
 437 
 438         return (p);
 439 }