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 } 270