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 2005 Sun Microsystems, Inc. All rights reserved.
23 * Use is subject to license terms.
24 */
25
26
27 /*
28 * double __k_lgamma(double x, int *signgamp);
29 *
30 * K.C. Ng, March, 1989.
31 *
32 * Part of the algorithm is based on W. Cody's lgamma function.
33 */
34
35 #include "libm.h"
36
37 #define D1_000 0x3FF00000
38 #define D2_000 0x40000000
39 #define D8_000 0x40200000
40 #define D0_999999 0x3FECCCCC
41 #define D0_731600 0x3FE76944
42 #define D0_231640 0x3FCDA661
43 #define D1_731631 0x3FFBB4C3
44 #define D1_231632 0x3FF3B4C4
45 #define D2_882E17 0x43900000 /* 2^58 */
46
47 static const double
48 one = 1.0,
49 zero = 0.0,
50 hln2pi = 0.9189385332046727417803297, /* log(2*pi)/2 */
51 pi = 3.1415926535897932384626434,
52 two52 = 4503599627370496.0; /* 43300000,00000000 (used by sin_pi) */
53
54 /*
55 * From Sunpro: Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
56 *
57 * The original implementation did not match mpfr or OpenLibm testsuite results.
58 */
59 static const double
60 half= 5.00000000000000000000e-01,
61 a0 = 7.72156649015328655494e-02,
62 a1 = 3.22467033424113591611e-01,
63 a2 = 6.73523010531292681824e-02,
64 a3 = 2.05808084325167332806e-02,
65 a4 = 7.38555086081402883957e-03,
66 a5 = 2.89051383673415629091e-03,
67 a6 = 1.19270763183362067845e-03,
68 a7 = 5.10069792153511336608e-04,
69 a8 = 2.20862790713908385557e-04,
70 a9 = 1.08011567247583939954e-04,
71 a10 = 2.52144565451257326939e-05,
72 a11 = 4.48640949618915160150e-05,
73 tc = 1.46163214496836224576e+00,
74 tf = -1.21486290535849611461e-01,
75 tt = -3.63867699703950536541e-18, /* tt = -(tail of tf) */
76 t0 = 4.83836122723810047042e-01,
77 t1 = -1.47587722994593911752e-01,
78 t2 = 6.46249402391333854778e-02,
79 t3 = -3.27885410759859649565e-02,
80 t4 = 1.79706750811820387126e-02,
81 t5 = -1.03142241298341437450e-02,
82 t6 = 6.10053870246291332635e-03,
83 t7 = -3.68452016781138256760e-03,
84 t8 = 2.25964780900612472250e-03,
85 t9 = -1.40346469989232843813e-03,
86 t10 = 8.81081882437654011382e-04,
87 t11 = -5.38595305356740546715e-04,
88 t12 = 3.15632070903625950361e-04,
89 t13 = -3.12754168375120860518e-04,
90 t14 = 3.35529192635519073543e-04,
91 u0 = -7.72156649015328655494e-02,
92 u1 = 6.32827064025093366517e-01,
93 u2 = 1.45492250137234768737e+00,
94 u3 = 9.77717527963372745603e-01,
95 u4 = 2.28963728064692451092e-01,
96 u5 = 1.33810918536787660377e-02,
97 v1 = 2.45597793713041134822e+00,
98 v2 = 2.12848976379893395361e+00,
99 v3 = 7.69285150456672783825e-01,
100 v4 = 1.04222645593369134254e-01,
101 v5 = 3.21709242282423911810e-03,
102 s0 = -7.72156649015328655494e-02,
103 s1 = 2.14982415960608852501e-01,
104 s2 = 3.25778796408930981787e-01,
105 s3 = 1.46350472652464452805e-01,
106 s4 = 2.66422703033638609560e-02,
107 s5 = 1.84028451407337715652e-03,
108 s6 = 3.19475326584100867617e-05,
109 r1 = 1.39200533467621045958e+00,
110 r2 = 7.21935547567138069525e-01,
111 r3 = 1.71933865632803078993e-01,
112 r4 = 1.86459191715652901344e-02,
113 r5 = 7.77942496381893596434e-04,
114 r6 = 7.32668430744625636189e-06,
115 w0 = 4.18938533204672725052e-01,
116 w1 = 8.33333333333329678849e-02,
117 w2 = -2.77777777728775536470e-03,
118 w3 = 7.93650558643019558500e-04,
119 w4 = -5.95187557450339963135e-04,
120 w5 = 8.36339918996282139126e-04,
121 w6 = -1.63092934096575273989e-03;
122
123 /*
124 * Return sin(pi*x). We assume x is finite and negative, and if it
125 * is an integer, then the sign of the zero returned doesn't matter.
126 */
127 static double
sin_pi(double x)128 sin_pi(double x) {
129 double y, z;
130 int n;
131
132 y = -x;
133 if (y <= 0.25)
134 return (__k_sin(pi * x, 0.0));
135 if (y >= two52)
136 return (zero);
137 z = floor(y);
138 if (y == z)
139 return (zero);
140
141 /* argument reduction: set y = |x| mod 2 */
142 y *= 0.5;
143 y = 2.0 * (y - floor(y));
144
145 /* now floor(y * 4) tells which octant y is in */
146 n = (int)(y * 4.0);
147 switch (n) {
148 case 0:
149 y = __k_sin(pi * y, 0.0);
150 break;
151 case 1:
152 case 2:
153 y = __k_cos(pi * (0.5 - y), 0.0);
154 break;
155 case 3:
156 case 4:
157 y = __k_sin(pi * (1.0 - y), 0.0);
158 break;
159 case 5:
160 case 6:
161 y = -__k_cos(pi * (y - 1.5), 0.0);
162 break;
163 default:
164 y = __k_sin(pi * (y - 2.0), 0.0);
165 break;
166 }
167 return (-y);
168 }
169
170 static double
neg(double z,int * signgamp)171 neg(double z, int *signgamp) {
172 double t, p;
173
174 /*
175 * written by K.C. Ng, Feb 2, 1989.
176 *
177 * Since
178 * -z*G(-z)*G(z) = pi/sin(pi*z),
179 * we have
180 * G(-z) = -pi/(sin(pi*z)*G(z)*z)
181 * = pi/(sin(pi*(-z))*G(z)*z)
182 * Algorithm
183 * z = |z|
184 * t = sin_pi(z); ...note that when z>2**52, z is an int
185 * and hence t=0.
186 *
187 * if(t==0.0) return 1.0/0.0;
188 * if(t< 0.0) *signgamp = -1; else t= -t;
189 * if(z+1.0==1.0) ...tiny z
190 * return -log(z);
191 * else
192 * return log(pi/(t*z))-__k_lgamma(z, signgamp);
193 */
194
195 t = sin_pi(z); /* t := sin(pi*z) */
196 if (t == zero) /* return 1.0/0.0 = +INF */
197 return (one / fabs(t));
198 z = -z;
199 p = z + one;
200 if (p == one)
201 p = -log(z);
202 else
203 p = log(pi / (fabs(t) * z)) - __k_lgamma(z, signgamp);
204 if (t < zero)
205 *signgamp = -1;
206 return (p);
207 }
208
209 double
__k_lgamma(double x,int * signgamp)210 __k_lgamma(double x, int *signgamp) {
211 double t, p, q, y;
212 double r, w, z, p1, p2, p3;
213 int i;
214 int32_t hx,lx;
215
216 /* purge off +-inf, +-0, NaN and negative arguments */
217 if (!finite(x))
218 return (x * x);
219 *signgamp = 1;
220 if (signbit(x))
221 return (neg(x, signgamp));
222
223 EXTRACT_WORDS(hx,lx,x);
224 hx &= 0x7fffffff;
225
226 if((((hx - D1_000)|lx) == 0) || (((hx - D2_000)|lx) == 0)) {
227 r = 0; /* purge off 1 and 2 */
228 } else if(hx < D2_000) {
229 if(hx <= D0_999999) { /* lgamma(x) = lgamma(x+1)-log(x) */
230 r = -log(x);
231 if(hx >= D0_731600) {y = one - x; i = 0;}
232 else if(hx >= D0_231640) {y = x - (tc - one); i = 1;}
233 else {y = x; i = 2;}
234 } else {
235 r = zero;
236 if(hx >= D1_731631) {y = 2.0 - x; i = 0;}
237 else if(hx >= D1_231632) {y = x - tc; i = 1;}
238 else {y = x - one; i = 2;}
239 }
240 switch(i) {
241 case 0:
242 z = y*y;
243 p1 = a0 + z*(a2 + z*(a4 + z*(a6 + z*(a8 + z*a10))));
244 p2 = z*(a1 + z*(a3 + z*(a5 + z*(a7 + z*(a9 + z*a11)))));
245 p = y*p1 + p2;
246 r += (p - 0.5*y); break;
247 case 1:
248 z = y*y;
249 w = z*y;
250 p1 = t0 + w*(t3 + w*(t6 + w*(t9 + w*t12)));
251 p2 = t1 + w*(t4 + w*(t7 + w*(t10 + w*t13)));
252 p3 = t2 + w*(t5 + w*(t8 + w*(t11 + w*t14)));
253 p = z*p1 - (tt - w*(p2 + y*p3));
254 r += (tf + p); break;
255 case 2:
256 p1 = y*(u0 + y*(u1 + y*(u2 + y*(u3 + y*(u4 + y*u5)))));
257 p2 = one + y*(v1 + y*(v2 + y*(v3 + y*(v4 + y*v5))));
258 r += (-0.5*y + p1/p2);
259 }
260 } else if (hx < D8_000) {
261 i = (int)x;
262 y = x - (double)i;
263 p = y*(s0 + y*(s1 + y*(s2 + y*(s3 + y*(s4 + y*(s5 + y*s6))))));
264 q = one + y*(r1 + y*(r2 + y*(r3 + y*(r4 + y*(r5 + y*r6)))));
265 r = half*y + p/q;
266 z = one; /* lgamma(1+s) = log(s) + lgamma(s) */
267 switch(i) {
268 case 7: z *= (y + 6.0); /* FALLTHRU */
269 case 6: z *= (y + 5.0); /* FALLTHRU */
270 case 5: z *= (y + 4.0); /* FALLTHRU */
271 case 4: z *= (y + 3.0); /* FALLTHRU */
272 case 3: z *= (y + 2.0); /* FALLTHRU */
273 r += log(z); break;
274 }
275 } else if (hx < D2_882E17) {
276 t = log(x);
277 z = one/x;
278 y = z*z;
279 w = w0 + z*(w1 + y*(w2 + y*(w3 + y*(w4 + y*(w5 + y*w6)))));
280 r = (x - half)*(t - one) + w;
281 } else {
282 r = x*(log(x) - one);
283 }
284 return r;
285 }
286