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 }