Mercurial > hg > CbC > CbC_gcc

comparison libquadmath/math/lgammaq.c @ 68:561a7518be6b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update gcc-4.6
author Nobuyasu Oshiro <dimolto@cr.ie.u-ryukyu.ac.jp>
date Sun, 21 Aug 2011 07:07:55 +0900
parents
children 04ced10e8804
comparison
equal deleted inserted replaced
67:f6334be47118 68:561a7518be6b
1 /* lgammal
2 *
3 * Natural logarithm of gamma function
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * __float128 x, y, lgammal();
10 * extern int sgngam;
11 *
12 * y = lgammal(x);
13 *
14 *
15 *
16 * DESCRIPTION:
17 *
18 * Returns the base e (2.718...) logarithm of the absolute
19 * value of the gamma function of the argument.
20 * The sign (+1 or -1) of the gamma function is returned in a
21 * global (extern) variable named sgngam.
22 *
23 * The positive domain is partitioned into numerous segments for approximation.
24 * For x > 10,
25 * log gamma(x) = (x - 0.5) log(x) - x + log sqrt(2 pi) + 1/x R(1/x^2)
26 * Near the minimum at x = x0 = 1.46... the approximation is
27 * log gamma(x0 + z) = log gamma(x0) + z^2 P(z)/Q(z)
28 * for small z.
29 * Elsewhere between 0 and 10,
30 * log gamma(n + z) = log gamma(n) + z P(z)/Q(z)
31 * for various selected n and small z.
32 *
33 * The cosecant reflection formula is employed for negative arguments.
34 *
35 *
36 *
37 * ACCURACY:
38 *
39 *
40 * arithmetic domain # trials peak rms
41 * Relative error:
42 * IEEE 10, 30 100000 3.9e-34 9.8e-35
43 * IEEE 0, 10 100000 3.8e-34 5.3e-35
44 * Absolute error:
45 * IEEE -10, 0 100000 8.0e-34 8.0e-35
46 * IEEE -30, -10 100000 4.4e-34 1.0e-34
47 * IEEE -100, 100 100000 1.0e-34
48 *
49 * The absolute error criterion is the same as relative error
50 * when the function magnitude is greater than one but it is absolute
51 * when the magnitude is less than one.
52 *
53 */
54
55 /* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
56
57 This library is free software; you can redistribute it and/or
58 modify it under the terms of the GNU Lesser General Public
59 License as published by the Free Software Foundation; either
60 version 2.1 of the License, or (at your option) any later version.
61
62 This library is distributed in the hope that it will be useful,
63 but WITHOUT ANY WARRANTY; without even the implied warranty of
64 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
65 Lesser General Public License for more details.
66
67 You should have received a copy of the GNU Lesser General Public
68 License along with this library; if not, write to the Free Software
69 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
70
71 #include "quadmath-imp.h"
72
73 static const __float128 PIQ = 3.1415926535897932384626433832795028841972E0Q;
74 static const __float128 MAXLGM = 1.0485738685148938358098967157129705071571E4928Q;
75 static const __float128 one = 1.0Q;
76 static const __float128 zero = 0.0Q;
77 static const __float128 huge = 1.0e4000Q;
78
79 /* log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x P(1/x^2)
80 1/x <= 0.0741 (x >= 13.495...)
81 Peak relative error 1.5e-36 */
82 static const __float128 ls2pi = 9.1893853320467274178032973640561763986140E-1Q;
83 #define NRASY 12
84 static const __float128 RASY[NRASY + 1] =
85 {
86 8.333333333333333333333333333310437112111E-2Q,
87 -2.777777777777777777777774789556228296902E-3Q,
88 7.936507936507936507795933938448586499183E-4Q,
89 -5.952380952380952041799269756378148574045E-4Q,
90 8.417508417507928904209891117498524452523E-4Q,
91 -1.917526917481263997778542329739806086290E-3Q,
92 6.410256381217852504446848671499409919280E-3Q,
93 -2.955064066900961649768101034477363301626E-2Q,
94 1.796402955865634243663453415388336954675E-1Q,
95 -1.391522089007758553455753477688592767741E0Q,
96 1.326130089598399157988112385013829305510E1Q,
97 -1.420412699593782497803472576479997819149E2Q,
98 1.218058922427762808938869872528846787020E3Q
99 };
100
101
102 /* log gamma(x+13) = log gamma(13) + x P(x)/Q(x)
103 -0.5 <= x <= 0.5
104 12.5 <= x+13 <= 13.5
105 Peak relative error 1.1e-36 */
106 static const __float128 lgam13a = 1.9987213134765625E1Q;
107 static const __float128 lgam13b = 1.3608962611495173623870550785125024484248E-6Q;
108 #define NRN13 7
109 static const __float128 RN13[NRN13 + 1] =
110 {
111 8.591478354823578150238226576156275285700E11Q,
112 2.347931159756482741018258864137297157668E11Q,
113 2.555408396679352028680662433943000804616E10Q,
114 1.408581709264464345480765758902967123937E9Q,
115 4.126759849752613822953004114044451046321E7Q,
116 6.133298899622688505854211579222889943778E5Q,
117 3.929248056293651597987893340755876578072E3Q,
118 6.850783280018706668924952057996075215223E0Q
119 };
120 #define NRD13 6
121 static const __float128 RD13[NRD13 + 1] =
122 {
123 3.401225382297342302296607039352935541669E11Q,
124 8.756765276918037910363513243563234551784E10Q,
125 8.873913342866613213078554180987647243903E9Q,
126 4.483797255342763263361893016049310017973E8Q,
127 1.178186288833066430952276702931512870676E7Q,
128 1.519928623743264797939103740132278337476E5Q,
129 7.989298844938119228411117593338850892311E2Q
130 /* 1.0E0Q */
131 };
132
133
134 /* log gamma(x+12) = log gamma(12) + x P(x)/Q(x)
135 -0.5 <= x <= 0.5
136 11.5 <= x+12 <= 12.5
137 Peak relative error 4.1e-36 */
138 static const __float128 lgam12a = 1.75023040771484375E1Q;
139 static const __float128 lgam12b = 3.7687254483392876529072161996717039575982E-6Q;
140 #define NRN12 7
141 static const __float128 RN12[NRN12 + 1] =
142 {
143 4.709859662695606986110997348630997559137E11Q,
144 1.398713878079497115037857470168777995230E11Q,
145 1.654654931821564315970930093932954900867E10Q,
146 9.916279414876676861193649489207282144036E8Q,
147 3.159604070526036074112008954113411389879E7Q,
148 5.109099197547205212294747623977502492861E5Q,
149 3.563054878276102790183396740969279826988E3Q,
150 6.769610657004672719224614163196946862747E0Q
151 };
152 #define NRD12 6
153 static const __float128 RD12[NRD12 + 1] =
154 {
155 1.928167007860968063912467318985802726613E11Q,
156 5.383198282277806237247492369072266389233E10Q,
157 5.915693215338294477444809323037871058363E9Q,
158 3.241438287570196713148310560147925781342E8Q,
159 9.236680081763754597872713592701048455890E6Q,
160 1.292246897881650919242713651166596478850E5Q,
161 7.366532445427159272584194816076600211171E2Q
162 /* 1.0E0Q */
163 };
164
165
166 /* log gamma(x+11) = log gamma(11) + x P(x)/Q(x)
167 -0.5 <= x <= 0.5
168 10.5 <= x+11 <= 11.5
169 Peak relative error 1.8e-35 */
170 static const __float128 lgam11a = 1.5104400634765625E1Q;
171 static const __float128 lgam11b = 1.1938309890295225709329251070371882250744E-5Q;
172 #define NRN11 7
173 static const __float128 RN11[NRN11 + 1] =
174 {
175 2.446960438029415837384622675816736622795E11Q,
176 7.955444974446413315803799763901729640350E10Q,
177 1.030555327949159293591618473447420338444E10Q,
178 6.765022131195302709153994345470493334946E8Q,
179 2.361892792609204855279723576041468347494E7Q,
180 4.186623629779479136428005806072176490125E5Q,
181 3.202506022088912768601325534149383594049E3Q,
182 6.681356101133728289358838690666225691363E0Q
183 };
184 #define NRD11 6
185 static const __float128 RD11[NRD11 + 1] =
186 {
187 1.040483786179428590683912396379079477432E11Q,
188 3.172251138489229497223696648369823779729E10Q,
189 3.806961885984850433709295832245848084614E9Q,
190 2.278070344022934913730015420611609620171E8Q,
191 7.089478198662651683977290023829391596481E6Q,
192 1.083246385105903533237139380509590158658E5Q,
193 6.744420991491385145885727942219463243597E2Q
194 /* 1.0E0Q */
195 };
196
197
198 /* log gamma(x+10) = log gamma(10) + x P(x)/Q(x)
199 -0.5 <= x <= 0.5
200 9.5 <= x+10 <= 10.5
201 Peak relative error 5.4e-37 */
202 static const __float128 lgam10a = 1.280181884765625E1Q;
203 static const __float128 lgam10b = 8.6324252196112077178745667061642811492557E-6Q;
204 #define NRN10 7
205 static const __float128 RN10[NRN10 + 1] =
206 {
207 -1.239059737177249934158597996648808363783E14Q,
208 -4.725899566371458992365624673357356908719E13Q,
209 -7.283906268647083312042059082837754850808E12Q,
210 -5.802855515464011422171165179767478794637E11Q,
211 -2.532349691157548788382820303182745897298E10Q,
212 -5.884260178023777312587193693477072061820E8Q,
213 -6.437774864512125749845840472131829114906E6Q,
214 -2.350975266781548931856017239843273049384E4Q
215 };
216 #define NRD10 7
217 static const __float128 RD10[NRD10 + 1] =
218 {
219 -5.502645997581822567468347817182347679552E13Q,
220 -1.970266640239849804162284805400136473801E13Q,
221 -2.819677689615038489384974042561531409392E12Q,
222 -2.056105863694742752589691183194061265094E11Q,
223 -8.053670086493258693186307810815819662078E9Q,
224 -1.632090155573373286153427982504851867131E8Q,
225 -1.483575879240631280658077826889223634921E6Q,
226 -4.002806669713232271615885826373550502510E3Q
227 /* 1.0E0Q */
228 };
229
230
231 /* log gamma(x+9) = log gamma(9) + x P(x)/Q(x)
232 -0.5 <= x <= 0.5
233 8.5 <= x+9 <= 9.5
234 Peak relative error 3.6e-36 */
235 static const __float128 lgam9a = 1.06045989990234375E1Q;
236 static const __float128 lgam9b = 3.9037218127284172274007216547549861681400E-6Q;
237 #define NRN9 7
238 static const __float128 RN9[NRN9 + 1] =
239 {
240 -4.936332264202687973364500998984608306189E13Q,
241 -2.101372682623700967335206138517766274855E13Q,
242 -3.615893404644823888655732817505129444195E12Q,
243 -3.217104993800878891194322691860075472926E11Q,
244 -1.568465330337375725685439173603032921399E10Q,
245 -4.073317518162025744377629219101510217761E8Q,
246 -4.983232096406156139324846656819246974500E6Q,
247 -2.036280038903695980912289722995505277253E4Q
248 };
249 #define NRD9 7
250 static const __float128 RD9[NRD9 + 1] =
251 {
252 -2.306006080437656357167128541231915480393E13Q,
253 -9.183606842453274924895648863832233799950E12Q,
254 -1.461857965935942962087907301194381010380E12Q,
255 -1.185728254682789754150068652663124298303E11Q,
256 -5.166285094703468567389566085480783070037E9Q,
257 -1.164573656694603024184768200787835094317E8Q,
258 -1.177343939483908678474886454113163527909E6Q,
259 -3.529391059783109732159524500029157638736E3Q
260 /* 1.0E0Q */
261 };
262
263
264 /* log gamma(x+8) = log gamma(8) + x P(x)/Q(x)
265 -0.5 <= x <= 0.5
266 7.5 <= x+8 <= 8.5
267 Peak relative error 2.4e-37 */
268 static const __float128 lgam8a = 8.525146484375E0Q;
269 static const __float128 lgam8b = 1.4876690414300165531036347125050759667737E-5Q;
270 #define NRN8 8
271 static const __float128 RN8[NRN8 + 1] =
272 {
273 6.600775438203423546565361176829139703289E11Q,
274 3.406361267593790705240802723914281025800E11Q,
275 7.222460928505293914746983300555538432830E10Q,
276 8.102984106025088123058747466840656458342E9Q,
277 5.157620015986282905232150979772409345927E8Q,
278 1.851445288272645829028129389609068641517E7Q,
279 3.489261702223124354745894067468953756656E5Q,
280 2.892095396706665774434217489775617756014E3Q,
281 6.596977510622195827183948478627058738034E0Q
282 };
283 #define NRD8 7
284 static const __float128 RD8[NRD8 + 1] =
285 {
286 3.274776546520735414638114828622673016920E11Q,
287 1.581811207929065544043963828487733970107E11Q,
288 3.108725655667825188135393076860104546416E10Q,
289 3.193055010502912617128480163681842165730E9Q,
290 1.830871482669835106357529710116211541839E8Q,
291 5.790862854275238129848491555068073485086E6Q,
292 9.305213264307921522842678835618803553589E4Q,
293 6.216974105861848386918949336819572333622E2Q
294 /* 1.0E0Q */
295 };
296
297
298 /* log gamma(x+7) = log gamma(7) + x P(x)/Q(x)
299 -0.5 <= x <= 0.5
300 6.5 <= x+7 <= 7.5
301 Peak relative error 3.2e-36 */
302 static const __float128 lgam7a = 6.5792388916015625E0Q;
303 static const __float128 lgam7b = 1.2320408538495060178292903945321122583007E-5Q;
304 #define NRN7 8
305 static const __float128 RN7[NRN7 + 1] =
306 {
307 2.065019306969459407636744543358209942213E11Q,
308 1.226919919023736909889724951708796532847E11Q,
309 2.996157990374348596472241776917953749106E10Q,
310 3.873001919306801037344727168434909521030E9Q,
311 2.841575255593761593270885753992732145094E8Q,
312 1.176342515359431913664715324652399565551E7Q,
313 2.558097039684188723597519300356028511547E5Q,
314 2.448525238332609439023786244782810774702E3Q,
315 6.460280377802030953041566617300902020435E0Q
316 };
317 #define NRD7 7
318 static const __float128 RD7[NRD7 + 1] =
319 {
320 1.102646614598516998880874785339049304483E11Q,
321 6.099297512712715445879759589407189290040E10Q,
322 1.372898136289611312713283201112060238351E10Q,
323 1.615306270420293159907951633566635172343E9Q,
324 1.061114435798489135996614242842561967459E8Q,
325 3.845638971184305248268608902030718674691E6Q,
326 7.081730675423444975703917836972720495507E4Q,
327 5.423122582741398226693137276201344096370E2Q
328 /* 1.0E0Q */
329 };
330
331
332 /* log gamma(x+6) = log gamma(6) + x P(x)/Q(x)
333 -0.5 <= x <= 0.5
334 5.5 <= x+6 <= 6.5
335 Peak relative error 6.2e-37 */
336 static const __float128 lgam6a = 4.7874908447265625E0Q;
337 static const __float128 lgam6b = 8.9805548349424770093452324304839959231517E-7Q;
338 #define NRN6 8
339 static const __float128 RN6[NRN6 + 1] =
340 {
341 -3.538412754670746879119162116819571823643E13Q,
342 -2.613432593406849155765698121483394257148E13Q,
343 -8.020670732770461579558867891923784753062E12Q,
344 -1.322227822931250045347591780332435433420E12Q,
345 -1.262809382777272476572558806855377129513E11Q,
346 -7.015006277027660872284922325741197022467E9Q,
347 -2.149320689089020841076532186783055727299E8Q,
348 -3.167210585700002703820077565539658995316E6Q,
349 -1.576834867378554185210279285358586385266E4Q
350 };
351 #define NRD6 8
352 static const __float128 RD6[NRD6 + 1] =
353 {
354 -2.073955870771283609792355579558899389085E13Q,
355 -1.421592856111673959642750863283919318175E13Q,
356 -4.012134994918353924219048850264207074949E12Q,
357 -6.013361045800992316498238470888523722431E11Q,
358 -5.145382510136622274784240527039643430628E10Q,
359 -2.510575820013409711678540476918249524123E9Q,
360 -6.564058379709759600836745035871373240904E7Q,
361 -7.861511116647120540275354855221373571536E5Q,
362 -2.821943442729620524365661338459579270561E3Q
363 /* 1.0E0Q */
364 };
365
366
367 /* log gamma(x+5) = log gamma(5) + x P(x)/Q(x)
368 -0.5 <= x <= 0.5
369 4.5 <= x+5 <= 5.5
370 Peak relative error 3.4e-37 */
371 static const __float128 lgam5a = 3.17803955078125E0Q;
372 static const __float128 lgam5b = 1.4279566695619646941601297055408873990961E-5Q;
373 #define NRN5 9
374 static const __float128 RN5[NRN5 + 1] =
375 {
376 2.010952885441805899580403215533972172098E11Q,
377 1.916132681242540921354921906708215338584E11Q,
378 7.679102403710581712903937970163206882492E10Q,
379 1.680514903671382470108010973615268125169E10Q,
380 2.181011222911537259440775283277711588410E9Q,
381 1.705361119398837808244780667539728356096E8Q,
382 7.792391565652481864976147945997033946360E6Q,
383 1.910741381027985291688667214472560023819E5Q,
384 2.088138241893612679762260077783794329559E3Q,
385 6.330318119566998299106803922739066556550E0Q
386 };
387 #define NRD5 8
388 static const __float128 RD5[NRD5 + 1] =
389 {
390 1.335189758138651840605141370223112376176E11Q,
391 1.174130445739492885895466097516530211283E11Q,
392 4.308006619274572338118732154886328519910E10Q,
393 8.547402888692578655814445003283720677468E9Q,
394 9.934628078575618309542580800421370730906E8Q,
395 6.847107420092173812998096295422311820672E7Q,
396 2.698552646016599923609773122139463150403E6Q,
397 5.526516251532464176412113632726150253215E4Q,
398 4.772343321713697385780533022595450486932E2Q
399 /* 1.0E0Q */
400 };
401
402
403 /* log gamma(x+4) = log gamma(4) + x P(x)/Q(x)
404 -0.5 <= x <= 0.5
405 3.5 <= x+4 <= 4.5
406 Peak relative error 6.7e-37 */
407 static const __float128 lgam4a = 1.791748046875E0Q;
408 static const __float128 lgam4b = 1.1422353055000812477358380702272722990692E-5Q;
409 #define NRN4 9
410 static const __float128 RN4[NRN4 + 1] =
411 {
412 -1.026583408246155508572442242188887829208E13Q,
413 -1.306476685384622809290193031208776258809E13Q,
414 -7.051088602207062164232806511992978915508E12Q,
415 -2.100849457735620004967624442027793656108E12Q,
416 -3.767473790774546963588549871673843260569E11Q,
417 -4.156387497364909963498394522336575984206E10Q,
418 -2.764021460668011732047778992419118757746E9Q,
419 -1.036617204107109779944986471142938641399E8Q,
420 -1.895730886640349026257780896972598305443E6Q,
421 -1.180509051468390914200720003907727988201E4Q
422 };
423 #define NRD4 9
424 static const __float128 RD4[NRD4 + 1] =
425 {
426 -8.172669122056002077809119378047536240889E12Q,
427 -9.477592426087986751343695251801814226960E12Q,
428 -4.629448850139318158743900253637212801682E12Q,
429 -1.237965465892012573255370078308035272942E12Q,
430 -1.971624313506929845158062177061297598956E11Q,
431 -1.905434843346570533229942397763361493610E10Q,
432 -1.089409357680461419743730978512856675984E9Q,
433 -3.416703082301143192939774401370222822430E7Q,
434 -4.981791914177103793218433195857635265295E5Q,
435 -2.192507743896742751483055798411231453733E3Q
436 /* 1.0E0Q */
437 };
438
439
440 /* log gamma(x+3) = log gamma(3) + x P(x)/Q(x)
441 -0.25 <= x <= 0.5
442 2.75 <= x+3 <= 3.5
443 Peak relative error 6.0e-37 */
444 static const __float128 lgam3a = 6.93145751953125E-1Q;
445 static const __float128 lgam3b = 1.4286068203094172321214581765680755001344E-6Q;
446
447 #define NRN3 9
448 static const __float128 RN3[NRN3 + 1] =
449 {
450 -4.813901815114776281494823863935820876670E11Q,
451 -8.425592975288250400493910291066881992620E11Q,
452 -6.228685507402467503655405482985516909157E11Q,
453 -2.531972054436786351403749276956707260499E11Q,
454 -6.170200796658926701311867484296426831687E10Q,
455 -9.211477458528156048231908798456365081135E9Q,
456 -8.251806236175037114064561038908691305583E8Q,
457 -4.147886355917831049939930101151160447495E7Q,
458 -1.010851868928346082547075956946476932162E6Q,
459 -8.333374463411801009783402800801201603736E3Q
460 };
461 #define NRD3 9
462 static const __float128 RD3[NRD3 + 1] =
463 {
464 -5.216713843111675050627304523368029262450E11Q,
465 -8.014292925418308759369583419234079164391E11Q,
466 -5.180106858220030014546267824392678611990E11Q,
467 -1.830406975497439003897734969120997840011E11Q,
468 -3.845274631904879621945745960119924118925E10Q,
469 -4.891033385370523863288908070309417710903E9Q,
470 -3.670172254411328640353855768698287474282E8Q,
471 -1.505316381525727713026364396635522516989E7Q,
472 -2.856327162923716881454613540575964890347E5Q,
473 -1.622140448015769906847567212766206894547E3Q
474 /* 1.0E0Q */
475 };
476
477
478 /* log gamma(x+2.5) = log gamma(2.5) + x P(x)/Q(x)
479 -0.125 <= x <= 0.25
480 2.375 <= x+2.5 <= 2.75 */
481 static const __float128 lgam2r5a = 2.8466796875E-1Q;
482 static const __float128 lgam2r5b = 1.4901722919159632494669682701924320137696E-5Q;
483 #define NRN2r5 8
484 static const __float128 RN2r5[NRN2r5 + 1] =
485 {
486 -4.676454313888335499356699817678862233205E9Q,
487 -9.361888347911187924389905984624216340639E9Q,
488 -7.695353600835685037920815799526540237703E9Q,
489 -3.364370100981509060441853085968900734521E9Q,
490 -8.449902011848163568670361316804900559863E8Q,
491 -1.225249050950801905108001246436783022179E8Q,
492 -9.732972931077110161639900388121650470926E6Q,
493 -3.695711763932153505623248207576425983573E5Q,
494 -4.717341584067827676530426007495274711306E3Q
495 };
496 #define NRD2r5 8
497 static const __float128 RD2r5[NRD2r5 + 1] =
498 {
499 -6.650657966618993679456019224416926875619E9Q,
500 -1.099511409330635807899718829033488771623E10Q,
501 -7.482546968307837168164311101447116903148E9Q,
502 -2.702967190056506495988922973755870557217E9Q,
503 -5.570008176482922704972943389590409280950E8Q,
504 -6.536934032192792470926310043166993233231E7Q,
505 -4.101991193844953082400035444146067511725E6Q,
506 -1.174082735875715802334430481065526664020E5Q,
507 -9.932840389994157592102947657277692978511E2Q
508 /* 1.0E0Q */
509 };
510
511
512 /* log gamma(x+2) = x P(x)/Q(x)
513 -0.125 <= x <= +0.375
514 1.875 <= x+2 <= 2.375
515 Peak relative error 4.6e-36 */
516 #define NRN2 9
517 static const __float128 RN2[NRN2 + 1] =
518 {
519 -3.716661929737318153526921358113793421524E9Q,
520 -1.138816715030710406922819131397532331321E10Q,
521 -1.421017419363526524544402598734013569950E10Q,
522 -9.510432842542519665483662502132010331451E9Q,
523 -3.747528562099410197957514973274474767329E9Q,
524 -8.923565763363912474488712255317033616626E8Q,
525 -1.261396653700237624185350402781338231697E8Q,
526 -9.918402520255661797735331317081425749014E6Q,
527 -3.753996255897143855113273724233104768831E5Q,
528 -4.778761333044147141559311805999540765612E3Q
529 };
530 #define NRD2 9
531 static const __float128 RD2[NRD2 + 1] =
532 {
533 -8.790916836764308497770359421351673950111E9Q,
534 -2.023108608053212516399197678553737477486E10Q,
535 -1.958067901852022239294231785363504458367E10Q,
536 -1.035515043621003101254252481625188704529E10Q,
537 -3.253884432621336737640841276619272224476E9Q,
538 -6.186383531162456814954947669274235815544E8Q,
539 -6.932557847749518463038934953605969951466E7Q,
540 -4.240731768287359608773351626528479703758E6Q,
541 -1.197343995089189188078944689846348116630E5Q,
542 -1.004622911670588064824904487064114090920E3Q
543 /* 1.0E0 */
544 };
545
546
547 /* log gamma(x+1.75) = log gamma(1.75) + x P(x)/Q(x)
548 -0.125 <= x <= +0.125
549 1.625 <= x+1.75 <= 1.875
550 Peak relative error 9.2e-37 */
551 static const __float128 lgam1r75a = -8.441162109375E-2Q;
552 static const __float128 lgam1r75b = 1.0500073264444042213965868602268256157604E-5Q;
553 #define NRN1r75 8
554 static const __float128 RN1r75[NRN1r75 + 1] =
555 {
556 -5.221061693929833937710891646275798251513E7Q,
557 -2.052466337474314812817883030472496436993E8Q,
558 -2.952718275974940270675670705084125640069E8Q,
559 -2.132294039648116684922965964126389017840E8Q,
560 -8.554103077186505960591321962207519908489E7Q,
561 -1.940250901348870867323943119132071960050E7Q,
562 -2.379394147112756860769336400290402208435E6Q,
563 -1.384060879999526222029386539622255797389E5Q,
564 -2.698453601378319296159355612094598695530E3Q
565 };
566 #define NRD1r75 8
567 static const __float128 RD1r75[NRD1r75 + 1] =
568 {
569 -2.109754689501705828789976311354395393605E8Q,
570 -5.036651829232895725959911504899241062286E8Q,
571 -4.954234699418689764943486770327295098084E8Q,
572 -2.589558042412676610775157783898195339410E8Q,
573 -7.731476117252958268044969614034776883031E7Q,
574 -1.316721702252481296030801191240867486965E7Q,
575 -1.201296501404876774861190604303728810836E6Q,
576 -5.007966406976106636109459072523610273928E4Q,
577 -6.155817990560743422008969155276229018209E2Q
578 /* 1.0E0Q */
579 };
580
581
582 /* log gamma(x+x0) = y0 + x^2 P(x)/Q(x)
583 -0.0867 <= x <= +0.1634
584 1.374932... <= x+x0 <= 1.625032...
585 Peak relative error 4.0e-36 */
586 static const __float128 x0a = 1.4616241455078125Q;
587 static const __float128 x0b = 7.9994605498412626595423257213002588621246E-6Q;
588 static const __float128 y0a = -1.21490478515625E-1Q;
589 static const __float128 y0b = 4.1879797753919044854428223084178486438269E-6Q;
590 #define NRN1r5 8
591 static const __float128 RN1r5[NRN1r5 + 1] =
592 {
593 6.827103657233705798067415468881313128066E5Q,
594 1.910041815932269464714909706705242148108E6Q,
595 2.194344176925978377083808566251427771951E6Q,
596 1.332921400100891472195055269688876427962E6Q,
597 4.589080973377307211815655093824787123508E5Q,
598 8.900334161263456942727083580232613796141E4Q,
599 9.053840838306019753209127312097612455236E3Q,
600 4.053367147553353374151852319743594873771E2Q,
601 5.040631576303952022968949605613514584950E0Q
602 };
603 #define NRD1r5 8
604 static const __float128 RD1r5[NRD1r5 + 1] =
605 {
606 1.411036368843183477558773688484699813355E6Q,
607 4.378121767236251950226362443134306184849E6Q,
608 5.682322855631723455425929877581697918168E6Q,
609 3.999065731556977782435009349967042222375E6Q,
610 1.653651390456781293163585493620758410333E6Q,
611 4.067774359067489605179546964969435858311E5Q,
612 5.741463295366557346748361781768833633256E4Q,
613 4.226404539738182992856094681115746692030E3Q,
614 1.316980975410327975566999780608618774469E2Q,
615 /* 1.0E0Q */
616 };
617
618
619 /* log gamma(x+1.25) = log gamma(1.25) + x P(x)/Q(x)
620 -.125 <= x <= +.125
621 1.125 <= x+1.25 <= 1.375
622 Peak relative error = 4.9e-36 */
623 static const __float128 lgam1r25a = -9.82818603515625E-2Q;
624 static const __float128 lgam1r25b = 1.0023929749338536146197303364159774377296E-5Q;
625 #define NRN1r25 9
626 static const __float128 RN1r25[NRN1r25 + 1] =
627 {
628 -9.054787275312026472896002240379580536760E4Q,
629 -8.685076892989927640126560802094680794471E4Q,
630 2.797898965448019916967849727279076547109E5Q,
631 6.175520827134342734546868356396008898299E5Q,
632 5.179626599589134831538516906517372619641E5Q,
633 2.253076616239043944538380039205558242161E5Q,
634 5.312653119599957228630544772499197307195E4Q,
635 6.434329437514083776052669599834938898255E3Q,
636 3.385414416983114598582554037612347549220E2Q,
637 4.907821957946273805080625052510832015792E0Q
638 };
639 #define NRD1r25 8
640 static const __float128 RD1r25[NRD1r25 + 1] =
641 {
642 3.980939377333448005389084785896660309000E5Q,
643 1.429634893085231519692365775184490465542E6Q,
644 2.145438946455476062850151428438668234336E6Q,
645 1.743786661358280837020848127465970357893E6Q,
646 8.316364251289743923178092656080441655273E5Q,
647 2.355732939106812496699621491135458324294E5Q,
648 3.822267399625696880571810137601310855419E4Q,
649 3.228463206479133236028576845538387620856E3Q,
650 1.152133170470059555646301189220117965514E2Q
651 /* 1.0E0Q */
652 };
653
654
655 /* log gamma(x + 1) = x P(x)/Q(x)
656 0.0 <= x <= +0.125
657 1.0 <= x+1 <= 1.125
658 Peak relative error 1.1e-35 */
659 #define NRN1 8
660 static const __float128 RN1[NRN1 + 1] =
661 {
662 -9.987560186094800756471055681088744738818E3Q,
663 -2.506039379419574361949680225279376329742E4Q,
664 -1.386770737662176516403363873617457652991E4Q,
665 1.439445846078103202928677244188837130744E4Q,
666 2.159612048879650471489449668295139990693E4Q,
667 1.047439813638144485276023138173676047079E4Q,
668 2.250316398054332592560412486630769139961E3Q,
669 1.958510425467720733041971651126443864041E2Q,
670 4.516830313569454663374271993200291219855E0Q
671 };
672 #define NRD1 7
673 static const __float128 RD1[NRD1 + 1] =
674 {
675 1.730299573175751778863269333703788214547E4Q,
676 6.807080914851328611903744668028014678148E4Q,
677 1.090071629101496938655806063184092302439E5Q,
678 9.124354356415154289343303999616003884080E4Q,
679 4.262071638655772404431164427024003253954E4Q,
680 1.096981664067373953673982635805821283581E4Q,
681 1.431229503796575892151252708527595787588E3Q,
682 7.734110684303689320830401788262295992921E1Q
683 /* 1.0E0 */
684 };
685
686
687 /* log gamma(x + 1) = x P(x)/Q(x)
688 -0.125 <= x <= 0
689 0.875 <= x+1 <= 1.0
690 Peak relative error 7.0e-37 */
691 #define NRNr9 8
692 static const __float128 RNr9[NRNr9 + 1] =
693 {
694 4.441379198241760069548832023257571176884E5Q,
695 1.273072988367176540909122090089580368732E6Q,
696 9.732422305818501557502584486510048387724E5Q,
697 -5.040539994443998275271644292272870348684E5Q,
698 -1.208719055525609446357448132109723786736E6Q,
699 -7.434275365370936547146540554419058907156E5Q,
700 -2.075642969983377738209203358199008185741E5Q,
701 -2.565534860781128618589288075109372218042E4Q,
702 -1.032901669542994124131223797515913955938E3Q,
703 };
704 #define NRDr9 8
705 static const __float128 RDr9[NRDr9 + 1] =
706 {
707 -7.694488331323118759486182246005193998007E5Q,
708 -3.301918855321234414232308938454112213751E6Q,
709 -5.856830900232338906742924836032279404702E6Q,
710 -5.540672519616151584486240871424021377540E6Q,
711 -3.006530901041386626148342989181721176919E6Q,
712 -9.350378280513062139466966374330795935163E5Q,
713 -1.566179100031063346901755685375732739511E5Q,
714 -1.205016539620260779274902967231510804992E4Q,
715 -2.724583156305709733221564484006088794284E2Q
716 /* 1.0E0 */
717 };
718
719
720 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
721
722 static __float128
723 neval (__float128 x, const __float128 *p, int n)
724 {
725 __float128 y;
726
727 p += n;
728 y = *p--;
729 do
730 {
731 y = y * x + *p--;
732 }
733 while (--n > 0);
734 return y;
735 }
736
737
738 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
739
740 static __float128
741 deval (__float128 x, const __float128 *p, int n)
742 {
743 __float128 y;
744
745 p += n;
746 y = x + *p--;
747 do
748 {
749 y = y * x + *p--;
750 }
751 while (--n > 0);
752 return y;
753 }
754
755
756 __float128
757 lgammaq (__float128 x)
758 {
759 __float128 p, q, w, z, nx;
760 int i, nn, sign;
761
762 sign = 1;
763
764 if (! finiteq (x))
765 return x * x;
766
767 if (x == 0.0Q)
768 {
769 if (signbitq (x))
770 sign = -1;
771 }
772
773 if (x < 0.0Q)
774 {
775 q = -x;
776 p = floorq (q);
777 if (p == q)
778 return (one / (p - p));
779 i = p;
780 if ((i & 1) == 0)
781 sign = -1;
782 else
783 sign = 1;
784 z = q - p;
785 if (z > 0.5Q)
786 {
787 p += 1.0Q;
788 z = p - q;
789 }
790 z = q * sinq (PIQ * z);
791 if (z == 0.0Q)
792 return (sign * huge * huge);
793 w = lgammaq (q);
794 z = logq (PIQ / z) - w;
795 return (z);
796 }
797
798 if (x < 13.5Q)
799 {
800 p = 0.0Q;
801 nx = floorq (x + 0.5Q);
802 nn = nx;
803 switch (nn)
804 {
805 case 0:
806 /* log gamma (x + 1) = log(x) + log gamma(x) */
807 if (x <= 0.125)
808 {
809 p = x * neval (x, RN1, NRN1) / deval (x, RD1, NRD1);
810 }
811 else if (x <= 0.375)
812 {
813 z = x - 0.25Q;
814 p = z * neval (z, RN1r25, NRN1r25) / deval (z, RD1r25, NRD1r25);
815 p += lgam1r25b;
816 p += lgam1r25a;
817 }
818 else if (x <= 0.625)
819 {
820 z = x + (1.0Q - x0a);
821 z = z - x0b;
822 p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
823 p = p * z * z;
824 p = p + y0b;
825 p = p + y0a;
826 }
827 else if (x <= 0.875)
828 {
829 z = x - 0.75Q;
830 p = z * neval (z, RN1r75, NRN1r75) / deval (z, RD1r75, NRD1r75);
831 p += lgam1r75b;
832 p += lgam1r75a;
833 }
834 else
835 {
836 z = x - 1.0Q;
837 p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
838 }
839 p = p - logq (x);
840 break;
841
842 case 1:
843 if (x < 0.875Q)
844 {
845 if (x <= 0.625)
846 {
847 z = x + (1.0Q - x0a);
848 z = z - x0b;
849 p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
850 p = p * z * z;
851 p = p + y0b;
852 p = p + y0a;
853 }
854 else if (x <= 0.875)
855 {
856 z = x - 0.75Q;
857 p = z * neval (z, RN1r75, NRN1r75)
858 / deval (z, RD1r75, NRD1r75);
859 p += lgam1r75b;
860 p += lgam1r75a;
861 }
862 else
863 {
864 z = x - 1.0Q;
865 p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
866 }
867 p = p - logq (x);
868 }
869 else if (x < 1.0Q)
870 {
871 z = x - 1.0Q;
872 p = z * neval (z, RNr9, NRNr9) / deval (z, RDr9, NRDr9);
873 }
874 else if (x == 1.0Q)
875 p = 0.0Q;
876 else if (x <= 1.125Q)
877 {
878 z = x - 1.0Q;
879 p = z * neval (z, RN1, NRN1) / deval (z, RD1, NRD1);
880 }
881 else if (x <= 1.375)
882 {
883 z = x - 1.25Q;
884 p = z * neval (z, RN1r25, NRN1r25) / deval (z, RD1r25, NRD1r25);
885 p += lgam1r25b;
886 p += lgam1r25a;
887 }
888 else
889 {
890 /* 1.375 <= x+x0 <= 1.625 */
891 z = x - x0a;
892 z = z - x0b;
893 p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
894 p = p * z * z;
895 p = p + y0b;
896 p = p + y0a;
897 }
898 break;
899
900 case 2:
901 if (x < 1.625Q)
902 {
903 z = x - x0a;
904 z = z - x0b;
905 p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
906 p = p * z * z;
907 p = p + y0b;
908 p = p + y0a;
909 }
910 else if (x < 1.875Q)
911 {
912 z = x - 1.75Q;
913 p = z * neval (z, RN1r75, NRN1r75) / deval (z, RD1r75, NRD1r75);
914 p += lgam1r75b;
915 p += lgam1r75a;
916 }
917 else if (x == 2.0Q)
918 p = 0.0Q;
919 else if (x < 2.375Q)
920 {
921 z = x - 2.0Q;
922 p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
923 }
924 else
925 {
926 z = x - 2.5Q;
927 p = z * neval (z, RN2r5, NRN2r5) / deval (z, RD2r5, NRD2r5);
928 p += lgam2r5b;
929 p += lgam2r5a;
930 }
931 break;
932
933 case 3:
934 if (x < 2.75)
935 {
936 z = x - 2.5Q;
937 p = z * neval (z, RN2r5, NRN2r5) / deval (z, RD2r5, NRD2r5);
938 p += lgam2r5b;
939 p += lgam2r5a;
940 }
941 else
942 {
943 z = x - 3.0Q;
944 p = z * neval (z, RN3, NRN3) / deval (z, RD3, NRD3);
945 p += lgam3b;
946 p += lgam3a;
947 }
948 break;
949
950 case 4:
951 z = x - 4.0Q;
952 p = z * neval (z, RN4, NRN4) / deval (z, RD4, NRD4);
953 p += lgam4b;
954 p += lgam4a;
955 break;
956
957 case 5:
958 z = x - 5.0Q;
959 p = z * neval (z, RN5, NRN5) / deval (z, RD5, NRD5);
960 p += lgam5b;
961 p += lgam5a;
962 break;
963
964 case 6:
965 z = x - 6.0Q;
966 p = z * neval (z, RN6, NRN6) / deval (z, RD6, NRD6);
967 p += lgam6b;
968 p += lgam6a;
969 break;
970
971 case 7:
972 z = x - 7.0Q;
973 p = z * neval (z, RN7, NRN7) / deval (z, RD7, NRD7);
974 p += lgam7b;
975 p += lgam7a;
976 break;
977
978 case 8:
979 z = x - 8.0Q;
980 p = z * neval (z, RN8, NRN8) / deval (z, RD8, NRD8);
981 p += lgam8b;
982 p += lgam8a;
983 break;
984
985 case 9:
986 z = x - 9.0Q;
987 p = z * neval (z, RN9, NRN9) / deval (z, RD9, NRD9);
988 p += lgam9b;
989 p += lgam9a;
990 break;
991
992 case 10:
993 z = x - 10.0Q;
994 p = z * neval (z, RN10, NRN10) / deval (z, RD10, NRD10);
995 p += lgam10b;
996 p += lgam10a;
997 break;
998
999 case 11:
1000 z = x - 11.0Q;
1001 p = z * neval (z, RN11, NRN11) / deval (z, RD11, NRD11);
1002 p += lgam11b;
1003 p += lgam11a;
1004 break;
1005
1006 case 12:
1007 z = x - 12.0Q;
1008 p = z * neval (z, RN12, NRN12) / deval (z, RD12, NRD12);
1009 p += lgam12b;
1010 p += lgam12a;
1011 break;
1012
1013 case 13:
1014 z = x - 13.0Q;
1015 p = z * neval (z, RN13, NRN13) / deval (z, RD13, NRD13);
1016 p += lgam13b;
1017 p += lgam13a;
1018 break;
1019 }
1020 return p;
1021 }
1022
1023 if (x > MAXLGM)
1024 return (sign * huge * huge);
1025
1026 q = ls2pi - x;
1027 q = (x - 0.5Q) * logq (x) + q;
1028 if (x > 1.0e18Q)
1029 return (q);
1030
1031 p = 1.0Q / (x * x);
1032 q += neval (p, RASY, NRASY) / x;
1033 return (q);
1034 }