usr/src/lib/libm/common/C/__lgamma.c

   1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   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  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  23  */
  24 /*
  25  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
  26  * Use is subject to license terms.
  27  */
  28
  29 /*
  30  * double __k_lgamma(double x, int *signgamp);
  31  *
  32  * K.C. Ng, March, 1989.
  33  *
  34  * Part of the algorithm is based on W. Cody's lgamma function.
  35  */
  36
  37 #include "libm.h"
  38
  39 static const double
  40 one     = 1.0,
  41 zero    = 0.0,
  42 hln2pi  = 0.9189385332046727417803297,  /* log(2*pi)/2 */
  43 pi      = 3.1415926535897932384626434,
  44 two52   = 4503599627370496.0,           /* 43300000,00000000 (used by sin_pi) */
  45 /*
  46  * Numerator and denominator coefficients for rational minimax Approximation
  47  * P/Q over (0.5,1.5).
  48  */
  49 D1 =    -5.772156649015328605195174e-1,
  50 p7 =     4.945235359296727046734888e0,
  51 p6 =     2.018112620856775083915565e2,
  52 p5 =     2.290838373831346393026739e3,
  53 p4 =     1.131967205903380828685045e4,
  54 p3 =     2.855724635671635335736389e4,
  55 p2 =     3.848496228443793359990269e4,
  56 p1 =     2.637748787624195437963534e4,
  57 p0 =     7.225813979700288197698961e3,
  58 q7 =     6.748212550303777196073036e1,
  59 q6 =     1.113332393857199323513008e3,
  60 q5 =     7.738757056935398733233834e3,
  61 q4 =     2.763987074403340708898585e4,
  62 q3 =     5.499310206226157329794414e4,
  63 q2 =     6.161122180066002127833352e4,
  64 q1 =     3.635127591501940507276287e4,
  65 q0 =     8.785536302431013170870835e3,
  66 /*
  67  * Numerator and denominator coefficients for rational minimax Approximation
  68  * G/H over (1.5,4.0).
  69  */
  70 D2 =     4.227843350984671393993777e-1,
  71 g7 =     4.974607845568932035012064e0,
  72 g6 =     5.424138599891070494101986e2,
  73 g5 =     1.550693864978364947665077e4,
  74 g4 =     1.847932904445632425417223e5,
  75 g3 =     1.088204769468828767498470e6,
  76 g2 =     3.338152967987029735917223e6,
  77 g1 =     5.106661678927352456275255e6,
  78 g0 =     3.074109054850539556250927e6,
  79 h7 =     1.830328399370592604055942e2,
  80 h6 =     7.765049321445005871323047e3,
  81 h5 =     1.331903827966074194402448e5,
  82 h4 =     1.136705821321969608938755e6,
  83 h3 =     5.267964117437946917577538e6,
  84 h2 =     1.346701454311101692290052e7,
  85 h1 =     1.782736530353274213975932e7,
  86 h0 =     9.533095591844353613395747e6,
  87 /*
  88  * Numerator and denominator coefficients for rational minimax Approximation
  89  * U/V over (4.0,12.0).
  90  */
  91 D4 =     1.791759469228055000094023e0,
  92 u7 =     1.474502166059939948905062e4,
  93 u6 =     2.426813369486704502836312e6,
  94 u5 =     1.214755574045093227939592e8,
  95 u4 =     2.663432449630976949898078e9,
  96 u3 =     2.940378956634553899906876e10,
  97 u2 =     1.702665737765398868392998e11,
  98 u1 =     4.926125793377430887588120e11,
  99 u0 =     5.606251856223951465078242e11,
 100 v7 =     2.690530175870899333379843e3,
 101 v6 =     6.393885654300092398984238e5,
 102 v5 =     4.135599930241388052042842e7,
 103 v4 =     1.120872109616147941376570e9,
 104 v3 =     1.488613728678813811542398e10,
 105 v2 =     1.016803586272438228077304e11,
 106 v1 =     3.417476345507377132798597e11,
 107 v0 =     4.463158187419713286462081e11,
 108 /*
 109  * Coefficients for minimax approximation over (12, INF).
 110  */
 111 c5 =    -1.910444077728e-03,
 112 c4 =     8.4171387781295e-04,
 113 c3 =    -5.952379913043012e-04,
 114 c2 =     7.93650793500350248e-04,
 115 c1 =    -2.777777777777681622553e-03,
 116 c0 =     8.333333333333333331554247e-02,
 117 c6 =     5.7083835261e-03;
 118
 119 /*
 120  * Return sin(pi*x).  We assume x is finite and negative, and if it
 121  * is an integer, then the sign of the zero returned doesn't matter.
 122  */
 123 static double
 124 sin_pi(double x) {
 125         double  y, z;
 126         int     n;
 127
 128         y = -x;
 129         if (y <= 0.25)
 130                 return (__k_sin(pi * x, 0.0));
 131         if (y >= two52)
 132                 return (zero);
 133         z = floor(y);
 134         if (y == z)
 135                 return (zero);
 136
 137         /* argument reduction: set y = |x| mod 2 */
 138         y *= 0.5;
 139         y = 2.0 * (y - floor(y));
 140
 141         /* now floor(y * 4) tells which octant y is in */
 142         n = (int)(y * 4.0);
 143         switch (n) {
 144         case 0:
 145                 y = __k_sin(pi * y, 0.0);
 146                 break;
 147         case 1:
 148         case 2:
 149                 y = __k_cos(pi * (0.5 - y), 0.0);
 150                 break;
 151         case 3:
 152         case 4:
 153                 y = __k_sin(pi * (1.0 - y), 0.0);
 154                 break;
 155         case 5:
 156         case 6:
 157                 y = -__k_cos(pi * (y - 1.5), 0.0);
 158                 break;
 159         default:
 160                 y = __k_sin(pi * (y - 2.0), 0.0);
 161                 break;
 162         }
 163         return (-y);
 164 }
 165
 166 static double
 167 neg(double z, int *signgamp) {
 168         double  t, p;
 169
 170         /*
 171          * written by K.C. Ng,  Feb 2, 1989.
 172          *
 173          * Since
 174          *              -z*G(-z)*G(z) = pi/sin(pi*z),
 175          * we have
 176          *      G(-z) = -pi/(sin(pi*z)*G(z)*z)
 177          *            =  pi/(sin(pi*(-z))*G(z)*z)
 178          * Algorithm
 179          *              z = |z|
 180          *              t = sin_pi(z); ...note that when z>2**52, z is an int
 181          *              and hence t=0.
 182          *
 183          *              if (t == 0.0) return 1.0/0.0;
 184          *              if (t< 0.0) *signgamp = -1; else t= -t;
 185          *              if (z+1.0 == 1.0)       ...tiny z
 186          *                  return -log(z);
 187          *              else
 188          *                  return log(pi/(t*z))-__k_lgamma(z, signgamp);
 189          */
 190
 191         t = sin_pi(z);                  /* t := sin(pi*z) */
 192         if (t == zero)                  /* return 1.0/0.0 = +INF */
 193                 return (one / fabs(t));
 194         z = -z;
 195         p = z + one;
 196         if (p == one)
 197                 p = -log(z);
 198         else
 199                 p = log(pi / (fabs(t) * z)) - __k_lgamma(z, signgamp);
 200         if (t < zero)
 201                 *signgamp = -1;
 202         return (p);
 203 }
 204
 205 double
 206 __k_lgamma(double x, int *signgamp) {
 207         double  t, p, q, cr, y;
 208
 209         /* purge off +-inf, NaN and negative arguments */
 210         if (!finite(x))
 211                 return (x * x);
 212         *signgamp = 1;
 213         if (signbit(x))
 214                 return (neg(x, signgamp));
 215
 216         /* lgamma(x) ~ log(1/x) for really tiny x */
 217         t = one + x;
 218         if (t == one) {
 219                 if (x == zero)
 220                         return (one / x);
 221                 return (-log(x));
 222         }
 223
 224         /* for tiny < x < inf */
 225         if (x <= 1.5) {
 226                 if (x < 0.6796875) {
 227                         cr = -log(x);
 228                         y = x;
 229                 } else {
 230                         cr = zero;
 231                         y = x - one;
 232                 }
 233
 234                 if (x <= 0.5 || x >= 0.6796875) {
 235                         if (x == one)
 236                                 return (zero);
 237                         p = p0+y*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
 238                         q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
 239                             (q7+y)))))));
 240                         return (cr+y*(D1+y*(p/q)));
 241                 } else {
 242                         y = x - one;
 243                         p = g0+y*(g1+y*(g2+y*(g3+y*(g4+y*(g5+y*(g6+y*g7))))));
 244                         q = h0+y*(h1+y*(h2+y*(h3+y*(h4+y*(h5+y*(h6+y*
 245                             (h7+y)))))));
 246                         return (cr+y*(D2+y*(p/q)));
 247                 }
 248         } else if (x <= 4.0) {
 249                 if (x == 2.0)
 250                         return (zero);
 251                 y = x - 2.0;
 252                 p = g0+y*(g1+y*(g2+y*(g3+y*(g4+y*(g5+y*(g6+y*g7))))));
 253                 q = h0+y*(h1+y*(h2+y*(h3+y*(h4+y*(h5+y*(h6+y*(h7+y)))))));
 254                 return (y*(D2+y*(p/q)));
 255         } else if (x <= 12.0) {
 256                 y = x - 4.0;
 257                 p = u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*(u6+y*u7))))));
 258                 q = v0+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*(v6+y*(v7-y)))))));
 259                 return (D4+y*(p/q));
 260         } else if (x <= 1.0e17) {               /* x ~< 2**(prec+3) */
 261                 t = one / x;
 262                 y = t * t;
 263                 p = hln2pi+t*(c0+y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*c6))))));
 264                 q = log(x);
 265                 return (x*(q-one)-(0.5*q-p));
 266         } else {                        /* may overflow */
 267                 return (x * (log(x) - 1.0));
 268         }
 269 }