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