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
76 static const _Float128 PIL = L(3.1415926535897932384626433832795028841972E0);
77 static const _Float128 MAXLGM = L(1.0485738685148938358098967157129705071571E4928);
78 static const _Float128 one = 1;
79 static 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 */
84 static const _Float128 ls2pi = L(9.1893853320467274178032973640561763986140E-1);
85 #define NRASY 12
86 static 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 */
108 static const _Float128 lgam13a = L(1.9987213134765625E1);
109 static const _Float128 lgam13b = L(1.3608962611495173623870550785125024484248E-6);
110 #define NRN13 7
111 static 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
123 static 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 */
140 static const _Float128 lgam12a = L(1.75023040771484375E1);
141 static const _Float128 lgam12b = L(3.7687254483392876529072161996717039575982E-6);
142 #define NRN12 7
143 static 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
155 static 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 */
172 static const _Float128 lgam11a = L(1.5104400634765625E1);
173 static const _Float128 lgam11b = L(1.1938309890295225709329251070371882250744E-5);
174 #define NRN11 7
175 static 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
187 static 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 */
204 static const _Float128 lgam10a = L(1.280181884765625E1);
205 static const _Float128 lgam10b = L(8.6324252196112077178745667061642811492557E-6);
206 #define NRN10 7
207 static 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
219 static 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 */
237 static const _Float128 lgam9a = L(1.06045989990234375E1);
238 static const _Float128 lgam9b = L(3.9037218127284172274007216547549861681400E-6);
239 #define NRN9 7
240 static 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
252 static 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 */
270 static const _Float128 lgam8a = L(8.525146484375E0);
271 static const _Float128 lgam8b = L(1.4876690414300165531036347125050759667737E-5);
272 #define NRN8 8
273 static 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
286 static 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 */
304 static const _Float128 lgam7a = L(6.5792388916015625E0);
305 static const _Float128 lgam7b = L(1.2320408538495060178292903945321122583007E-5);
306 #define NRN7 8
307 static 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
320 static 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 */
338 static const _Float128 lgam6a = L(4.7874908447265625E0);
339 static const _Float128 lgam6b = L(8.9805548349424770093452324304839959231517E-7);
340 #define NRN6 8
341 static 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
354 static 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 */
373 static const _Float128 lgam5a = L(3.17803955078125E0);
374 static const _Float128 lgam5b = L(1.4279566695619646941601297055408873990961E-5);
375 #define NRN5 9
376 static 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
390 static 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 */
409 static const _Float128 lgam4a = L(1.791748046875E0);
410 static const _Float128 lgam4b = L(1.1422353055000812477358380702272722990692E-5);
411 #define NRN4 9
412 static 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
426 static 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 */
446 static const _Float128 lgam3a = L(6.93145751953125E-1);
447 static const _Float128 lgam3b = L(1.4286068203094172321214581765680755001344E-6);
448
449 #define NRN3 9
450 static 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
464 static 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 */
483 static const _Float128 lgam2r5a = L(2.8466796875E-1);
484 static const _Float128 lgam2r5b = L(1.4901722919159632494669682701924320137696E-5);
485 #define NRN2r5 8
486 static 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
499 static 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
519 static 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
533 static 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 */
553 static const _Float128 lgam1r75a = L(-8.441162109375E-2);
554 static const _Float128 lgam1r75b = L(1.0500073264444042213965868602268256157604E-5);
555 #define NRN1r75 8
556 static 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
569 static 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 */
588 static const _Float128 x0a = L(1.4616241455078125);
589 static const _Float128 x0b = L(7.9994605498412626595423257213002588621246E-6);
590 static const _Float128 y0a = L(-1.21490478515625E-1);
591 static const _Float128 y0b = L(4.1879797753919044854428223084178486438269E-6);
592 #define NRN1r5 8
593 static 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
606 static 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 */
625 static const _Float128 lgam1r25a = L(-9.82818603515625E-2);
626 static const _Float128 lgam1r25b = L(1.0023929749338536146197303364159774377296E-5);
627 #define NRN1r25 9
628 static 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
642 static 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
662 static 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
675 static 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
694 static 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
707 static 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
724 static _Float128
725 neval (_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
742 static _Float128
743 deval (_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 }
1043 libm_alias_finite (__ieee754_lgammal_r, __lgammal_r)