Print this page
11210 libm should be cstyle(1ONBLD) clean
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/m9x/tgammaf.c
+++ new/usr/src/lib/libm/common/m9x/tgammaf.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
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.
|
↓ open down ↓ |
14 lines elided |
↑ open up ↑ |
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 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 24 */
25 +
25 26 /*
26 27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 28 * Use is subject to license terms.
28 29 */
29 30
30 31 #pragma weak __tgammaf = tgammaf
31 32
32 33 /*
33 34 * True gamma function
34 35 *
35 36 * float tgammaf(float x)
36 37 *
37 38 * Algorithm: see tgamma.c
38 39 *
39 40 * Maximum error observed: 0.87ulp (both positive and negative arguments)
|
↓ open down ↓ |
5 lines elided |
↑ open up ↑ |
40 41 */
41 42
42 43 #include "libm.h"
43 44 #include <math.h>
44 45 #if defined(__SUNPRO_C)
45 46 #include <sunmath.h>
46 47 #endif
47 48 #include <sys/isa_defs.h>
48 49
49 50 #if defined(_BIG_ENDIAN)
50 -#define HIWORD 0
51 -#define LOWORD 1
51 +#define HIWORD 0
52 +#define LOWORD 1
52 53 #else
53 -#define HIWORD 1
54 -#define LOWORD 0
54 +#define HIWORD 1
55 +#define LOWORD 0
55 56 #endif
56 -#define __HI(x) ((int *) &x)[HIWORD]
57 -#define __LO(x) ((unsigned *) &x)[LOWORD]
57 +#define __HI(x) ((int *)&x)[HIWORD]
58 +#define __LO(x) ((unsigned *)&x)[LOWORD]
58 59
59 60 /* Coefficients for primary intervals GTi() */
60 61 static const double cr[] = {
61 62 /* p1 */
62 63 +7.09087253435088360271451613398019280077561279443e-0001,
63 64 -5.17229560788652108545141978238701790105241761089e-0001,
64 65 +5.23403394528150789405825222323770647162337764327e-0001,
65 66 -4.54586308717075010784041566069480411732634814899e-0001,
66 67 +4.20596490915239085459964590559256913498190955233e-0001,
67 68 -3.57307589712377520978332185838241458642142185789e-0001,
68 -
69 69 /* p2 */
70 70 +4.28486983980295198166056119223984284434264344578e-0001,
71 71 -1.30704539487709138528680121627899735386650103914e-0001,
72 72 +1.60856285038051955072861219352655851542955430871e-0001,
73 73 -9.22285161346010583774458802067371182158937943507e-0002,
74 74 +7.19240511767225260740890292605070595560626179357e-0002,
75 75 -4.88158265593355093703112238534484636193260459574e-0002,
76 -
77 76 /* p3 */
78 77 +3.82409531118807759081121479786092134814808872880e-0001,
79 78 +2.65309888180188647956400403013495759365167853426e-0002,
80 79 +8.06815109775079171923561169415370309376296739835e-0002,
81 80 -1.54821591666137613928840890835174351674007764799e-0002,
82 81 +1.76308239242717268530498313416899188157165183405e-0002,
83 -
84 82 /* GZi and TZi */
85 83 +0.9382046279096824494097535615803269576988, /* GZ1 */
86 84 +0.8856031944108887002788159005825887332080, /* GZ2 */
87 85 +0.9367814114636523216188468970808378497426, /* GZ3 */
88 - -0.3517214357852935791015625, /* TZ1 */
89 - +0.280530631542205810546875, /* TZ3 */
86 + -0.3517214357852935791015625, /* TZ1 */
87 + +0.280530631542205810546875, /* TZ3 */
90 88 };
91 89
92 -#define P10 cr[0]
93 -#define P11 cr[1]
94 -#define P12 cr[2]
95 -#define P13 cr[3]
96 -#define P14 cr[4]
97 -#define P15 cr[5]
98 -#define P20 cr[6]
99 -#define P21 cr[7]
100 -#define P22 cr[8]
101 -#define P23 cr[9]
102 -#define P24 cr[10]
103 -#define P25 cr[11]
104 -#define P30 cr[12]
105 -#define P31 cr[13]
106 -#define P32 cr[14]
107 -#define P33 cr[15]
108 -#define P34 cr[16]
109 -#define GZ1 cr[17]
110 -#define GZ2 cr[18]
111 -#define GZ3 cr[19]
112 -#define TZ1 cr[20]
113 -#define TZ3 cr[21]
90 +#define P10 cr[0]
91 +#define P11 cr[1]
92 +#define P12 cr[2]
93 +#define P13 cr[3]
94 +#define P14 cr[4]
95 +#define P15 cr[5]
96 +#define P20 cr[6]
97 +#define P21 cr[7]
98 +#define P22 cr[8]
99 +#define P23 cr[9]
100 +#define P24 cr[10]
101 +#define P25 cr[11]
102 +#define P30 cr[12]
103 +#define P31 cr[13]
104 +#define P32 cr[14]
105 +#define P33 cr[15]
106 +#define P34 cr[16]
107 +#define GZ1 cr[17]
108 +#define GZ2 cr[18]
109 +#define GZ3 cr[19]
110 +#define TZ1 cr[20]
111 +#define TZ3 cr[21]
114 112
115 113 /* compute gamma(y) for y in GT1 = [1.0000, 1.2845] */
116 114 static double
117 -GT1(double y) {
115 +GT1(double y)
116 +{
118 117 double z, r;
119 118
120 119 z = y * y;
121 120 r = TZ1 * y + z * ((P10 + y * P11 + z * P12) + (z * y) * (P13 + y *
122 - P14 + z * P15));
121 + P14 + z * P15));
123 122 return (GZ1 + r);
124 123 }
125 124
126 125 /* compute gamma(y) for y in GT2 = [1.2844, 1.6374] */
127 126 static double
128 -GT2(double y) {
127 +GT2(double y)
128 +{
129 129 double z;
130 130
131 131 z = y * y;
132 132 return (GZ2 + z * ((P20 + y * P21 + z * P22) + (z * y) * (P23 + y *
133 - P24 + z * P25)));
133 + P24 + z * P25)));
134 134 }
135 135
136 136 /* compute gamma(y) for y in GT3 = [1.6373, 2.0000] */
137 137 static double
138 -GT3(double y) {
139 -double z, r;
138 +GT3(double y)
139 +{
140 + double z, r;
140 141
141 142 z = y * y;
142 143 r = TZ3 * y + z * ((P30 + y * P31 + z * P32) + (z * y) * (P33 + y *
143 - P34));
144 + P34));
144 145 return (GZ3 + r);
145 146 }
146 147
147 -/* INDENT OFF */
148 148 static const double c[] = {
149 -+1.0,
150 -+2.0,
151 -+0.5,
152 -+1.0e-300,
153 -+6.666717231848518054693623697539230e-0001, /* A1=T3[0] */
154 -+8.33333330959694065245736888749042811909994573178e-0002, /* GP[0] */
155 --2.77765545601667179767706600890361535225507762168e-0003, /* GP[1] */
156 -+7.77830853479775281781085278324621033523037489883e-0004, /* GP[2] */
157 -+4.18938533204672741744150788368695779923320328369e-0001, /* hln2pi */
158 -+2.16608493924982901946e-02, /* ln2_32 */
159 -+4.61662413084468283841e+01, /* invln2_32 */
160 -+5.00004103388988968841156421415669985414073453720e-0001, /* Et1 */
161 -+1.66667656752800761782778277828110208108687545908e-0001, /* Et2 */
149 + +1.0,
150 + +2.0,
151 + +0.5,
152 + +1.0e-300,
153 + +6.666717231848518054693623697539230e-0001, /* A1=T3[0] */
154 + +8.33333330959694065245736888749042811909994573178e-0002, /* GP[0] */
155 + -2.77765545601667179767706600890361535225507762168e-0003, /* GP[1] */
156 + +7.77830853479775281781085278324621033523037489883e-0004, /* GP[2] */
157 + +4.18938533204672741744150788368695779923320328369e-0001, /* hln2pi */
158 + +2.16608493924982901946e-02, /* ln2_32 */
159 + +4.61662413084468283841e+01, /* invln2_32 */
160 + +5.00004103388988968841156421415669985414073453720e-0001, /* Et1 */
161 + +1.66667656752800761782778277828110208108687545908e-0001, /* Et2 */
162 162 };
163 163
164 164 #define one c[0]
165 165 #define two c[1]
166 166 #define half c[2]
167 167 #define tiny c[3]
168 168 #define A1 c[4]
169 169 #define GP0 c[5]
170 170 #define GP1 c[6]
171 171 #define GP2 c[7]
172 172 #define hln2pi c[8]
173 173 #define ln2_32 c[9]
174 174 #define invln2_32 c[10]
175 175 #define Et1 c[11]
176 176 #define Et2 c[12]
177 177
178 178 /* S[j] = 2**(j/32.) for the final computation of exp(w) */
179 179 static const double S[] = {
180 -+1.00000000000000000000e+00, /* 3FF0000000000000 */
181 -+1.02189714865411662714e+00, /* 3FF059B0D3158574 */
182 -+1.04427378242741375480e+00, /* 3FF0B5586CF9890F */
183 -+1.06714040067682369717e+00, /* 3FF11301D0125B51 */
184 -+1.09050773266525768967e+00, /* 3FF172B83C7D517B */
185 -+1.11438674259589243221e+00, /* 3FF1D4873168B9AA */
186 -+1.13878863475669156458e+00, /* 3FF2387A6E756238 */
187 -+1.16372485877757747552e+00, /* 3FF29E9DF51FDEE1 */
188 -+1.18920711500272102690e+00, /* 3FF306FE0A31B715 */
189 -+1.21524735998046895524e+00, /* 3FF371A7373AA9CB */
190 -+1.24185781207348400201e+00, /* 3FF3DEA64C123422 */
191 -+1.26905095719173321989e+00, /* 3FF44E086061892D */
192 -+1.29683955465100964055e+00, /* 3FF4BFDAD5362A27 */
193 -+1.32523664315974132322e+00, /* 3FF5342B569D4F82 */
194 -+1.35425554693689265129e+00, /* 3FF5AB07DD485429 */
195 -+1.38390988196383202258e+00, /* 3FF6247EB03A5585 */
196 -+1.41421356237309514547e+00, /* 3FF6A09E667F3BCD */
197 -+1.44518080697704665027e+00, /* 3FF71F75E8EC5F74 */
198 -+1.47682614593949934623e+00, /* 3FF7A11473EB0187 */
199 -+1.50916442759342284141e+00, /* 3FF82589994CCE13 */
200 -+1.54221082540794074411e+00, /* 3FF8ACE5422AA0DB */
201 -+1.57598084510788649659e+00, /* 3FF93737B0CDC5E5 */
202 -+1.61049033194925428347e+00, /* 3FF9C49182A3F090 */
203 -+1.64575547815396494578e+00, /* 3FFA5503B23E255D */
204 -+1.68179283050742900407e+00, /* 3FFAE89F995AD3AD */
205 -+1.71861929812247793414e+00, /* 3FFB7F76F2FB5E47 */
206 -+1.75625216037329945351e+00, /* 3FFC199BDD85529C */
207 -+1.79470907500310716820e+00, /* 3FFCB720DCEF9069 */
208 -+1.83400808640934243066e+00, /* 3FFD5818DCFBA487 */
209 -+1.87416763411029996256e+00, /* 3FFDFC97337B9B5F */
210 -+1.91520656139714740007e+00, /* 3FFEA4AFA2A490DA */
211 -+1.95714412417540017941e+00, /* 3FFF50765B6E4540 */
180 + +1.00000000000000000000e+00, /* 3FF0000000000000 */
181 + +1.02189714865411662714e+00, /* 3FF059B0D3158574 */
182 + +1.04427378242741375480e+00, /* 3FF0B5586CF9890F */
183 + +1.06714040067682369717e+00, /* 3FF11301D0125B51 */
184 + +1.09050773266525768967e+00, /* 3FF172B83C7D517B */
185 + +1.11438674259589243221e+00, /* 3FF1D4873168B9AA */
186 + +1.13878863475669156458e+00, /* 3FF2387A6E756238 */
187 + +1.16372485877757747552e+00, /* 3FF29E9DF51FDEE1 */
188 + +1.18920711500272102690e+00, /* 3FF306FE0A31B715 */
189 + +1.21524735998046895524e+00, /* 3FF371A7373AA9CB */
190 + +1.24185781207348400201e+00, /* 3FF3DEA64C123422 */
191 + +1.26905095719173321989e+00, /* 3FF44E086061892D */
192 + +1.29683955465100964055e+00, /* 3FF4BFDAD5362A27 */
193 + +1.32523664315974132322e+00, /* 3FF5342B569D4F82 */
194 + +1.35425554693689265129e+00, /* 3FF5AB07DD485429 */
195 + +1.38390988196383202258e+00, /* 3FF6247EB03A5585 */
196 + +1.41421356237309514547e+00, /* 3FF6A09E667F3BCD */
197 + +1.44518080697704665027e+00, /* 3FF71F75E8EC5F74 */
198 + +1.47682614593949934623e+00, /* 3FF7A11473EB0187 */
199 + +1.50916442759342284141e+00, /* 3FF82589994CCE13 */
200 + +1.54221082540794074411e+00, /* 3FF8ACE5422AA0DB */
201 + +1.57598084510788649659e+00, /* 3FF93737B0CDC5E5 */
202 + +1.61049033194925428347e+00, /* 3FF9C49182A3F090 */
203 + +1.64575547815396494578e+00, /* 3FFA5503B23E255D */
204 + +1.68179283050742900407e+00, /* 3FFAE89F995AD3AD */
205 + +1.71861929812247793414e+00, /* 3FFB7F76F2FB5E47 */
206 + +1.75625216037329945351e+00, /* 3FFC199BDD85529C */
207 + +1.79470907500310716820e+00, /* 3FFCB720DCEF9069 */
208 + +1.83400808640934243066e+00, /* 3FFD5818DCFBA487 */
209 + +1.87416763411029996256e+00, /* 3FFDFC97337B9B5F */
210 + +1.91520656139714740007e+00, /* 3FFEA4AFA2A490DA */
211 + +1.95714412417540017941e+00, /* 3FFF50765B6E4540 */
212 212 };
213 -/* INDENT ON */
214 213
215 -/* INDENT OFF */
214 +
216 215 /*
217 216 * return tgammaf(x) in double for 8<x<=35.040096283... using Stirling's formula
218 217 * log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + (1/x)*P(1/(x*x))
219 218 */
219 +
220 220 /*
221 221 * compute ss = log(x)-1
222 222 *
223 223 * log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y, 1<=y<2,
224 224 * j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
225 225 * T1(n-3) = n*log(2)-1, n=3,4,5
226 226 * T2(j) = log(z[j]),
227 227 * T3(s) = 2s + A1*s^3
228 228 * Note
229 229 * (1) Remez error for T3(s) is bounded by 2**(-35.8)
230 230 * (see mpremez/work/Log/tgamma_log_2_outr1)
231 231 */
232 -
233 -static const double T1[] = { /* T1[j]=(j+3)*log(2)-1 */
234 -+1.079441541679835928251696364375e+00,
235 -+1.772588722239781237668928485833e+00,
236 -+2.465735902799726547086160607291e+00,
232 +static const double T1[] = { /* T1[j]=(j+3)*log(2)-1 */
233 + +1.079441541679835928251696364375e+00,
234 + +1.772588722239781237668928485833e+00,
235 + +2.465735902799726547086160607291e+00,
237 236 };
238 237
239 -static const double T2[] = { /* T2[j]=log(1+j/64+1/128) */
240 -+7.782140442054948947462900061137e-03,
241 -+2.316705928153437822879916096229e-02,
242 -+3.831886430213659919375532512380e-02,
243 -+5.324451451881228286587019378653e-02,
244 -+6.795066190850774939456527777263e-02,
245 -+8.244366921107459126816006866831e-02,
246 -+9.672962645855111229557105648746e-02,
247 -+1.108143663402901141948061693232e-01,
248 -+1.247034785009572358634065153809e-01,
249 -+1.384023228591191356853258736016e-01,
250 -+1.519160420258419750718034248969e-01,
251 -+1.652495728953071628756114492772e-01,
252 -+1.784076574728182971194002415109e-01,
253 -+1.913948529996294546092988075613e-01,
254 -+2.042155414286908915038203861962e-01,
255 -+2.168739383006143596190895257443e-01,
256 -+2.293741010648458299914807250461e-01,
257 -+2.417199368871451681443075159135e-01,
258 -+2.539152099809634441373232979066e-01,
259 -+2.659635484971379413391259265375e-01,
260 -+2.778684510034563061863500329234e-01,
261 -+2.896332925830426768788930555257e-01,
262 -+3.012613305781617810128755382338e-01,
263 -+3.127557100038968883862465596883e-01,
264 -+3.241194686542119760906707604350e-01,
265 -+3.353555419211378302571795798142e-01,
266 -+3.464667673462085809184621884258e-01,
267 -+3.574558889218037742260094901409e-01,
268 -+3.683255611587076530482301540504e-01,
269 -+3.790783529349694583908533456310e-01,
270 -+3.897167511400252133704636040035e-01,
271 -+4.002431641270127069293251019951e-01,
272 -+4.106599249852683859343062031758e-01,
273 -+4.209692946441296361288671615068e-01,
274 -+4.311734648183713408591724789556e-01,
275 -+4.412745608048752294894964416613e-01,
276 -+4.512746441394585851446923830790e-01,
277 -+4.611757151221701663679999255979e-01,
278 -+4.709797152187910125468978560564e-01,
279 -+4.806885293457519076766184554480e-01,
280 -+4.903039880451938381503461596457e-01,
281 -+4.998278695564493298213314152470e-01,
282 -+5.092619017898079468040749192283e-01,
283 -+5.186077642080456321529769963648e-01,
284 -+5.278670896208423851138922177783e-01,
285 -+5.370414658968836545667292441538e-01,
286 -+5.461324375981356503823972092312e-01,
287 -+5.551415075405015927154803595159e-01,
288 -+5.640701382848029660713842900902e-01,
289 -+5.729197535617855090927567266263e-01,
290 -+5.816917396346224825206107537254e-01,
291 -+5.903874466021763746419167081236e-01,
292 -+5.990081896460833993816000244617e-01,
293 -+6.075552502245417955010851527911e-01,
294 -+6.160298772155140196475659281967e-01,
295 -+6.244332880118935010425387440547e-01,
296 -+6.327666695710378295457864685036e-01,
297 -+6.410311794209312910556013344054e-01,
298 -+6.492279466251098188908399699053e-01,
299 -+6.573580727083600301418900232459e-01,
300 -+6.654226325450904489500926100067e-01,
301 -+6.734226752121667202979603888010e-01,
302 -+6.813592248079030689480715595681e-01,
303 -+6.892332812388089803249143378146e-01,
238 +static const double T2[] = { /* T2[j]=log(1+j/64+1/128) */
239 + +7.782140442054948947462900061137e-03,
240 + +2.316705928153437822879916096229e-02,
241 + +3.831886430213659919375532512380e-02,
242 + +5.324451451881228286587019378653e-02,
243 + +6.795066190850774939456527777263e-02,
244 + +8.244366921107459126816006866831e-02,
245 + +9.672962645855111229557105648746e-02,
246 + +1.108143663402901141948061693232e-01,
247 + +1.247034785009572358634065153809e-01,
248 + +1.384023228591191356853258736016e-01,
249 + +1.519160420258419750718034248969e-01,
250 + +1.652495728953071628756114492772e-01,
251 + +1.784076574728182971194002415109e-01,
252 + +1.913948529996294546092988075613e-01,
253 + +2.042155414286908915038203861962e-01,
254 + +2.168739383006143596190895257443e-01,
255 + +2.293741010648458299914807250461e-01,
256 + +2.417199368871451681443075159135e-01,
257 + +2.539152099809634441373232979066e-01,
258 + +2.659635484971379413391259265375e-01,
259 + +2.778684510034563061863500329234e-01,
260 + +2.896332925830426768788930555257e-01,
261 + +3.012613305781617810128755382338e-01,
262 + +3.127557100038968883862465596883e-01,
263 + +3.241194686542119760906707604350e-01,
264 + +3.353555419211378302571795798142e-01,
265 + +3.464667673462085809184621884258e-01,
266 + +3.574558889218037742260094901409e-01,
267 + +3.683255611587076530482301540504e-01,
268 + +3.790783529349694583908533456310e-01,
269 + +3.897167511400252133704636040035e-01,
270 + +4.002431641270127069293251019951e-01,
271 + +4.106599249852683859343062031758e-01,
272 + +4.209692946441296361288671615068e-01,
273 + +4.311734648183713408591724789556e-01,
274 + +4.412745608048752294894964416613e-01,
275 + +4.512746441394585851446923830790e-01,
276 + +4.611757151221701663679999255979e-01,
277 + +4.709797152187910125468978560564e-01,
278 + +4.806885293457519076766184554480e-01,
279 + +4.903039880451938381503461596457e-01,
280 + +4.998278695564493298213314152470e-01,
281 + +5.092619017898079468040749192283e-01,
282 + +5.186077642080456321529769963648e-01,
283 + +5.278670896208423851138922177783e-01,
284 + +5.370414658968836545667292441538e-01,
285 + +5.461324375981356503823972092312e-01,
286 + +5.551415075405015927154803595159e-01,
287 + +5.640701382848029660713842900902e-01,
288 + +5.729197535617855090927567266263e-01,
289 + +5.816917396346224825206107537254e-01,
290 + +5.903874466021763746419167081236e-01,
291 + +5.990081896460833993816000244617e-01,
292 + +6.075552502245417955010851527911e-01,
293 + +6.160298772155140196475659281967e-01,
294 + +6.244332880118935010425387440547e-01,
295 + +6.327666695710378295457864685036e-01,
296 + +6.410311794209312910556013344054e-01,
297 + +6.492279466251098188908399699053e-01,
298 + +6.573580727083600301418900232459e-01,
299 + +6.654226325450904489500926100067e-01,
300 + +6.734226752121667202979603888010e-01,
301 + +6.813592248079030689480715595681e-01,
302 + +6.892332812388089803249143378146e-01,
304 303 };
305 -/* INDENT ON */
306 304
307 305 static double
308 -large_gam(double x) {
306 +large_gam(double x)
307 +{
309 308 double ss, zz, z, t1, t2, w, y, u;
310 309 unsigned lx;
311 310 int k, ix, j, m;
312 311
313 312 ix = __HI(x);
314 313 lx = __LO(x);
315 314 m = (ix >> 20) - 0x3ff; /* exponent of x, range:3-5 */
316 315 ix = (ix & 0x000fffff) | 0x3ff00000; /* y = scale x to [1,2] */
317 316 __HI(y) = ix;
318 317 __LO(y) = lx;
319 318 __HI(z) = (ix & 0xffffc000) | 0x2000; /* z[j]=1+j/64+1/128 */
320 319 __LO(z) = 0;
321 320 j = (ix >> 14) & 0x3f;
322 321 t1 = y + z;
323 322 t2 = y - z;
324 323 u = t2 / t1;
325 324 ss = T1[m - 3] + T2[j] + u * (two + A1 * (u * u));
326 - /* ss = log(x)-1 */
325 + /* ss = log(x)-1 */
326 +
327 327 /*
328 328 * compute ww = (x-.5)*(log(x)-1) + .5*(log(2pi)-1) + 1/x*(P(1/x^2)))
329 329 * where ss = log(x) - 1
330 330 */
331 331 z = one / x;
332 332 zz = z * z;
333 333 w = ((x - half) * ss + hln2pi) + z * (GP0 + zz * GP1 + (zz * zz) * GP2);
334 - k = (int) (w * invln2_32 + half);
334 + k = (int)(w * invln2_32 + half);
335 335
336 336 /* compute the exponential of w */
337 337 j = k & 0x1f;
338 338 m = k >> 5;
339 - z = w - (double) k *ln2_32;
339 + z = w - (double)k * ln2_32;
340 340 zz = S[j] * (one + z + (z * z) * (Et1 + z * Et2));
341 341 __HI(zz) += m << 20;
342 342 return (zz);
343 343 }
344 -/* INDENT OFF */
344 +
345 +
345 346 /*
346 347 * kpsin(x)= sin(pi*x)/pi
347 348 * 3 5 7 9
348 349 * = x+ks[0]*x +ks[1]*x +ks[2]*x +ks[3]*x
349 350 */
350 351 static const double ks[] = {
351 --1.64493404985645811354476665052005342839447790544e+0000,
352 -+8.11740794458351064092797249069438269367389272270e-0001,
353 --1.90703144603551216933075809162889536878854055202e-0001,
354 -+2.55742333994264563281155312271481108635575331201e-0002,
352 + -1.64493404985645811354476665052005342839447790544e+0000,
353 + +8.11740794458351064092797249069438269367389272270e-0001,
354 + -1.90703144603551216933075809162889536878854055202e-0001,
355 + +2.55742333994264563281155312271481108635575331201e-0002,
355 356 };
356 -/* INDENT ON */
357 357
358 358 static double
359 -kpsin(double x) {
359 +kpsin(double x)
360 +{
360 361 double z;
361 362
362 363 z = x * x;
363 364 return (x + (x * z) * ((ks[0] + z * ks[1]) + (z * z) * (ks[2] + z *
364 - ks[3])));
365 + ks[3])));
365 366 }
366 367
367 -/* INDENT OFF */
368 +
368 369 /*
369 370 * kpcos(x)= cos(pi*x)/pi
370 371 * 2 4 6
371 372 * = kc[0]+kc[1]*x +kc[2]*x +kc[3]*x
372 373 */
373 374 static const double kc[] = {
374 -+3.18309886183790671537767526745028724068919291480e-0001,
375 --1.57079581447762568199467875065854538626594937791e+0000,
376 -+1.29183528092558692844073004029568674027807393862e+0000,
377 --4.20232949771307685981015914425195471602739075537e-0001,
375 + +3.18309886183790671537767526745028724068919291480e-0001,
376 + -1.57079581447762568199467875065854538626594937791e+0000,
377 + +1.29183528092558692844073004029568674027807393862e+0000,
378 + -4.20232949771307685981015914425195471602739075537e-0001,
378 379 };
379 -/* INDENT ON */
380 380
381 381 static double
382 -kpcos(double x) {
382 +kpcos(double x)
383 +{
383 384 double z;
384 385
385 386 z = x * x;
386 387 return (kc[0] + z * (kc[1] + z * kc[2] + (z * z) * kc[3]));
387 388 }
388 389
389 -/* INDENT OFF */
390 -static const double
391 -t0z1 = 0.134861805732790769689793935774652917006,
392 -t0z2 = 0.461632144968362341262659542325721328468,
393 -t0z3 = 0.819773101100500601787868704921606996312;
394 - /* 1.134861805732790769689793935774652917006 */
395 -/* INDENT ON */
390 +static const double t0z1 = 0.134861805732790769689793935774652917006,
391 + t0z2 = 0.461632144968362341262659542325721328468,
392 + t0z3 = 0.819773101100500601787868704921606996312;
393 +
394 +/*
395 + * 1.134861805732790769689793935774652917006
396 + */
396 397
397 398 /*
398 399 * gamma(x+i) for 0 <= x < 1
399 400 */
400 401 static double
401 -gam_n(int i, double x) {
402 +gam_n(int i, double x)
403 +{
402 404 double rr = 0.0L, yy;
403 405 double z1, z2;
404 406
405 407 /* compute yy = gamma(x+1) */
406 408 if (x > 0.2845) {
407 409 if (x > 0.6374)
408 410 yy = GT3(x - t0z3);
409 411 else
410 412 yy = GT2(x - t0z2);
411 - } else
413 + } else {
412 414 yy = GT1(x - t0z1);
415 + }
413 416
414 417 /* compute gamma(x+i) = (x+i-1)*...*(x+1)*yy, 0<i<8 */
415 418 switch (i) {
416 - case 0: /* yy/x */
419 + case 0: /* yy/x */
417 420 rr = yy / x;
418 421 break;
419 - case 1: /* yy */
422 + case 1: /* yy */
420 423 rr = yy;
421 424 break;
422 - case 2: /* (x+1)*yy */
425 + case 2: /* (x+1)*yy */
423 426 rr = (x + one) * yy;
424 427 break;
425 - case 3: /* (x+2)*(x+1)*yy */
428 + case 3: /* (x+2)*(x+1)*yy */
426 429 rr = (x + one) * (x + two) * yy;
427 430 break;
428 431
429 - case 4: /* (x+1)*(x+3)*(x+2)*yy */
432 + case 4: /* (x+1)*(x+3)*(x+2)*yy */
430 433 rr = (x + one) * (x + two) * ((x + 3.0) * yy);
431 434 break;
432 - case 5: /* ((x+1)*(x+4)*(x+2)*(x+3))*yy */
435 + case 5: /* ((x+1)*(x+4)*(x+2)*(x+3))*yy */
433 436 z1 = (x + two) * (x + 3.0) * yy;
434 437 z2 = (x + one) * (x + 4.0);
435 438 rr = z1 * z2;
436 439 break;
437 - case 6: /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */
440 + case 6: /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */
438 441 z1 = (x + two) * (x + 3.0);
439 442 z2 = (x + 5.0) * yy;
440 443 rr = z1 * (z1 - two) * z2;
441 444 break;
442 - case 7: /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */
445 + case 7: /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */
443 446 z1 = (x + two) * (x + 3.0);
444 447 z2 = (x + 5.0) * (x + 6.0) * yy;
445 448 rr = z1 * (z1 - two) * z2;
446 449 break;
447 450 }
451 +
448 452 return (rr);
449 453 }
450 454
451 455 float
452 -tgammaf(float xf) {
456 +tgammaf(float xf)
457 +{
453 458 float zf;
454 459 double ss, ww;
455 460 double x, y, z;
456 461 int i, j, k, ix, hx, xk;
457 462
458 - hx = *(int *) &xf;
463 + hx = *(int *)&xf;
459 464 ix = hx & 0x7fffffff;
460 465
461 - x = (double) xf;
466 + x = (double)xf;
467 +
462 468 if (ix < 0x33800000)
463 - return (1.0F / xf); /* |x| < 2**-24 */
469 + return (1.0F / xf); /* |x| < 2**-24 */
464 470
465 471 if (ix >= 0x7f800000)
466 - return (xf * ((hx < 0)? 0.0F : xf)); /* +-Inf or NaN */
472 + return (xf * ((hx < 0) ? 0.0F : xf)); /* +-Inf or NaN */
467 473
468 - if (hx > 0x420C290F) /* x > 35.040096283... overflow */
469 - return (float)(x / tiny);
474 + if (hx > 0x420C290F) /* x > 35.040096283... overflow */
475 + return ((float)(x / tiny));
470 476
471 - if (hx >= 0x41000000) /* x >= 8 */
472 - return ((float) large_gam(x));
477 + if (hx >= 0x41000000) /* x >= 8 */
478 + return ((float)large_gam(x));
473 479
474 - if (hx > 0) { /* 0 < x < 8 */
475 - i = (int) xf;
476 - return ((float) gam_n(i, x - (double) i));
480 + if (hx > 0) { /* 0 < x < 8 */
481 + i = (int)xf;
482 + return ((float)gam_n(i, x - (double)i));
477 483 }
478 484
479 - /* negative x */
480 - /* INDENT OFF */
485 + /*
486 + * negative x
487 + */
488 +
481 489 /*
482 490 * compute xk =
483 491 * -2 ... x is an even int (-inf is considered even)
484 492 * -1 ... x is an odd int
485 493 * +0 ... x is not an int but chopped to an even int
486 494 * +1 ... x is not an int but chopped to an odd int
487 495 */
488 - /* INDENT ON */
489 496 xk = 0;
497 +
490 498 if (ix >= 0x4b000000) {
491 499 if (ix > 0x4b000000)
492 500 xk = -2;
493 501 else
494 502 xk = -2 + (ix & 1);
495 503 } else if (ix >= 0x3f800000) {
496 504 k = (ix >> 23) - 0x7f;
497 505 j = ix >> (23 - k);
506 +
498 507 if ((j << (23 - k)) == ix)
499 508 xk = -2 + (j & 1);
500 509 else
501 510 xk = j & 1;
502 511 }
512 +
503 513 if (xk < 0) {
504 514 /* 0/0 invalid NaN, ideally gamma(-n)= (-1)**(n+1) * inf */
505 515 zf = xf - xf;
506 516 return (zf / zf);
507 517 }
508 518
509 519 /* negative underflow thresold */
510 - if (ix > 0x4224000B) { /* x < -(41+11ulp) */
520 + if (ix > 0x4224000B) { /* x < -(41+11ulp) */
511 521 if (xk == 0)
512 522 z = -tiny;
513 523 else
514 524 z = tiny;
525 +
515 526 return ((float)z);
516 527 }
517 528
518 - /* INDENT OFF */
519 - /* now compute gamma(x) by -1/((sin(pi*y)/pi)*gamma(1+y)), y = -x */
529 + /*
530 + * now compute gamma(x) by -1/((sin(pi*y)/pi)*gamma(1+y)), y = -x
531 + */
532 +
520 533 /*
521 534 * First compute ss = -sin(pi*y)/pi , so that
522 535 * gamma(x) = 1/(ss*gamma(1+y))
523 536 */
524 - /* INDENT ON */
525 537 y = -x;
526 - j = (int) y;
527 - z = y - (double) j;
528 - if (z > 0.3183098861837906715377675)
538 + j = (int)y;
539 + z = y - (double)j;
540 +
541 + if (z > 0.3183098861837906715377675) {
529 542 if (z > 0.6816901138162093284622325)
530 543 ss = kpsin(one - z);
531 544 else
532 545 ss = kpcos(0.5 - z);
533 - else
546 + } else {
534 547 ss = kpsin(z);
548 + }
549 +
535 550 if (xk == 0)
536 551 ss = -ss;
537 552
538 553 /* Then compute ww = gamma(1+y) */
539 554 if (j < 7)
540 555 ww = gam_n(j + 1, z);
541 556 else
542 557 ww = large_gam(y + one);
543 558
544 559 /* return 1/(ss*ww) */
545 - return ((float) (one / (ww * ss)));
560 + return ((float)(one / (ww * ss)));
546 561 }
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX