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