Print this page
| Split |
Close |
| Expand all |
| Collapse all |
--- old/usr/src/lib/libm/common/m9x/tgammal.c
+++ new/usr/src/lib/libm/common/m9x/tgammal.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.
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 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 27 * Use is subject to license terms.
28 28 */
29 29
30 30 #if defined(ELFOBJ)
31 31 #pragma weak tgammal = __tgammal
32 32 #endif
33 33
34 34 #include "libm.h"
35 35 #include <sys/isa_defs.h>
36 36
37 37 #if defined(_BIG_ENDIAN)
38 38 #define H0_WORD(x) ((unsigned *) &x)[0]
39 39 #define H3_WORD(x) ((unsigned *) &x)[3]
40 40 #define CHOPPED(x) (long double) ((double) (x))
41 41 #else
42 42 #define H0_WORD(x) ((((int *) &x)[2] << 16) | \
43 43 (0x0000ffff & (((unsigned *) &x)[1] >> 15)))
44 44 #define H3_WORD(x) ((unsigned *) &x)[0]
45 45 #define CHOPPED(x) (long double) ((float) (x))
46 46 #endif
47 47
48 48 struct LDouble {
49 49 long double h, l;
50 50 };
51 51
52 52 /* INDENT OFF */
53 53 /* Primary interval GTi() */
54 54 static const long double P1[] = {
55 55 +0.709086836199777919037185741507610124611513720557L,
56 56 +4.45754781206489035827915969367354835667391606951e-0001L,
57 57 +3.21049298735832382311662273882632210062918153852e-0002L,
58 58 -5.71296796342106617651765245858289197369688864350e-0003L,
59 59 +6.04666892891998977081619174969855831606965352773e-0003L,
60 60 +8.99106186996888711939627812174765258822658645168e-0004L,
61 61 -6.96496846144407741431207008527018441810175568949e-0005L,
62 62 +1.52597046118984020814225409300131445070213882429e-0005L,
63 63 +5.68521076168495673844711465407432189190681541547e-0007L,
64 64 +3.30749673519634895220582062520286565610418952979e-0008L,
65 65 };
66 66 static const long double Q1[] = {
67 67 +1.0+0000L,
68 68 +1.35806511721671070408570853537257079579490650668e+0000L,
69 69 +2.97567810153429553405327140096063086994072952961e-0001L,
70 70 -1.52956835982588571502954372821681851681118097870e-0001L,
71 71 -2.88248519561420109768781615289082053597954521218e-0002L,
72 72 +1.03475311719937405219789948456313936302378395955e-0002L,
73 73 +4.12310203243891222368965360124391297374822742313e-0004L,
74 74 -3.12653708152290867248931925120380729518332507388e-0004L,
75 75 +2.36672170850409745237358105667757760527014332458e-0005L,
76 76 };
77 77 static const long double P2[] = {
78 78 +0.428486815855585429730209907810650135255270600668084114L,
79 79 +2.62768479103809762805691743305424077975230551176e-0001L,
80 80 +3.81187532685392297608310837995193946591425896150e-0002L,
81 81 +3.00063075891811043820666846129131255948527925381e-0003L,
82 82 +2.47315407812279164228398470797498649142513408654e-0003L,
83 83 +3.62838199917848372586173483147214880464782938664e-0004L,
84 84 +3.43991105975492623982725644046473030098172692423e-0006L,
85 85 +4.56902151569603272237014240794257659159045432895e-0006L,
86 86 +2.13734755837595695602045100675540011352948958453e-0007L,
87 87 +9.74123440547918230781670266967882492234877125358e-0009L,
88 88 };
89 89 static const long double Q2[] = {
90 90 +1.0L,
91 91 +9.18284118632506842664645516830761489700556179701e-0001L,
92 92 -6.41430858837830766045202076965923776189154874947e-0003L,
93 93 -1.24400885809771073213345747437964149775410921376e-0001L,
94 94 +4.69803798146251757538856567522481979624746875964e-0003L,
95 95 +7.18309447069495315914284705109868696262662082731e-0003L,
96 96 -8.75812626987894695112722600697653425786166399105e-0004L,
97 97 -1.23539972377769277995959339188431498626674835169e-0004L,
98 98 +3.10019017590151598732360097849672925448587547746e-0005L,
99 99 -1.77260223349332617658921874288026777465782364070e-0006L,
100 100 };
101 101 static const long double P3[] = {
102 102 +0.3824094797345675048502747661075355640070439388902L,
103 103 +3.42198093076618495415854906335908427159833377774e-0001L,
104 104 +9.63828189500585568303961406863153237440702754858e-0002L,
105 105 +8.76069421042696384852462044188520252156846768667e-0003L,
106 106 +1.86477890389161491224872014149309015261897537488e-0003L,
107 107 +8.16871354540309895879974742853701311541286944191e-0004L,
108 108 +6.83783483674600322518695090864659381650125625216e-0005L,
109 109 -1.10168269719261574708565935172719209272190828456e-0006L,
110 110 +9.66243228508380420159234853278906717065629721016e-0007L,
111 111 +2.31858885579177250541163820671121664974334728142e-0008L,
112 112 };
113 113 static const long double Q3[] = {
114 114 +1.0L,
115 115 +8.25479821168813634632437430090376252512793067339e-0001L,
116 116 -1.62251363073937769739639623669295110346015576320e-0002L,
117 117 -1.10621286905916732758745130629426559691187579852e-0001L,
118 118 +3.48309693970985612644446415789230015515365291459e-0003L,
119 119 +6.73553737487488333032431261131289672347043401328e-0003L,
120 120 -7.63222008393372630162743587811004613050245128051e-0004L,
121 121 -1.35792670669190631476784768961953711773073251336e-0004L,
122 122 +3.19610150954223587006220730065608156460205690618e-0005L,
123 123 -1.82096553862822346610109522015129585693354348322e-0006L,
124 124 };
125 125
126 126 static const long double
127 127 #if defined(__x86)
128 128 GZ1_h = 0.938204627909682449364570100414084663498215377L,
129 129 GZ1_l = 4.518346116624229420055327632718530617227944106e-20L,
130 130 GZ2_h = 0.885603194410888700264725126309883762587560340L,
131 131 GZ2_l = 1.409077427270497062039119290776508217077297169e-20L,
132 132 GZ3_h = 0.936781411463652321613537060640553022494714241L,
133 133 GZ3_l = 5.309836440284827247897772963887219035221996813e-21L,
134 134 #else
135 135 GZ1_h = 0.938204627909682449409753561580326910854647031L,
136 136 GZ1_l = 4.684412162199460089642452580902345976446297037e-35L,
137 137 GZ2_h = 0.885603194410888700278815900582588658192658794L,
138 138 GZ2_l = 7.501529273890253789219935569758713534641074860e-35L,
139 139 GZ3_h = 0.936781411463652321618846897080837818855399840L,
140 140 GZ3_l = 3.088721217404784363585591914529361687403776917e-35L,
141 141 #endif
142 142 TZ1 = -0.3517214357852935791015625L,
143 143 TZ3 = 0.280530631542205810546875L;
144 144 /* INDENT ON */
145 145
146 146 /* INDENT OFF */
147 147 /*
148 148 * compute gamma(y=yh+yl) for y in GT1 = [1.0000, 1.2845]
149 149 * ...assume yh got 53 or 24(i386) significant bits
150 150 */
151 151 /* INDENT ON */
152 152 static struct LDouble
153 153 GT1(long double yh, long double yl) {
154 154 long double t3, t4, y;
155 155 int i;
156 156 struct LDouble r;
157 157
158 158 y = yh + yl;
159 159 for (t4 = Q1[8], t3 = P1[8] + y * P1[9], i = 7; i >= 0; i--) {
160 160 t4 = t4 * y + Q1[i];
161 161 t3 = t3 * y + P1[i];
162 162 }
163 163 t3 = (y * y) * t3 / t4;
164 164 t3 += (TZ1 * yl + GZ1_l);
165 165 t4 = TZ1 * yh;
166 166 r.h = CHOPPED((t4 + GZ1_h + t3));
167 167 t3 += (t4 - (r.h - GZ1_h));
168 168 r.l = t3;
169 169 return (r);
170 170 }
171 171
172 172 /* INDENT OFF */
173 173 /*
174 174 * compute gamma(y=yh+yl) for y in GT2 = [1.2844, 1.6374]
175 175 * ...assume yh got 53 significant bits
176 176 */
177 177 /* INDENT ON */
178 178 static struct LDouble
179 179 GT2(long double yh, long double yl) {
180 180 long double t3, t4, y;
181 181 int i;
182 182 struct LDouble r;
183 183
184 184 y = yh + yl;
185 185 for (t4 = Q2[9], t3 = P2[9], i = 8; i >= 0; i--) {
186 186 t4 = t4 * y + Q2[i];
187 187 t3 = t3 * y + P2[i];
188 188 }
189 189 t3 = GZ2_l + (y * y) * t3 / t4;
190 190 r.h = CHOPPED((GZ2_h + t3));
191 191 r.l = t3 - (r.h - GZ2_h);
192 192 return (r);
193 193 }
194 194
195 195 /* INDENT OFF */
196 196 /*
197 197 * compute gamma(y=yh+yl) for y in GT3 = [1.6373, 2.0000]
198 198 * ...assume yh got 53 significant bits
199 199 */
200 200 /* INDENT ON */
201 201 static struct LDouble
202 202 GT3(long double yh, long double yl) {
203 203 long double t3, t4, y;
204 204 int i;
205 205 struct LDouble r;
206 206
207 207 y = yh + yl;
208 208 for (t4 = Q3[9], t3 = P3[9], i = 8; i >= 0; i--) {
209 209 t4 = t4 * y + Q3[i];
210 210 t3 = t3 * y + P3[i];
211 211 }
212 212 t3 = (y * y) * t3 / t4;
213 213 t3 += (TZ3 * yl + GZ3_l);
214 214 t4 = TZ3 * yh;
215 215 r.h = CHOPPED((t4 + GZ3_h + t3));
216 216 t3 += (t4 - (r.h - GZ3_h));
217 217 r.l = t3;
218 218 return (r);
219 219 }
220 220
221 221 /* INDENT OFF */
222 222 /* Hex value of GP[0] shoule be 3FB55555 55555555 */
223 223 static const long double GP[] = {
224 224 +0.083333333333333333333333333333333172839171301L,
225 225 -2.77777777777777777777777777492501211999399424104e-0003L,
226 226 +7.93650793650793650793635650541638236350020883243e-0004L,
227 227 -5.95238095238095238057299772679324503339241961704e-0004L,
228 228 +8.41750841750841696138422987977683524926142600321e-0004L,
229 229 -1.91752691752686682825032547823699662178842123308e-0003L,
230 230 +6.41025641022403480921891559356473451161279359322e-0003L,
231 231 -2.95506535798414019189819587455577003732808185071e-0002L,
232 232 +1.79644367229970031486079180060923073476568732136e-0001L,
233 233 -1.39243086487274662174562872567057200255649290646e+0000L,
234 234 +1.34025874044417962188677816477842265259608269775e+0001L,
235 235 -1.56803713480127469414495545399982508700748274318e+0002L,
236 236 +2.18739841656201561694927630335099313968924493891e+0003L,
237 237 -3.55249848644100338419187038090925410976237921269e+0004L,
238 238 +6.43464880437835286216768959439484376449179576452e+0005L,
239 239 -1.20459154385577014992600342782821389605893904624e+0007L,
240 240 +2.09263249637351298563934942349749718491071093210e+0008L,
241 241 -2.96247483183169219343745316433899599834685703457e+0009L,
242 242 +2.88984933605896033154727626086506756972327292981e+0010L,
243 243 -1.40960434146030007732838382416230610302678063984e+0011L, /* 19 */
244 244 };
245 245
246 246 static const long double T3[] = {
247 247 +0.666666666666666666666666666666666634567834260213L, /* T3[0] */
248 248 +0.400000000000000000000000000040853636176634934140L, /* T3[1] */
249 249 +0.285714285714285714285696975252753987869020263448L, /* T3[2] */
250 250 +0.222222222222222225593221101192317258554772129875L, /* T3[3] */
251 251 +0.181818181817850192105847183461778186703779262916L, /* T3[4] */
252 252 +0.153846169861348633757101285952333369222567014596L, /* T3[5] */
253 253 +0.133033462889260193922261296772841229985047571265L, /* T3[6] */
254 254 };
255 255
256 256 static const long double c[] = {
257 257 0.0L,
258 258 1.0L,
259 259 2.0L,
260 260 0.5L,
261 261 1.0e-4930L, /* tiny */
262 262 4.18937683105468750000e-01L, /* hln2pim1_h */
263 263 8.50099203991780329736405617639861397473637783412817152e-07L, /* hln2pim1_l */
264 264 0.418938533204672741780329736405617639861397473637783412817152L, /* hln2pim1 */
265 265 2.16608493865351192653179168701171875e-02L, /* ln2_32hi */
266 266 5.96317165397058692545083025235937919875797669127130e-12L, /* ln2_32lo */
267 267 46.16624130844682903551758979206054839765267053289554989233L, /* invln2_32 */
268 268 #if defined(__x86)
269 269 1.7555483429044629170023839037639845628291e+03L, /* overflow */
270 270 #else
271 271 1.7555483429044629170038892160702032034177e+03L, /* overflow */
272 272 #endif
273 273 };
274 274
275 275 #define zero c[0]
276 276 #define one c[1]
277 277 #define two c[2]
278 278 #define half c[3]
279 279 #define tiny c[4]
280 280 #define hln2pim1_h c[5]
281 281 #define hln2pim1_l c[6]
282 282 #define hln2pim1 c[7]
283 283 #define ln2_32hi c[8]
284 284 #define ln2_32lo c[9]
285 285 #define invln2_32 c[10]
286 286 #define overflow c[11]
287 287
288 288 /*
289 289 * |exp(r) - (1+r+Et0*r^2+...+Et10*r^12)| <= 2^(-128.88) for |r|<=ln2/64
290 290 */
291 291 static const long double Et[] = {
292 292 +5.0000000000000000000e-1L,
293 293 +1.66666666666666666666666666666828835166292152466e-0001L,
294 294 +4.16666666666666666666666666666693398646592712189e-0002L,
295 295 +8.33333333333333333333331748774512601775591115951e-0003L,
296 296 +1.38888888888888888888888845356011511394764753997e-0003L,
297 297 +1.98412698412698413237140350092993252684198882102e-0004L,
298 298 +2.48015873015873016080222025357442659895814371694e-0005L,
299 299 +2.75573192239028921114572986441972140933432317798e-0006L,
300 300 +2.75573192239448470555548102895526369739856219317e-0007L,
301 301 +2.50521677867683935940853997995937600214167232477e-0008L,
302 302 +2.08767928899010367374984448513685566514152147362e-0009L,
303 303 };
304 304
305 305 /*
306 306 * long double precision coefficients for computing log(x)-1 in tgamma.
307 307 * See "algorithm" for details
308 308 *
309 309 * log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y, 1<=y<2,
310 310 * j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
311 311 * T1(n) = T1[2n,2n+1] = n*log(2)-1,
312 312 * T2(j) = T2[2j,2j+1] = log(z[j]),
313 313 * T3(s) = 2s + T3[0]s^3 + T3[1]s^5 + T3[2]s^7 + ... + T3[6]s^15
314 314 * Note
315 315 * (1) the leading entries are truncated to 24 binary point.
316 316 * (2) Remez error for T3(s) is bounded by 2**(-136.54)
317 317 */
318 318 static const long double T1[] = {
319 319 -1.000000000000000000000000000000000000000000e+00L,
320 320 +0.000000000000000000000000000000000000000000e+00L,
321 321 -3.068528175354003906250000000000000000000000e-01L,
322 322 -1.904654299957767878541823431924500011926579e-09L,
323 323 +3.862943053245544433593750000000000000000000e-01L,
324 324 +5.579533617547508924291635313615100141107647e-08L,
325 325 +1.079441487789154052734375000000000000000000e+00L,
326 326 +5.389068187551732136437452970422650211661470e-08L,
327 327 +1.772588670253753662109375000000000000000000e+00L,
328 328 +5.198602757555955348583270627230200282215294e-08L,
329 329 +2.465735852718353271484375000000000000000000e+00L,
330 330 +5.008137327560178560729088284037750352769117e-08L,
331 331 +3.158883035182952880859375000000000000000000e+00L,
332 332 +4.817671897564401772874905940845299849351090e-08L,
333 333 +3.852030217647552490234375000000000000000000e+00L,
334 334 +4.627206467568624985020723597652849919904913e-08L,
335 335 +4.545177400112152099609375000000000000000000e+00L,
336 336 +4.436741037572848197166541254460399990458737e-08L,
337 337 +5.238324582576751708984375000000000000000000e+00L,
338 338 +4.246275607577071409312358911267950061012560e-08L,
339 339 +5.931471765041351318359375000000000000000000e+00L,
340 340 +4.055810177581294621458176568075500131566384e-08L,
341 341 };
342 342
343 343 /*
344 344 * T2[2i,2i+1] = log(1+i/64+1/128)
345 345 */
346 346 static const long double T2[] = {
347 347 +7.7821016311645507812500000000000000000000e-03L,
348 348 +3.8810890398166212900061136763678127453570e-08L,
349 349 +2.3167014122009277343750000000000000000000e-02L,
350 350 +4.5159525100885049160962289916579411752759e-08L,
351 351 +3.8318812847137451171875000000000000000000e-02L,
352 352 +5.1454999148021880325123797290345960518164e-08L,
353 353 +5.3244471549987792968750000000000000000000e-02L,
354 354 +4.2968824489897120193786528776939573415076e-08L,
355 355 +6.7950606346130371093750000000000000000000e-02L,
356 356 +5.5562377378300815277772629414034632394030e-08L,
357 357 +8.2443654537200927734375000000000000000000e-02L,
358 358 +1.4673873663533785068668307805914095366600e-08L,
359 359 +9.6729576587677001953125000000000000000000e-02L,
360 360 +4.9870874110342446056487463437015041543346e-08L,
361 361 +1.1081433296203613281250000000000000000000e-01L,
362 362 +3.3378253981382306169323211928098474801099e-08L,
363 363 +1.2470346689224243164062500000000000000000e-01L,
364 364 +1.1608714804222781515380863268491613205318e-08L,
365 365 +1.3840228319168090820312500000000000000000e-01L,
366 366 +3.9667438227482200873601649187393160823607e-08L,
367 367 +1.5191602706909179687500000000000000000000e-01L,
368 368 +1.4956750178196803424896884511327584958252e-08L,
369 369 +1.6524952650070190429687500000000000000000e-01L,
370 370 +4.6394605258578736449277240313729237989366e-08L,
371 371 +1.7840760946273803710937500000000000000000e-01L,
372 372 +4.8010080260010025241510941968354682199540e-08L,
373 373 +1.9139480590820312500000000000000000000000e-01L,
374 374 +4.7091426329609298807561308873447039132856e-08L,
375 375 +2.0421552658081054687500000000000000000000e-01L,
376 376 +1.4847880344628820386196239272213742113867e-08L,
377 377 +2.1687388420104980468750000000000000000000e-01L,
378 378 +5.4099564554931589525744347498478964801484e-08L,
379 379 +2.2937405109405517578125000000000000000000e-01L,
380 380 +4.9970790654210230725046139871550961365282e-08L,
381 381 +2.4171990156173706054687500000000000000000e-01L,
382 382 +3.5325408107597432515913513900103385655073e-08L,
383 383 +2.5391519069671630859375000000000000000000e-01L,
384 384 +1.9284247135543573297906606667466299224747e-08L,
385 385 +2.6596349477767944335937500000000000000000e-01L,
386 386 +5.3719458497979750926537543389268821141517e-08L,
387 387 +2.7786844968795776367187500000000000000000e-01L,
388 388 +1.3154985425144750329234012330820349974537e-09L,
389 389 +2.8963327407836914062500000000000000000000e-01L,
390 390 +1.8504673536253893055525668970003860369760e-08L,
391 391 +3.0126130580902099609375000000000000000000e-01L,
392 392 +2.4769140784919125538233755492657352680723e-08L,
393 393 +3.1275570392608642578125000000000000000000e-01L,
394 394 +6.0778104626049965596883190321597861455475e-09L,
395 395 +3.2411944866180419921875000000000000000000e-01L,
396 396 +1.9992407776871920760434987352182336158873e-08L,
397 397 +3.3535552024841308593750000000000000000000e-01L,
398 398 +2.1672724744319679579814166199074433006807e-08L,
399 399 +3.4646672010421752929687500000000000000000e-01L,
400 400 +4.7241991051621587188425772950711830538414e-08L,
401 401 +3.5745584964752197265625000000000000000000e-01L,
402 402 +3.9274281801569759490140904474434669956562e-08L,
403 403 +3.6832553148269653320312500000000000000000e-01L,
404 404 +2.9676011119845105154050398826897178765758e-08L,
405 405 +3.7907832860946655273437500000000000000000e-01L,
406 406 +2.4325502905656478345631019858881408009210e-08L,
407 407 +3.8971674442291259765625000000000000000000e-01L,
408 408 +6.7171126157142136040035208670510556529487e-09L,
409 409 +4.0024316310882568359375000000000000000000e-01L,
410 410 +1.0181870233355751019951311700799406124957e-09L,
411 411 +4.1065990924835205078125000000000000000000e-01L,
412 412 +1.5736916335153056203175822787661567534220e-08L,
413 413 +4.2096924781799316406250000000000000000000e-01L,
414 414 +4.6826136472066367161506795972449857268707e-08L,
415 415 +4.3117344379425048828125000000000000000000e-01L,
416 416 +2.1024120852577922478955594998480144051225e-08L,
417 417 +4.4127452373504638671875000000000000000000e-01L,
418 418 +3.7069828842770746441661301225362605528786e-08L,
419 419 +4.5127463340759277343750000000000000000000e-01L,
420 420 +1.0731865811707192383079012478685922879010e-08L,
421 421 +4.6117568016052246093750000000000000000000e-01L,
422 422 +3.4961647705430499925597855358603099030515e-08L,
423 423 +4.7097969055175781250000000000000000000000e-01L,
424 424 +2.4667033200046897856056359251373510964634e-08L,
425 425 +4.8068851232528686523437500000000000000000e-01L,
426 426 +1.7020465042442243455448011551208861216878e-08L,
427 427 +4.9030393362045288085937500000000000000000e-01L,
428 428 +5.4424740957290971159645746860530583309571e-08L,
429 429 +4.9982786178588867187500000000000000000000e-01L,
430 430 +7.7705606579463314152470441415126573566105e-09L,
431 431 +5.0926184654235839843750000000000000000000e-01L,
432 432 +5.5247449548366574919228323824878565745713e-08L,
433 433 +5.1860773563385009765625000000000000000000e-01L,
434 434 +2.8574195534496726996364798698556235730848e-08L,
435 435 +5.2786707878112792968750000000000000000000e-01L,
436 436 +1.0839714455426392217778300963558522088193e-08L,
437 437 +5.3704142570495605468750000000000000000000e-01L,
438 438 +4.0191927599879229244153832299023744345999e-08L,
439 439 +5.4613238573074340820312500000000000000000e-01L,
440 440 +5.1867392242179272209231209163864971792889e-08L,
441 441 +5.5514144897460937500000000000000000000000e-01L,
442 442 +5.8565892217715480359515904050170125743178e-08L,
443 443 +5.6407010555267333984375000000000000000000e-01L,
444 444 +3.2732129626227634290090190711817681692354e-08L,
445 445 +5.7291972637176513671875000000000000000000e-01L,
446 446 +2.7190020372374006726626261068626400393936e-08L,
447 447 +5.8169168233871459960937500000000000000000e-01L,
448 448 +5.7295907882911235753725372340709967597394e-08L,
449 449 +5.9038740396499633789062500000000000000000e-01L,
450 450 +4.2637180036751291708123598757577783615014e-08L,
451 451 +5.9900814294815063476562500000000000000000e-01L,
452 452 +4.6697932764615975024461651502060474048774e-08L,
453 453 +6.0755521059036254882812500000000000000000e-01L,
454 454 +3.9634179246672960152791125371893149820625e-08L,
455 455 +6.1602985858917236328125000000000000000000e-01L,
456 456 +1.8626341656366315928196700650292529688219e-08L,
457 457 +6.2443327903747558593750000000000000000000e-01L,
458 458 +8.9744179151050387440546731199093039879228e-09L,
459 459 +6.3276666402816772460937500000000000000000e-01L,
460 460 +5.5428701049364114685035797584887586099726e-09L,
461 461 +6.4103114604949951171875000000000000000000e-01L,
462 462 +3.3371431779336851334405392546708949047361e-08L,
463 463 +6.4922791719436645507812500000000000000000e-01L,
464 464 +2.9430743363812714969905311122271269100885e-08L,
465 465 +6.5735805034637451171875000000000000000000e-01L,
466 466 +2.2361985518423140023245936165514147093250e-08L,
467 467 +6.6542261838912963867187500000000000000000e-01L,
468 468 +1.4155960810278217610006660181148303091649e-08L,
469 469 +6.7342263460159301757812500000000000000000e-01L,
470 470 +4.0610573702719835388801017264750843477878e-08L,
471 471 +6.8135917186737060546875000000000000000000e-01L,
472 472 +5.2940532463479321559568089441735584156689e-08L,
473 473 +6.8923324346542358398437500000000000000000e-01L,
474 474 +3.7773385396340539337814603903232796216537e-08L,
475 475 };
476 476
477 477 /*
478 478 * S[j],S_trail[j] = 2**(j/32.) for the final computation of exp(t+w)
479 479 */
480 480 static const long double S[] = {
481 481 #if defined(__x86)
482 482 +1.0000000000000000000000000e+00L,
483 483 +1.0218971486541166782081522e+00L,
484 484 +1.0442737824274138402382006e+00L,
485 485 +1.0671404006768236181297224e+00L,
486 486 +1.0905077326652576591003302e+00L,
487 487 +1.1143867425958925362894369e+00L,
488 488 +1.1387886347566916536971221e+00L,
489 489 +1.1637248587775775137938619e+00L,
490 490 +1.1892071150027210666875674e+00L,
491 491 +1.2152473599804688780476325e+00L,
492 492 +1.2418578120734840485256747e+00L,
493 493 +1.2690509571917332224885722e+00L,
494 494 +1.2968395546510096659215822e+00L,
495 495 +1.3252366431597412945939118e+00L,
496 496 +1.3542555469368927282668852e+00L,
497 497 +1.3839098819638319548151403e+00L,
498 498 +1.4142135623730950487637881e+00L,
499 499 +1.4451808069770466200253470e+00L,
500 500 +1.4768261459394993113155431e+00L,
501 501 +1.5091644275934227397133885e+00L,
502 502 +1.5422108254079408235859630e+00L,
503 503 +1.5759808451078864864006862e+00L,
504 504 +1.6104903319492543080837174e+00L,
505 505 +1.6457554781539648445110730e+00L,
506 506 +1.6817928305074290860378350e+00L,
507 507 +1.7186192981224779156032914e+00L,
508 508 +1.7562521603732994831094730e+00L,
509 509 +1.7947090750031071864148413e+00L,
510 510 +1.8340080864093424633989166e+00L,
511 511 +1.8741676341102999013002103e+00L,
512 512 +1.9152065613971472938202589e+00L,
513 513 +1.9571441241754002689657438e+00L,
514 514 #else
515 515 +1.00000000000000000000000000000000000e+00L,
516 516 +1.02189714865411667823448013478329942e+00L,
517 517 +1.04427378242741384032196647873992910e+00L,
518 518 +1.06714040067682361816952112099280918e+00L,
519 519 +1.09050773266525765920701065576070789e+00L,
520 520 +1.11438674259589253630881295691960313e+00L,
521 521 +1.13878863475669165370383028384151134e+00L,
522 522 +1.16372485877757751381357359909218536e+00L,
523 523 +1.18920711500272106671749997056047593e+00L,
524 524 +1.21524735998046887811652025133879836e+00L,
525 525 +1.24185781207348404859367746872659561e+00L,
526 526 +1.26905095719173322255441908103233805e+00L,
527 527 +1.29683955465100966593375411779245118e+00L,
528 528 +1.32523664315974129462953709549872168e+00L,
529 529 +1.35425554693689272829801474014070273e+00L,
530 530 +1.38390988196383195487265952726519287e+00L,
531 531 +1.41421356237309504880168872420969798e+00L,
532 532 +1.44518080697704662003700624147167095e+00L,
533 533 +1.47682614593949931138690748037404985e+00L,
534 534 +1.50916442759342273976601955103319352e+00L,
535 535 +1.54221082540794082361229186209073479e+00L,
536 536 +1.57598084510788648645527016018190504e+00L,
537 537 +1.61049033194925430817952066735740067e+00L,
538 538 +1.64575547815396484451875672472582254e+00L,
539 539 +1.68179283050742908606225095246642969e+00L,
540 540 +1.71861929812247791562934437645631244e+00L,
541 541 +1.75625216037329948311216061937531314e+00L,
542 542 +1.79470907500310718642770324212778174e+00L,
543 543 +1.83400808640934246348708318958828892e+00L,
544 544 +1.87416763411029990132999894995444645e+00L,
545 545 +1.91520656139714729387261127029583086e+00L,
546 546 +1.95714412417540026901832225162687149e+00L,
547 547 #endif
548 548 };
549 549 static const long double S_trail[] = {
550 550 #if defined(__x86)
551 551 +0.0000000000000000000000000e+00L,
552 552 +2.6327965667180882569382524e-20L,
553 553 +8.3765863521895191129661899e-20L,
554 554 +3.9798705777454504249209575e-20L,
555 555 +1.0668046596651558640993042e-19L,
556 556 +1.9376009847285360448117114e-20L,
557 557 +6.7081819456112953751277576e-21L,
558 558 +1.9711680502629186462729727e-20L,
559 559 +2.9932584438449523689104569e-20L,
560 560 +6.8887754153039109411061914e-20L,
561 561 +6.8002718741225378942847820e-20L,
562 562 +6.5846917376975403439742349e-20L,
563 563 +1.2171958727511372194876001e-20L,
564 564 +3.5625253228704087115438260e-20L,
565 565 +3.1129551559077560956309179e-20L,
566 566 +5.7519192396164779846216492e-20L,
567 567 +3.7900651177865141593101239e-20L,
568 568 +1.1659262405698741798080115e-20L,
569 569 +7.1364385105284695967172478e-20L,
570 570 +5.2631003710812203588788949e-20L,
571 571 +2.6328853788732632868460580e-20L,
572 572 +5.4583950085438242788190141e-20L,
573 573 +9.5803254376938269960718656e-20L,
574 574 +7.6837733983874245823512279e-21L,
575 575 +2.4415965910835093824202087e-20L,
576 576 +2.6052966871016580981769728e-20L,
577 577 +2.6876456344632553875309579e-21L,
578 578 +1.2861930155613700201703279e-20L,
579 579 +8.8166633394037485606572294e-20L,
580 580 +2.9788615389580190940837037e-20L,
581 581 +5.2352341619805098677422139e-20L,
582 582 +5.2578463064010463732242363e-20L,
583 583 #else
584 584 +0.00000000000000000000000000000000000e+00L,
585 585 +1.80506787420330954745573333054573786e-35L,
586 586 -9.37452029228042742195756741973083214e-35L,
587 587 -1.59696844729275877071290963023149997e-35L,
588 588 +9.11249341012502297851168610167248666e-35L,
589 589 -6.50422820697854828723037477525938871e-35L,
590 590 -8.14846884452585113732569176748815532e-35L,
591 591 -5.06621457672180031337233074514290335e-35L,
592 592 -1.35983097468881697374987563824591912e-35L,
593 593 +9.49742763556319647030771056643324660e-35L,
594 594 -3.28317052317699860161506596533391526e-36L,
595 595 -5.01723570938719041029018653045842895e-35L,
596 596 -2.39147479768910917162283430160264014e-35L,
597 597 -8.35057135763390881529889073794408385e-36L,
598 598 +7.03675688907326504242173719067187644e-35L,
599 599 -5.18248485306464645753689301856695619e-35L,
600 600 +9.42224254862183206569211673639406488e-35L,
601 601 -3.96750082539886230916730613021641828e-35L,
602 602 +7.14352899156330061452327361509276724e-35L,
603 603 +1.15987125286798512424651783410044433e-35L,
604 604 +4.69693347835811549530973921320187447e-35L,
605 605 -3.38651317599500471079924198499981917e-35L,
606 606 -8.58731877429824706886865593510387445e-35L,
607 607 -9.60595154874935050318549936224606909e-35L,
608 608 +9.60973393212801278450755869714178581e-35L,
609 609 +6.37839792144002843924476144978084855e-35L,
610 610 +7.79243078569586424945646112516927770e-35L,
611 611 +7.36133776758845652413193083663393220e-35L,
612 612 -6.47299514791334723003521457561217053e-35L,
613 613 +8.58747441795369869427879806229522962e-35L,
614 614 +2.37181542282517483569165122830269098e-35L,
615 615 -3.02689168209611877300459737342190031e-37L,
616 616 #endif
617 617 };
618 618 /* INDENT ON */
619 619
620 620 /* INDENT OFF */
621 621 /*
622 622 * return tgamma(x) scaled by 2**-m for 8<x<=171.62... using Stirling's formula
623 623 * log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + (1/x)*P(1/(x*x))
624 624 * = L1 + L2 + L3,
625 625 */
626 626 /* INDENT ON */
627 627 static struct LDouble
628 628 large_gam(long double x, int *m) {
629 629 long double z, t1, t2, t3, z2, t5, w, y, u, r, v;
630 630 long double t24 = 16777216.0L, p24 = 1.0L / 16777216.0L;
631 631 int n2, j2, k, ix, j, i;
632 632 struct LDouble zz;
633 633 long double u2, ss_h, ss_l, r_h, w_h, w_l, t4;
634 634
635 635 /* INDENT OFF */
636 636 /*
637 637 * compute ss = ss.h+ss.l = log(x)-1 (see tgamma_log.h for details)
638 638 *
639 639 * log(x) - 1 = T1(n) + T2(j) + T3(s), where x = 2**n * y, 1<=y<2,
640 640 * j=[64*y], z[j]=1+j/64+1/128, s = (y-z[j])/(y+z[j]), and
641 641 * T1(n) = T1[2n,2n+1] = n*log(2)-1,
642 642 * T2(j) = T2[2j,2j+1] = log(z[j]),
643 643 * T3(s) = 2s + T3[0]s^3 + T3[1]s^5 + ... + T3[6]s^15
644 644 * Note
645 645 * (1) the leading entries are truncated to 24 binary point.
646 646 * (2) Remez error for T3(s) is bounded by 2**(-72.4)
647 647 * 2**(-24)
648 648 * _________V___________________
649 649 * T1(n): |_________|___________________|
650 650 * _______ ______________________
651 651 * T2(j): |_______|______________________|
652 652 * ____ _______________________
653 653 * 2s: |____|_______________________|
654 654 * __________________________
655 655 * + T3(s)-2s: |__________________________|
656 656 * -------------------------------------------
657 657 * [leading] + [Trailing]
658 658 */
659 659 /* INDENT ON */
660 660 ix = H0_WORD(x);
661 661 n2 = (ix >> 16) - 0x3fff; /* exponent of x, range:3-10 */
662 662 y = scalbnl(x, -n2); /* y = scale x to [1,2] */
663 663 n2 += n2; /* 2n */
664 664 j = (ix >> 10) & 0x3f; /* j */
665 665 z = 1.0078125L + (long double) j * 0.015625L; /* z[j]=1+j/64+1/128 */
666 666 j2 = j + j;
667 667 t1 = y + z;
668 668 t2 = y - z;
669 669 r = one / t1;
670 670 u = r * t2; /* u = (y-z)/(y+z) */
671 671 t1 = CHOPPED(t1);
672 672 t4 = T2[j2 + 1] + T1[n2 + 1];
673 673 z2 = u * u;
674 674 k = H0_WORD(u) & 0x7fffffff;
675 675 t3 = T2[j2] + T1[n2];
676 676 for (t5 = T3[6], i = 5; i >= 0; i--)
677 677 t5 = z2 * t5 + T3[i];
678 678 if ((k >> 16) < 0x3fec) { /* |u|<2**-19 */
679 679 t2 = t4 + u * (two + z2 * t5);
680 680 } else {
681 681 t5 = t4 + (u * z2) * t5;
682 682 u2 = u + u;
683 683 v = (long double) ((int) (u2 * t24)) * p24;
684 684 t2 = t5 + r * ((two * t2 - v * t1) - v * (y - (t1 - z)));
685 685 t3 += v;
686 686 }
687 687 ss_h = CHOPPED((t2 + t3));
688 688 ss_l = t2 - (ss_h - t3);
689 689 /* INDENT OFF */
690 690 /*
691 691 * compute ww = (x-.5)*(log(x)-1) + .5*(log(2pi)-1) + 1/x*(P(1/x^2)))
692 692 * where ss = log(x) - 1 in already in extra precision
693 693 */
694 694 /* INDENT ON */
695 695 z = one / x;
696 696 r = x - half;
697 697 r_h = CHOPPED((r));
698 698 w_h = r_h * ss_h + hln2pim1_h;
699 699 z2 = z * z;
700 700 w = (r - r_h) * ss_h + r * ss_l;
701 701 t1 = GP[19];
702 702 for (i = 18; i > 0; i--)
703 703 t1 = z2 * t1 + GP[i];
704 704 w += hln2pim1_l;
705 705 w_l = z * (GP[0] + z2 * t1) + w;
706 706 k = (int) ((w_h + w_l) * invln2_32 + half);
707 707
708 708 /* compute the exponential of w_h+w_l */
709 709
710 710 j = k & 0x1f;
711 711 *m = k >> 5;
712 712 t3 = (long double) k;
713 713
714 714 /* perform w - k*ln2_32 (represent as w_h - w_l) */
715 715 t1 = w_h - t3 * ln2_32hi;
716 716 t2 = t3 * ln2_32lo;
717 717 w = t2 - w_l;
718 718 w_h = t1 - w;
719 719 w_l = w - (t1 - w_h);
720 720
721 721 /* compute exp(w_h-w_l) */
722 722 z = w_h - w_l;
723 723 for (t1 = Et[10], i = 9; i >= 0; i--)
724 724 t1 = z * t1 + Et[i];
725 725 t3 = w_h - (w_l - (z * z) * t1); /* t3 = expm1(z) */
726 726 zz.l = S_trail[j] * (one + t3) + S[j] * t3;
727 727 zz.h = S[j];
728 728 return (zz);
729 729 }
730 730
731 731 /* INDENT OFF */
732 732 /*
733 733 * kpsin(x)= sin(pi*x)/pi
734 734 * 3 5 7 9 11 27
735 735 * = x+ks[0]*x +ks[1]*x +ks[2]*x +ks[3]*x +ks[4]*x + ... + ks[12]*x
736 736 */
737 737 static const long double ks[] = {
738 738 -1.64493406684822643647241516664602518705158902870e+0000L,
739 739 +8.11742425283353643637002772405874238094995726160e-0001L,
740 740 -1.90751824122084213696472111835337366232282723933e-0001L,
741 741 +2.61478478176548005046532613563241288115395517084e-0002L,
742 742 -2.34608103545582363750893072647117829448016479971e-0003L,
743 743 +1.48428793031071003684606647212534027556262040158e-0004L,
744 744 -6.97587366165638046518462722252768122615952898698e-0006L,
745 745 +2.53121740413702536928659271747187500934840057929e-0007L,
746 746 -7.30471182221385990397683641695766121301933621956e-0009L,
747 747 +1.71653847451163495739958249695549313987973589884e-0010L,
748 748 -3.34813314714560776122245796929054813458341420565e-0012L,
749 749 +5.50724992262622033449487808306969135431411753047e-0014L,
750 750 -7.67678132753577998601234393215802221104236979928e-0016L,
751 751 };
752 752 /* INDENT ON */
753 753
754 754 /*
755 755 * assume x is not tiny and positive
756 756 */
757 757 static struct LDouble
758 758 kpsin(long double x) {
759 759 long double z, t1, t2;
760 760 struct LDouble xx;
761 761 int i;
762 762
763 763 z = x * x;
764 764 xx.h = x;
765 765 for (t2 = ks[12], i = 11; i > 0; i--)
766 766 t2 = z * t2 + ks[i];
767 767 t1 = z * x;
768 768 t2 *= z * t1;
769 769 xx.l = t1 * ks[0] + t2;
770 770 return (xx);
771 771 }
772 772
773 773 /* INDENT OFF */
774 774 /*
775 775 * kpcos(x)= cos(pi*x)/pi
776 776 * 2 4 6 8 10 12
777 777 * = 1/pi +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x +kc[4]*x +kc[5]*x
778 778 *
779 779 * 2 4 6 8 10 22
780 780 * = 1/pi - pi/2*x +kc[0]*x +kc[1]*x +kc[2]*x +kc[3]*x +...+kc[9]*x
781 781 *
782 782 * -pi/2*x*x = (npi_2_h + npi_2_l) * (x_f+x_l)*(x_f+x_l)
783 783 * = npi_2_h*(x_f+x_l)*(x_f+x_l) + npi_2_l*x*x
784 784 * = npi_2_h*x_f*x_f + npi_2_h*(x*x-x_f*x_f) + npi_2_l*x*x
785 785 * = npi_2_h*x_f*x_f + npi_2_h*(x+x_f)*(x-x_f) + npi_2_l*x*x
786 786 * Here x_f = (long double) (float)x
787 787 * Note that pi/2(in hex) =
788 788 * 1.921FB54442D18469898CC51701B839A252049C1114CF98E804177D4C76273644A29
789 789 * npi_2_h = -pi/2 chopped to 25 bits = -1.921FB50000000000000000000000000 =
790 790 * -1.570796310901641845703125000000000 and
791 791 * npi_2_l =
792 792 * -0.0000004442D18469898CC51701B839A252049C1114CF98E804177D4C76273644A29 =
793 793 * -.0000000158932547735281966916397514420985846996875529104874722961539 =
794 794 * -1.5893254773528196691639751442098584699687552910487472296153e-8
795 795 * 1/pi(in hex) =
796 796 * .517CC1B727220A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D5EF5DE2B
797 797 * will be splitted into:
798 798 * one_pi_h = 1/pi chopped to 48 bits = .517CC1B727220000000000... and
799 799 * one_pi_l = .0000000000000A94FE13ABE8FA9A6EE06DB14ACC9E21C820FF28B1D5EF5DE2B
800 800 */
801 801
802 802 static const long double
803 803 #if defined(__x86)
804 804 one_pi_h = 0.3183098861481994390487670898437500L, /* 31 bits */
805 805 one_pi_l = 3.559123248900043690127872406891929148e-11L,
806 806 #else
807 807 one_pi_h = 0.31830988618379052468299050815403461456298828125L,
808 808 one_pi_l = 1.46854777018590994109505931010230912897495334688117e-16L,
809 809 #endif
810 810 npi_2_h = -1.570796310901641845703125000000000L,
811 811 npi_2_l = -1.5893254773528196691639751442098584699687552910e-8L;
812 812
813 813 static const long double kc[] = {
814 814 +1.29192819501249250731151312779548918765320728489e+0000L,
815 815 -4.25027339979557573976029596929319207009444090366e-0001L,
816 816 +7.49080661650990096109672954618317623888421628613e-0002L,
817 817 -8.21458866111282287985539464173976555436050215120e-0003L,
818 818 +6.14202578809529228503205255165761204750211603402e-0004L,
819 819 -3.33073432691149607007217330302595267179545908740e-0005L,
820 820 +1.36970959047832085796809745461530865597993680204e-0006L,
821 821 -4.41780774262583514450246512727201806217271097336e-0008L,
822 822 +1.14741409212381858820016567664488123478660705759e-0009L,
823 823 -2.44261236114707374558437500654381006300502749632e-0011L,
824 824 };
825 825 /* INDENT ON */
826 826
827 827 /*
828 828 * assume x is not tiny and positive
829 829 */
830 830 static struct LDouble
831 831 kpcos(long double x) {
832 832 long double z, t1, t2, t3, t4, x4, x8;
833 833 int i;
834 834 struct LDouble xx;
835 835
836 836 z = x * x;
837 837 xx.h = one_pi_h;
838 838 t1 = (long double) ((float) x);
839 839 x4 = z * z;
840 840 t2 = npi_2_l * z + npi_2_h * (x + t1) * (x - t1);
841 841 for (i = 8, t3 = kc[9]; i >= 0; i--)
842 842 t3 = z * t3 + kc[i];
843 843 t3 = one_pi_l + x4 * t3;
844 844 t4 = t1 * t1 * npi_2_h;
845 845 x8 = t2 + t3;
846 846 xx.l = x8 + t4;
847 847 return (xx);
848 848 }
849 849
850 850 /* INDENT OFF */
851 851 static const long double
852 852 /* 0.13486180573279076968979393577465291700642511139552429398233 */
853 853 #if defined(__x86)
854 854 t0z1 = 0.1348618057327907696779385054997035808810L,
855 855 t0z1_l = 1.1855430274949336125392717150257379614654e-20L,
856 856 #else
857 857 t0z1 = 0.1348618057327907696897939357746529168654L,
858 858 t0z1_l = 1.4102088588676879418739164486159514674310e-37L,
859 859 #endif
860 860 /* 0.46163214496836234126265954232572132846819620400644635129599 */
861 861 #if defined(__x86)
862 862 t0z2 = 0.4616321449683623412538115843295472018326L,
863 863 t0z2_l = 8.84795799617412663558532305039261747030640e-21L,
864 864 #else
865 865 t0z2 = 0.46163214496836234126265954232572132343318L,
866 866 t0z2_l = 5.03501162329616380465302666480916271611101e-36L,
867 867 #endif
868 868 /* 0.81977310110050060178786870492160699631174407846245179119586 */
869 869 #if defined(__x86)
870 870 t0z3 = 0.81977310110050060178773362329351925836817L,
871 871 t0z3_l = 1.350816280877379435658077052534574556256230e-22L
872 872 #else
873 873 t0z3 = 0.8197731011005006017878687049216069516957449L,
|
↓ open down ↓ |
873 lines elided |
↑ open up ↑ |
874 874 t0z3_l = 4.461599916947014419045492615933551648857380e-35L
875 875 #endif
876 876 ;
877 877 /* INDENT ON */
878 878
879 879 /*
880 880 * gamma(x+i) for 0 <= x < 1
881 881 */
882 882 static struct LDouble
883 883 gam_n(int i, long double x) {
884 - struct LDouble rr, yy;
884 + struct LDouble rr = {0.0L, 0.0L}, yy;
885 885 long double r1, r2, t2, z, xh, xl, yh, yl, zh, z1, z2, zl, x5, wh, wl;
886 886
887 887 /* compute yy = gamma(x+1) */
888 888 if (x > 0.2845L) {
889 889 if (x > 0.6374L) {
890 890 r1 = x - t0z3;
891 891 r2 = CHOPPED((r1 - t0z3_l));
892 892 t2 = r1 - r2;
893 893 yy = GT3(r2, t2 - t0z3_l);
894 894 } else {
895 895 r1 = x - t0z2;
896 896 r2 = CHOPPED((r1 - t0z2_l));
897 897 t2 = r1 - r2;
898 898 yy = GT2(r2, t2 - t0z2_l);
899 899 }
900 900 } else {
901 901 r1 = x - t0z1;
902 902 r2 = CHOPPED((r1 - t0z1_l));
903 903 t2 = r1 - r2;
904 904 yy = GT1(r2, t2 - t0z1_l);
905 905 }
906 906 /* compute gamma(x+i) = (x+i-1)*...*(x+1)*yy, 0<i<8 */
907 907 switch (i) {
908 908 case 0: /* yy/x */
909 909 r1 = one / x;
910 910 xh = CHOPPED((x)); /* x is not tiny */
911 911 rr.h = CHOPPED(((yy.h + yy.l) * r1));
912 912 rr.l = r1 * (yy.h - rr.h * xh) - ((r1 * rr.h) * (x - xh) -
913 913 r1 * yy.l);
914 914 break;
915 915 case 1: /* yy */
916 916 rr.h = yy.h;
917 917 rr.l = yy.l;
918 918 break;
919 919 case 2: /* (x+1)*yy */
920 920 z = x + one; /* may not be exact */
921 921 zh = CHOPPED((z));
922 922 rr.h = zh * yy.h;
923 923 rr.l = z * yy.l + (x - (zh - one)) * yy.h;
924 924 break;
925 925 case 3: /* (x+2)*(x+1)*yy */
926 926 z1 = x + one;
927 927 z2 = x + 2.0L;
928 928 z = z1 * z2;
929 929 xh = CHOPPED((z));
930 930 zh = CHOPPED((z1));
931 931 xl = (x - (zh - one)) * (z2 + zh) - (xh - zh * (zh + one));
932 932
933 933 rr.h = xh * yy.h;
934 934 rr.l = z * yy.l + xl * yy.h;
935 935 break;
936 936
937 937 case 4: /* (x+1)*(x+3)*(x+2)*yy */
938 938 z1 = x + 2.0L;
939 939 z2 = (x + one) * (x + 3.0L);
940 940 zh = CHOPPED(z1);
941 941 zl = x - (zh - 2.0L);
942 942 xh = CHOPPED(z2);
943 943 xl = zl * (zh + z1) - (xh - (zh * zh - one));
944 944
945 945 /* wh+wl=(x+2)*yy */
946 946 wh = CHOPPED((z1 * (yy.h + yy.l)));
947 947 wl = (zl * yy.h + z1 * yy.l) - (wh - zh * yy.h);
948 948
949 949 rr.h = xh * wh;
950 950 rr.l = z2 * wl + xl * wh;
951 951
952 952 break;
953 953 case 5: /* ((x+1)*(x+4)*(x+2)*(x+3))*yy */
954 954 z1 = x + 2.0L;
955 955 z2 = x + 3.0L;
956 956 z = z1 * z2;
957 957 zh = CHOPPED((z1));
958 958 yh = CHOPPED((z));
959 959 yl = (x - (zh - 2.0L)) * (z2 + zh) - (yh - zh * (zh + one));
960 960 z2 = z - 2.0L;
961 961 z *= z2;
962 962 xh = CHOPPED((z));
963 963 xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L));
964 964 rr.h = xh * yy.h;
965 965 rr.l = z * yy.l + xl * yy.h;
966 966 break;
967 967 case 6: /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5))*yy */
968 968 z1 = x + 2.0L;
969 969 z2 = x + 3.0L;
970 970 z = z1 * z2;
971 971 zh = CHOPPED((z1));
972 972 yh = CHOPPED((z));
973 973 z1 = x - (zh - 2.0L);
974 974 yl = z1 * (z2 + zh) - (yh - zh * (zh + one));
975 975 z2 = z - 2.0L;
976 976 x5 = x + 5.0L;
977 977 z *= z2;
978 978 xh = CHOPPED(z);
979 979 zh += 3.0;
980 980 xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L));
981 981 /* xh+xl=(x+1)*...*(x+4) */
982 982 /* wh+wl=(x+5)*yy */
983 983 wh = CHOPPED((x5 * (yy.h + yy.l)));
984 984 wl = (z1 * yy.h + x5 * yy.l) - (wh - zh * yy.h);
985 985 rr.h = wh * xh;
986 986 rr.l = z * wl + xl * wh;
987 987 break;
988 988 case 7: /* ((x+1)*(x+2)*(x+3)*(x+4)*(x+5)*(x+6))*yy */
989 989 z1 = x + 3.0L;
990 990 z2 = x + 4.0L;
991 991 z = z2 * z1;
992 992 zh = CHOPPED((z1));
993 993 yh = CHOPPED((z)); /* yh+yl = (x+3)(x+4) */
994 994 yl = (x - (zh - 3.0L)) * (z2 + zh) - (yh - (zh * (zh + one)));
995 995 z1 = x + 6.0L;
996 996 z2 = z - 2.0L; /* z2 = (x+2)*(x+5) */
997 997 z *= z2;
998 998 xh = CHOPPED((z));
999 999 xl = yl * (z2 + yh) - (xh - yh * (yh - 2.0L));
1000 1000 /* xh+xl=(x+2)*...*(x+5) */
1001 1001 /* wh+wl=(x+1)(x+6)*yy */
1002 1002 z2 -= 4.0L; /* z2 = (x+1)(x+6) */
1003 1003 wh = CHOPPED((z2 * (yy.h + yy.l)));
1004 1004 wl = (z2 * yy.l + yl * yy.h) - (wh - (yh - 6.0L) * yy.h);
1005 1005 rr.h = wh * xh;
1006 1006 rr.l = z * wl + xl * wh;
1007 1007 }
1008 1008 return (rr);
1009 1009 }
1010 1010
1011 1011 long double
1012 1012 tgammal(long double x) {
1013 1013 struct LDouble ss, ww;
1014 1014 long double t, t1, t2, t3, t4, t5, w, y, z, z1, z2, z3, z5;
1015 1015 int i, j, m, ix, hx, xk;
1016 1016 unsigned lx;
1017 1017
1018 1018 hx = H0_WORD(x);
1019 1019 lx = H3_WORD(x);
1020 1020 ix = hx & 0x7fffffff;
1021 1021 y = x;
1022 1022 if (ix < 0x3f8e0000) { /* x < 2**-113 */
1023 1023 return (one / x);
1024 1024 }
|
↓ open down ↓ |
130 lines elided |
↑ open up ↑ |
1025 1025 if (ix >= 0x7fff0000)
1026 1026 return (x * ((hx < 0)? zero : x)); /* Inf or NaN */
1027 1027 if (x > overflow) /* overflow threshold */
1028 1028 return (x * 1.0e4932L);
1029 1029 if (hx >= 0x40020000) { /* x >= 8 */
1030 1030 ww = large_gam(x, &m);
1031 1031 w = ww.h + ww.l;
1032 1032 return (scalbnl(w, m));
1033 1033 }
1034 1034
1035 - if (hx > 0) { /* x from 0 to 8 */
1035 + if (hx > 0) { /* 0 < x < 8 */
1036 1036 i = (int) x;
1037 1037 ww = gam_n(i, x - (long double) i);
1038 1038 return (ww.h + ww.l);
1039 1039 }
1040 1040 /* INDENT OFF */
1041 1041 /* negative x */
1042 1042 /*
1043 1043 * compute xk =
1044 1044 * -2 ... x is an even int (-inf is considered an even #)
1045 1045 * -1 ... x is an odd int
1046 1046 * +0 ... x is not an int but chopped to an even int
1047 1047 * +1 ... x is not an int but chopped to an odd int
1048 1048 */
1049 1049 /* INDENT ON */
1050 1050 xk = 0;
1051 1051 #if defined(__x86)
1052 1052 if (ix >= 0x403e0000) { /* x >= 2**63 } */
1053 1053 if (ix >= 0x403f0000)
1054 1054 xk = -2;
1055 1055 else
1056 1056 xk = -2 + (lx & 1);
1057 1057 #else
1058 1058 if (ix >= 0x406f0000) { /* x >= 2**112 */
1059 1059 if (ix >= 0x40700000)
1060 1060 xk = -2;
1061 1061 else
1062 1062 xk = -2 + (lx & 1);
1063 1063 #endif
1064 1064 } else if (ix >= 0x3fff0000) {
1065 1065 w = -x;
1066 1066 t1 = floorl(w);
1067 1067 t2 = t1 * half;
1068 1068 t3 = floorl(t2);
1069 1069 if (t1 == w) {
1070 1070 if (t2 == t3)
1071 1071 xk = -2;
1072 1072 else
1073 1073 xk = -1;
1074 1074 } else {
1075 1075 if (t2 == t3)
1076 1076 xk = 0;
1077 1077 else
1078 1078 xk = 1;
1079 1079 }
1080 1080 }
1081 1081
1082 1082 if (xk < 0) {
1083 1083 /* return NaN. Ideally gamma(-n)= (-1)**(n+1) * inf */
1084 1084 return (x - x) / (x - x);
1085 1085 }
1086 1086
1087 1087 /*
1088 1088 * negative underflow thresold -(1774+9ulp)
1089 1089 */
1090 1090 if (x < -1774.0000000000000000000000000000017749370L) {
1091 1091 z = tiny / x;
1092 1092 if (xk == 1)
1093 1093 z = -z;
1094 1094 return (z * tiny);
1095 1095 }
1096 1096
1097 1097 /* INDENT OFF */
1098 1098 /*
1099 1099 * now compute gamma(x) by -1/((sin(pi*y)/pi)*gamma(1+y)), y = -x
1100 1100 */
1101 1101 /*
1102 1102 * First compute ss = -sin(pi*y)/pi so that
1103 1103 * gamma(x) = 1/(ss*gamma(1+y))
1104 1104 */
1105 1105 /* INDENT ON */
1106 1106 y = -x;
1107 1107 j = (int) y;
1108 1108 z = y - (long double) j;
1109 1109 if (z > 0.3183098861837906715377675L)
1110 1110 if (z > 0.6816901138162093284622325L)
1111 1111 ss = kpsin(one - z);
1112 1112 else
1113 1113 ss = kpcos(0.5L - z);
1114 1114 else
1115 1115 ss = kpsin(z);
1116 1116 if (xk == 0) {
1117 1117 ss.h = -ss.h;
1118 1118 ss.l = -ss.l;
1119 1119 }
1120 1120
1121 1121 /* Then compute ww = gamma(1+y), note that result scale to 2**m */
1122 1122 m = 0;
1123 1123 if (j < 7) {
1124 1124 ww = gam_n(j + 1, z);
1125 1125 } else {
1126 1126 w = y + one;
1127 1127 if ((lx & 1) == 0) { /* y+1 exact (note that y<184) */
1128 1128 ww = large_gam(w, &m);
1129 1129 } else {
1130 1130 t = w - one;
1131 1131 if (t == y) { /* y+one exact */
1132 1132 ww = large_gam(w, &m);
1133 1133 } else { /* use y*gamma(y) */
1134 1134 if (j == 7)
1135 1135 ww = gam_n(j, z);
1136 1136 else
1137 1137 ww = large_gam(y, &m);
1138 1138 t4 = ww.h + ww.l;
1139 1139 t1 = CHOPPED((y));
1140 1140 t2 = CHOPPED((t4));
1141 1141 /* t4 will not be too large */
1142 1142 ww.l = y * (ww.l - (t2 - ww.h)) + (y - t1) * t2;
1143 1143 ww.h = t1 * t2;
1144 1144 }
1145 1145 }
1146 1146 }
1147 1147
1148 1148 /* compute 1/(ss*ww) */
1149 1149 t3 = ss.h + ss.l;
1150 1150 t4 = ww.h + ww.l;
1151 1151 t1 = CHOPPED((t3));
1152 1152 t2 = CHOPPED((t4));
1153 1153 z1 = ss.l - (t1 - ss.h); /* (t1,z1) = ss */
1154 1154 z2 = ww.l - (t2 - ww.h); /* (t2,z2) = ww */
1155 1155 t3 = t3 * t4; /* t3 = ss*ww */
1156 1156 z3 = one / t3; /* z3 = 1/(ss*ww) */
1157 1157 t5 = t1 * t2;
1158 1158 z5 = z1 * t4 + t1 * z2; /* (t5,z5) = ss*ww */
1159 1159 t1 = CHOPPED((t3)); /* (t1,z1) = ss*ww */
1160 1160 z1 = z5 - (t1 - t5);
1161 1161 t2 = CHOPPED((z3)); /* leading 1/(ss*ww) */
1162 1162 z2 = z3 * (t2 * z1 - (one - t2 * t1));
1163 1163 z = t2 - z2;
1164 1164
1165 1165 return (scalbnl(z, -m));
1166 1166 }
|
↓ open down ↓ |
121 lines elided |
↑ open up ↑ |
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX