Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/C/__lgamma.c
+++ new/usr/src/lib/libm/common/C/__lgamma.c
1 1 /*
2 2 * CDDL HEADER START
3 3 *
4 4 * The contents of this file are subject to the terms of the
5 5 * Common Development and Distribution License (the "License").
6 6 * You may not use this file except in compliance with the License.
7 7 *
8 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 9 * or http://www.opensolaris.org/os/licensing.
10 10 * See the License for the specific language governing permissions
|
↓ open down ↓ |
10 lines elided |
↑ open up ↑ |
11 11 * and limitations under the License.
12 12 *
13 13 * When distributing Covered Code, include this CDDL HEADER in each
14 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 15 * If applicable, add the following below this CDDL HEADER, with the
16 16 * fields enclosed by brackets "[]" replaced with your own identifying
17 17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 18 *
19 19 * CDDL HEADER END
20 20 */
21 +
21 22 /*
22 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
23 24 */
25 +
24 26 /*
25 27 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
26 28 * Use is subject to license terms.
27 29 */
28 30
29 31 /*
30 32 * double __k_lgamma(double x, int *signgamp);
31 33 *
32 34 * K.C. Ng, March, 1989.
33 35 *
34 36 * Part of the algorithm is based on W. Cody's lgamma function.
35 37 */
36 38
37 39 #include "libm.h"
38 40
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) */
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 +
45 47 /*
46 48 * Numerator and denominator coefficients for rational minimax Approximation
47 49 * P/Q over (0.5,1.5).
48 50 */
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,
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 +
66 69 /*
67 70 * Numerator and denominator coefficients for rational minimax Approximation
68 71 * G/H over (1.5,4.0).
69 72 */
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,
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 +
87 91 /*
88 92 * Numerator and denominator coefficients for rational minimax Approximation
89 93 * U/V over (4.0,12.0).
90 94 */
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,
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 +
108 113 /*
109 114 * Coefficients for minimax approximation over (12, INF).
110 115 */
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;
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;
118 123
119 124 /*
120 125 * Return sin(pi*x). We assume x is finite and negative, and if it
121 126 * is an integer, then the sign of the zero returned doesn't matter.
122 127 */
123 128 static double
124 -sin_pi(double x) {
125 - double y, z;
126 - int n;
129 +sin_pi(double x)
130 +{
131 + double y, z;
132 + int n;
127 133
128 134 y = -x;
135 +
129 136 if (y <= 0.25)
130 137 return (__k_sin(pi * x, 0.0));
138 +
131 139 if (y >= two52)
132 140 return (zero);
141 +
133 142 z = floor(y);
143 +
134 144 if (y == z)
135 145 return (zero);
136 146
137 147 /* argument reduction: set y = |x| mod 2 */
138 148 y *= 0.5;
139 149 y = 2.0 * (y - floor(y));
140 150
141 151 /* now floor(y * 4) tells which octant y is in */
142 152 n = (int)(y * 4.0);
153 +
143 154 switch (n) {
144 155 case 0:
145 156 y = __k_sin(pi * y, 0.0);
146 157 break;
147 158 case 1:
148 159 case 2:
149 160 y = __k_cos(pi * (0.5 - y), 0.0);
150 161 break;
151 162 case 3:
152 163 case 4:
153 164 y = __k_sin(pi * (1.0 - y), 0.0);
154 165 break;
155 166 case 5:
156 167 case 6:
157 168 y = -__k_cos(pi * (y - 1.5), 0.0);
158 169 break;
159 170 default:
160 171 y = __k_sin(pi * (y - 2.0), 0.0);
161 172 break;
162 173 }
174 +
163 175 return (-y);
164 176 }
165 177
166 178 static double
167 -neg(double z, int *signgamp) {
168 - double t, p;
179 +neg(double z, int *signgamp)
180 +{
181 + double t, p;
169 182
183 + /* BEGIN CSTYLED */
170 184 /*
171 185 * written by K.C. Ng, Feb 2, 1989.
172 186 *
173 187 * Since
174 188 * -z*G(-z)*G(z) = pi/sin(pi*z),
175 189 * we have
176 - * G(-z) = -pi/(sin(pi*z)*G(z)*z)
177 - * = pi/(sin(pi*(-z))*G(z)*z)
190 + * G(-z) = -pi/(sin(pi*z)*G(z)*z)
191 + * = pi/(sin(pi*(-z))*G(z)*z)
178 192 * Algorithm
179 193 * z = |z|
180 194 * t = sin_pi(z); ...note that when z>2**52, z is an int
181 195 * and hence t=0.
182 196 *
183 197 * if (t == 0.0) return 1.0/0.0;
184 198 * if (t< 0.0) *signgamp = -1; else t= -t;
185 199 * if (z+1.0 == 1.0) ...tiny z
186 200 * return -log(z);
187 201 * else
188 202 * return log(pi/(t*z))-__k_lgamma(z, signgamp);
189 203 */
204 + /* END CSTYLED */
190 205
191 206 t = sin_pi(z); /* t := sin(pi*z) */
207 +
192 208 if (t == zero) /* return 1.0/0.0 = +INF */
193 209 return (one / fabs(t));
210 +
194 211 z = -z;
195 212 p = z + one;
213 +
196 214 if (p == one)
197 215 p = -log(z);
198 216 else
199 217 p = log(pi / (fabs(t) * z)) - __k_lgamma(z, signgamp);
218 +
200 219 if (t < zero)
201 220 *signgamp = -1;
221 +
202 222 return (p);
203 223 }
204 224
205 225 double
206 -__k_lgamma(double x, int *signgamp) {
207 - double t, p, q, cr, y;
226 +__k_lgamma(double x, int *signgamp)
227 +{
228 + double t, p, q, cr, y;
208 229
209 230 /* purge off +-inf, NaN and negative arguments */
210 231 if (!finite(x))
211 232 return (x * x);
233 +
212 234 *signgamp = 1;
235 +
213 236 if (signbit(x))
214 237 return (neg(x, signgamp));
215 238
216 239 /* lgamma(x) ~ log(1/x) for really tiny x */
217 240 t = one + x;
241 +
218 242 if (t == one) {
219 243 if (x == zero)
220 244 return (one / x);
245 +
221 246 return (-log(x));
222 247 }
223 248
224 249 /* for tiny < x < inf */
225 250 if (x <= 1.5) {
226 251 if (x < 0.6796875) {
227 252 cr = -log(x);
228 253 y = x;
229 254 } else {
230 255 cr = zero;
231 256 y = x - one;
232 257 }
233 258
234 259 if (x <= 0.5 || x >= 0.6796875) {
235 260 if (x == one)
236 261 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)));
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)));
241 268 } else {
242 269 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)));
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)));
247 275 }
248 276 } else if (x <= 4.0) {
249 277 if (x == 2.0)
250 278 return (zero);
279 +
251 280 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)));
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)));
255 286 } else if (x <= 12.0) {
256 287 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) */
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) */
261 294 t = one / x;
262 295 y = t * t;
263 - p = hln2pi+t*(c0+y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*c6))))));
296 + p = hln2pi + t * (c0 + y * (c1 + y * (c2 + y * (c3 + y * (c4 +
297 + y * (c5 + y * c6))))));
264 298 q = log(x);
265 - return (x*(q-one)-(0.5*q-p));
299 + return (x * (q - one) - (0.5 * q - p));
266 300 } else { /* may overflow */
267 301 return (x * (log(x) - 1.0));
268 302 }
269 303 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX