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