xref: /relibc/openlibm/ld128/e_lgammal_r.c (revision beb387c4b293b36e4c6a5bf2036887998454a7d6)
1 /*	$OpenBSD: e_lgammal.c,v 1.3 2011/07/09 05:29:06 martynas Exp $	*/
2 
3 /*
4  * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5  *
6  * Permission to use, copy, modify, and distribute this software for any
7  * purpose with or without fee is hereby granted, provided that the above
8  * copyright notice and this permission notice appear in all copies.
9  *
10  * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11  * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12  * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13  * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14  * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15  * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16  * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17  */
18 
19 /*                                                      lgammal_r
20  *
21  *      Natural logarithm of gamma function
22  *
23  *
24  *
25  * SYNOPSIS:
26  *
27  * long double x, y, lgammal_r();
28  * int signgam;
29  *
30  * y = lgammal_r(x, &signgam);
31  *
32  *
33  *
34  * DESCRIPTION:
35  *
36  * Returns the base e (2.718...) logarithm of the absolute
37  * value of the gamma function of the argument.
38  * The sign (+1 or -1) of the gamma function is returned through signgamp.
39  *
40  * The positive domain is partitioned into numerous segments for approximation.
41  * For x > 10,
42  *   log gamma(x) = (x - 0.5) log(x) - x + log sqrt(2 pi) + 1/x R(1/x^2)
43  * Near the minimum at x = x0 = 1.46... the approximation is
44  *   log gamma(x0 + z) = log gamma(x0) + z^2 P(z)/Q(z)
45  * for small z.
46  * Elsewhere between 0 and 10,
47  *   log gamma(n + z) = log gamma(n) + z P(z)/Q(z)
48  * for various selected n and small z.
49  *
50  * The cosecant reflection formula is employed for negative arguments.
51  *
52  *
53  *
54  * ACCURACY:
55  *
56  *
57  * arithmetic      domain        # trials     peak         rms
58  *                                            Relative error:
59  *    IEEE         10, 30         100000     3.9e-34     9.8e-35
60  *    IEEE          0, 10         100000     3.8e-34     5.3e-35
61  *                                            Absolute error:
62  *    IEEE         -10, 0         100000     8.0e-34     8.0e-35
63  *    IEEE         -30, -10       100000     4.4e-34     1.0e-34
64  *    IEEE        -100, 100       100000                 1.0e-34
65  *
66  * The absolute error criterion is the same as relative error
67  * when the function magnitude is greater than one but it is absolute
68  * when the magnitude is less than one.
69  *
70  */
71 
72 #include <openlibm.h>
73 
74 #include "math_private.h"
75 
76 static const long double PIL = 3.1415926535897932384626433832795028841972E0L;
77 static const long double MAXLGM = 1.0485738685148938358098967157129705071571E4928L;
78 static const long double one = 1.0L;
79 static const long double huge = 1.0e4000L;
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 long double ls2pi = 9.1893853320467274178032973640561763986140E-1L;
85 #define NRASY 12
86 static const long double RASY[NRASY + 1] =
87 {
88   8.333333333333333333333333333310437112111E-2L,
89  -2.777777777777777777777774789556228296902E-3L,
90   7.936507936507936507795933938448586499183E-4L,
91  -5.952380952380952041799269756378148574045E-4L,
92   8.417508417507928904209891117498524452523E-4L,
93  -1.917526917481263997778542329739806086290E-3L,
94   6.410256381217852504446848671499409919280E-3L,
95  -2.955064066900961649768101034477363301626E-2L,
96   1.796402955865634243663453415388336954675E-1L,
97  -1.391522089007758553455753477688592767741E0L,
98   1.326130089598399157988112385013829305510E1L,
99  -1.420412699593782497803472576479997819149E2L,
100   1.218058922427762808938869872528846787020E3L
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 long double lgam13a = 1.9987213134765625E1L;
109 static const long double lgam13b = 1.3608962611495173623870550785125024484248E-6L;
110 #define NRN13 7
111 static const long double RN13[NRN13 + 1] =
112 {
113   8.591478354823578150238226576156275285700E11L,
114   2.347931159756482741018258864137297157668E11L,
115   2.555408396679352028680662433943000804616E10L,
116   1.408581709264464345480765758902967123937E9L,
117   4.126759849752613822953004114044451046321E7L,
118   6.133298899622688505854211579222889943778E5L,
119   3.929248056293651597987893340755876578072E3L,
120   6.850783280018706668924952057996075215223E0L
121 };
122 #define NRD13 6
123 static const long double RD13[NRD13 + 1] =
124 {
125   3.401225382297342302296607039352935541669E11L,
126   8.756765276918037910363513243563234551784E10L,
127   8.873913342866613213078554180987647243903E9L,
128   4.483797255342763263361893016049310017973E8L,
129   1.178186288833066430952276702931512870676E7L,
130   1.519928623743264797939103740132278337476E5L,
131   7.989298844938119228411117593338850892311E2L
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 long double lgam12a = 1.75023040771484375E1L;
141 static const long double lgam12b = 3.7687254483392876529072161996717039575982E-6L;
142 #define NRN12 7
143 static const long double RN12[NRN12 + 1] =
144 {
145   4.709859662695606986110997348630997559137E11L,
146   1.398713878079497115037857470168777995230E11L,
147   1.654654931821564315970930093932954900867E10L,
148   9.916279414876676861193649489207282144036E8L,
149   3.159604070526036074112008954113411389879E7L,
150   5.109099197547205212294747623977502492861E5L,
151   3.563054878276102790183396740969279826988E3L,
152   6.769610657004672719224614163196946862747E0L
153 };
154 #define NRD12 6
155 static const long double RD12[NRD12 + 1] =
156 {
157   1.928167007860968063912467318985802726613E11L,
158   5.383198282277806237247492369072266389233E10L,
159   5.915693215338294477444809323037871058363E9L,
160   3.241438287570196713148310560147925781342E8L,
161   9.236680081763754597872713592701048455890E6L,
162   1.292246897881650919242713651166596478850E5L,
163   7.366532445427159272584194816076600211171E2L
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 long double lgam11a = 1.5104400634765625E1L;
173 static const long double lgam11b = 1.1938309890295225709329251070371882250744E-5L;
174 #define NRN11 7
175 static const long double RN11[NRN11 + 1] =
176 {
177   2.446960438029415837384622675816736622795E11L,
178   7.955444974446413315803799763901729640350E10L,
179   1.030555327949159293591618473447420338444E10L,
180   6.765022131195302709153994345470493334946E8L,
181   2.361892792609204855279723576041468347494E7L,
182   4.186623629779479136428005806072176490125E5L,
183   3.202506022088912768601325534149383594049E3L,
184   6.681356101133728289358838690666225691363E0L
185 };
186 #define NRD11 6
187 static const long double RD11[NRD11 + 1] =
188 {
189   1.040483786179428590683912396379079477432E11L,
190   3.172251138489229497223696648369823779729E10L,
191   3.806961885984850433709295832245848084614E9L,
192   2.278070344022934913730015420611609620171E8L,
193   7.089478198662651683977290023829391596481E6L,
194   1.083246385105903533237139380509590158658E5L,
195   6.744420991491385145885727942219463243597E2L
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 long double lgam10a = 1.280181884765625E1L;
205 static const long double lgam10b = 8.6324252196112077178745667061642811492557E-6L;
206 #define NRN10 7
207 static const long double RN10[NRN10 + 1] =
208 {
209   -1.239059737177249934158597996648808363783E14L,
210   -4.725899566371458992365624673357356908719E13L,
211   -7.283906268647083312042059082837754850808E12L,
212   -5.802855515464011422171165179767478794637E11L,
213   -2.532349691157548788382820303182745897298E10L,
214   -5.884260178023777312587193693477072061820E8L,
215   -6.437774864512125749845840472131829114906E6L,
216   -2.350975266781548931856017239843273049384E4L
217 };
218 #define NRD10 7
219 static const long double RD10[NRD10 + 1] =
220 {
221   -5.502645997581822567468347817182347679552E13L,
222   -1.970266640239849804162284805400136473801E13L,
223   -2.819677689615038489384974042561531409392E12L,
224   -2.056105863694742752589691183194061265094E11L,
225   -8.053670086493258693186307810815819662078E9L,
226   -1.632090155573373286153427982504851867131E8L,
227   -1.483575879240631280658077826889223634921E6L,
228   -4.002806669713232271615885826373550502510E3L
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 long double lgam9a = 1.06045989990234375E1L;
238 static const long double lgam9b = 3.9037218127284172274007216547549861681400E-6L;
239 #define NRN9 7
240 static const long double RN9[NRN9 + 1] =
241 {
242   -4.936332264202687973364500998984608306189E13L,
243   -2.101372682623700967335206138517766274855E13L,
244   -3.615893404644823888655732817505129444195E12L,
245   -3.217104993800878891194322691860075472926E11L,
246   -1.568465330337375725685439173603032921399E10L,
247   -4.073317518162025744377629219101510217761E8L,
248   -4.983232096406156139324846656819246974500E6L,
249   -2.036280038903695980912289722995505277253E4L
250 };
251 #define NRD9 7
252 static const long double RD9[NRD9 + 1] =
253 {
254   -2.306006080437656357167128541231915480393E13L,
255   -9.183606842453274924895648863832233799950E12L,
256   -1.461857965935942962087907301194381010380E12L,
257   -1.185728254682789754150068652663124298303E11L,
258   -5.166285094703468567389566085480783070037E9L,
259   -1.164573656694603024184768200787835094317E8L,
260   -1.177343939483908678474886454113163527909E6L,
261   -3.529391059783109732159524500029157638736E3L
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 long double lgam8a = 8.525146484375E0L;
271 static const long double lgam8b = 1.4876690414300165531036347125050759667737E-5L;
272 #define NRN8 8
273 static const long double RN8[NRN8 + 1] =
274 {
275   6.600775438203423546565361176829139703289E11L,
276   3.406361267593790705240802723914281025800E11L,
277   7.222460928505293914746983300555538432830E10L,
278   8.102984106025088123058747466840656458342E9L,
279   5.157620015986282905232150979772409345927E8L,
280   1.851445288272645829028129389609068641517E7L,
281   3.489261702223124354745894067468953756656E5L,
282   2.892095396706665774434217489775617756014E3L,
283   6.596977510622195827183948478627058738034E0L
284 };
285 #define NRD8 7
286 static const long double RD8[NRD8 + 1] =
287 {
288   3.274776546520735414638114828622673016920E11L,
289   1.581811207929065544043963828487733970107E11L,
290   3.108725655667825188135393076860104546416E10L,
291   3.193055010502912617128480163681842165730E9L,
292   1.830871482669835106357529710116211541839E8L,
293   5.790862854275238129848491555068073485086E6L,
294   9.305213264307921522842678835618803553589E4L,
295   6.216974105861848386918949336819572333622E2L
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 long double lgam7a = 6.5792388916015625E0L;
305 static const long double lgam7b = 1.2320408538495060178292903945321122583007E-5L;
306 #define NRN7 8
307 static const long double RN7[NRN7 + 1] =
308 {
309   2.065019306969459407636744543358209942213E11L,
310   1.226919919023736909889724951708796532847E11L,
311   2.996157990374348596472241776917953749106E10L,
312   3.873001919306801037344727168434909521030E9L,
313   2.841575255593761593270885753992732145094E8L,
314   1.176342515359431913664715324652399565551E7L,
315   2.558097039684188723597519300356028511547E5L,
316   2.448525238332609439023786244782810774702E3L,
317   6.460280377802030953041566617300902020435E0L
318 };
319 #define NRD7 7
320 static const long double RD7[NRD7 + 1] =
321 {
322   1.102646614598516998880874785339049304483E11L,
323   6.099297512712715445879759589407189290040E10L,
324   1.372898136289611312713283201112060238351E10L,
325   1.615306270420293159907951633566635172343E9L,
326   1.061114435798489135996614242842561967459E8L,
327   3.845638971184305248268608902030718674691E6L,
328   7.081730675423444975703917836972720495507E4L,
329   5.423122582741398226693137276201344096370E2L
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 long double lgam6a = 4.7874908447265625E0L;
339 static const long double lgam6b = 8.9805548349424770093452324304839959231517E-7L;
340 #define NRN6 8
341 static const long double RN6[NRN6 + 1] =
342 {
343   -3.538412754670746879119162116819571823643E13L,
344   -2.613432593406849155765698121483394257148E13L,
345   -8.020670732770461579558867891923784753062E12L,
346   -1.322227822931250045347591780332435433420E12L,
347   -1.262809382777272476572558806855377129513E11L,
348   -7.015006277027660872284922325741197022467E9L,
349   -2.149320689089020841076532186783055727299E8L,
350   -3.167210585700002703820077565539658995316E6L,
351   -1.576834867378554185210279285358586385266E4L
352 };
353 #define NRD6 8
354 static const long double RD6[NRD6 + 1] =
355 {
356   -2.073955870771283609792355579558899389085E13L,
357   -1.421592856111673959642750863283919318175E13L,
358   -4.012134994918353924219048850264207074949E12L,
359   -6.013361045800992316498238470888523722431E11L,
360   -5.145382510136622274784240527039643430628E10L,
361   -2.510575820013409711678540476918249524123E9L,
362   -6.564058379709759600836745035871373240904E7L,
363   -7.861511116647120540275354855221373571536E5L,
364   -2.821943442729620524365661338459579270561E3L
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 long double lgam5a = 3.17803955078125E0L;
374 static const long double lgam5b = 1.4279566695619646941601297055408873990961E-5L;
375 #define NRN5 9
376 static const long double RN5[NRN5 + 1] =
377 {
378   2.010952885441805899580403215533972172098E11L,
379   1.916132681242540921354921906708215338584E11L,
380   7.679102403710581712903937970163206882492E10L,
381   1.680514903671382470108010973615268125169E10L,
382   2.181011222911537259440775283277711588410E9L,
383   1.705361119398837808244780667539728356096E8L,
384   7.792391565652481864976147945997033946360E6L,
385   1.910741381027985291688667214472560023819E5L,
386   2.088138241893612679762260077783794329559E3L,
387   6.330318119566998299106803922739066556550E0L
388 };
389 #define NRD5 8
390 static const long double RD5[NRD5 + 1] =
391 {
392   1.335189758138651840605141370223112376176E11L,
393   1.174130445739492885895466097516530211283E11L,
394   4.308006619274572338118732154886328519910E10L,
395   8.547402888692578655814445003283720677468E9L,
396   9.934628078575618309542580800421370730906E8L,
397   6.847107420092173812998096295422311820672E7L,
398   2.698552646016599923609773122139463150403E6L,
399   5.526516251532464176412113632726150253215E4L,
400   4.772343321713697385780533022595450486932E2L
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 long double lgam4a = 1.791748046875E0L;
410 static const long double lgam4b = 1.1422353055000812477358380702272722990692E-5L;
411 #define NRN4 9
412 static const long double RN4[NRN4 + 1] =
413 {
414   -1.026583408246155508572442242188887829208E13L,
415   -1.306476685384622809290193031208776258809E13L,
416   -7.051088602207062164232806511992978915508E12L,
417   -2.100849457735620004967624442027793656108E12L,
418   -3.767473790774546963588549871673843260569E11L,
419   -4.156387497364909963498394522336575984206E10L,
420   -2.764021460668011732047778992419118757746E9L,
421   -1.036617204107109779944986471142938641399E8L,
422   -1.895730886640349026257780896972598305443E6L,
423   -1.180509051468390914200720003907727988201E4L
424 };
425 #define NRD4 9
426 static const long double RD4[NRD4 + 1] =
427 {
428   -8.172669122056002077809119378047536240889E12L,
429   -9.477592426087986751343695251801814226960E12L,
430   -4.629448850139318158743900253637212801682E12L,
431   -1.237965465892012573255370078308035272942E12L,
432   -1.971624313506929845158062177061297598956E11L,
433   -1.905434843346570533229942397763361493610E10L,
434   -1.089409357680461419743730978512856675984E9L,
435   -3.416703082301143192939774401370222822430E7L,
436   -4.981791914177103793218433195857635265295E5L,
437   -2.192507743896742751483055798411231453733E3L
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 long double lgam3a = 6.93145751953125E-1L;
447 static const long double lgam3b = 1.4286068203094172321214581765680755001344E-6L;
448 
449 #define NRN3 9
450 static const long double RN3[NRN3 + 1] =
451 {
452   -4.813901815114776281494823863935820876670E11L,
453   -8.425592975288250400493910291066881992620E11L,
454   -6.228685507402467503655405482985516909157E11L,
455   -2.531972054436786351403749276956707260499E11L,
456   -6.170200796658926701311867484296426831687E10L,
457   -9.211477458528156048231908798456365081135E9L,
458   -8.251806236175037114064561038908691305583E8L,
459   -4.147886355917831049939930101151160447495E7L,
460   -1.010851868928346082547075956946476932162E6L,
461   -8.333374463411801009783402800801201603736E3L
462 };
463 #define NRD3 9
464 static const long double RD3[NRD3 + 1] =
465 {
466   -5.216713843111675050627304523368029262450E11L,
467   -8.014292925418308759369583419234079164391E11L,
468   -5.180106858220030014546267824392678611990E11L,
469   -1.830406975497439003897734969120997840011E11L,
470   -3.845274631904879621945745960119924118925E10L,
471   -4.891033385370523863288908070309417710903E9L,
472   -3.670172254411328640353855768698287474282E8L,
473   -1.505316381525727713026364396635522516989E7L,
474   -2.856327162923716881454613540575964890347E5L,
475   -1.622140448015769906847567212766206894547E3L
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 long double lgam2r5a = 2.8466796875E-1L;
484 static const long double lgam2r5b = 1.4901722919159632494669682701924320137696E-5L;
485 #define NRN2r5 8
486 static const long double RN2r5[NRN2r5 + 1] =
487 {
488   -4.676454313888335499356699817678862233205E9L,
489   -9.361888347911187924389905984624216340639E9L,
490   -7.695353600835685037920815799526540237703E9L,
491   -3.364370100981509060441853085968900734521E9L,
492   -8.449902011848163568670361316804900559863E8L,
493   -1.225249050950801905108001246436783022179E8L,
494   -9.732972931077110161639900388121650470926E6L,
495   -3.695711763932153505623248207576425983573E5L,
496   -4.717341584067827676530426007495274711306E3L
497 };
498 #define NRD2r5 8
499 static const long double RD2r5[NRD2r5 + 1] =
500 {
501   -6.650657966618993679456019224416926875619E9L,
502   -1.099511409330635807899718829033488771623E10L,
503   -7.482546968307837168164311101447116903148E9L,
504   -2.702967190056506495988922973755870557217E9L,
505   -5.570008176482922704972943389590409280950E8L,
506   -6.536934032192792470926310043166993233231E7L,
507   -4.101991193844953082400035444146067511725E6L,
508   -1.174082735875715802334430481065526664020E5L,
509   -9.932840389994157592102947657277692978511E2L
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 long double RN2[NRN2 + 1] =
520 {
521   -3.716661929737318153526921358113793421524E9L,
522   -1.138816715030710406922819131397532331321E10L,
523   -1.421017419363526524544402598734013569950E10L,
524   -9.510432842542519665483662502132010331451E9L,
525   -3.747528562099410197957514973274474767329E9L,
526   -8.923565763363912474488712255317033616626E8L,
527   -1.261396653700237624185350402781338231697E8L,
528   -9.918402520255661797735331317081425749014E6L,
529   -3.753996255897143855113273724233104768831E5L,
530   -4.778761333044147141559311805999540765612E3L
531 };
532 #define NRD2 9
533 static const long double RD2[NRD2 + 1] =
534 {
535   -8.790916836764308497770359421351673950111E9L,
536   -2.023108608053212516399197678553737477486E10L,
537   -1.958067901852022239294231785363504458367E10L,
538   -1.035515043621003101254252481625188704529E10L,
539   -3.253884432621336737640841276619272224476E9L,
540   -6.186383531162456814954947669274235815544E8L,
541   -6.932557847749518463038934953605969951466E7L,
542   -4.240731768287359608773351626528479703758E6L,
543   -1.197343995089189188078944689846348116630E5L,
544   -1.004622911670588064824904487064114090920E3L
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 long double lgam1r75a = -8.441162109375E-2L;
554 static const long double lgam1r75b = 1.0500073264444042213965868602268256157604E-5L;
555 #define NRN1r75 8
556 static const long double RN1r75[NRN1r75 + 1] =
557 {
558   -5.221061693929833937710891646275798251513E7L,
559   -2.052466337474314812817883030472496436993E8L,
560   -2.952718275974940270675670705084125640069E8L,
561   -2.132294039648116684922965964126389017840E8L,
562   -8.554103077186505960591321962207519908489E7L,
563   -1.940250901348870867323943119132071960050E7L,
564   -2.379394147112756860769336400290402208435E6L,
565   -1.384060879999526222029386539622255797389E5L,
566   -2.698453601378319296159355612094598695530E3L
567 };
568 #define NRD1r75 8
569 static const long double RD1r75[NRD1r75 + 1] =
570 {
571   -2.109754689501705828789976311354395393605E8L,
572   -5.036651829232895725959911504899241062286E8L,
573   -4.954234699418689764943486770327295098084E8L,
574   -2.589558042412676610775157783898195339410E8L,
575   -7.731476117252958268044969614034776883031E7L,
576   -1.316721702252481296030801191240867486965E7L,
577   -1.201296501404876774861190604303728810836E6L,
578   -5.007966406976106636109459072523610273928E4L,
579   -6.155817990560743422008969155276229018209E2L
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 long double x0a = 1.4616241455078125L;
589 static const long double x0b = 7.9994605498412626595423257213002588621246E-6L;
590 static const long double y0a = -1.21490478515625E-1L;
591 static const long double y0b = 4.1879797753919044854428223084178486438269E-6L;
592 #define NRN1r5 8
593 static const long double RN1r5[NRN1r5 + 1] =
594 {
595   6.827103657233705798067415468881313128066E5L,
596   1.910041815932269464714909706705242148108E6L,
597   2.194344176925978377083808566251427771951E6L,
598   1.332921400100891472195055269688876427962E6L,
599   4.589080973377307211815655093824787123508E5L,
600   8.900334161263456942727083580232613796141E4L,
601   9.053840838306019753209127312097612455236E3L,
602   4.053367147553353374151852319743594873771E2L,
603   5.040631576303952022968949605613514584950E0L
604 };
605 #define NRD1r5 8
606 static const long double RD1r5[NRD1r5 + 1] =
607 {
608   1.411036368843183477558773688484699813355E6L,
609   4.378121767236251950226362443134306184849E6L,
610   5.682322855631723455425929877581697918168E6L,
611   3.999065731556977782435009349967042222375E6L,
612   1.653651390456781293163585493620758410333E6L,
613   4.067774359067489605179546964969435858311E5L,
614   5.741463295366557346748361781768833633256E4L,
615   4.226404539738182992856094681115746692030E3L,
616   1.316980975410327975566999780608618774469E2L,
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 long double lgam1r25a = -9.82818603515625E-2L;
626 static const long double lgam1r25b = 1.0023929749338536146197303364159774377296E-5L;
627 #define NRN1r25 9
628 static const long double RN1r25[NRN1r25 + 1] =
629 {
630   -9.054787275312026472896002240379580536760E4L,
631   -8.685076892989927640126560802094680794471E4L,
632   2.797898965448019916967849727279076547109E5L,
633   6.175520827134342734546868356396008898299E5L,
634   5.179626599589134831538516906517372619641E5L,
635   2.253076616239043944538380039205558242161E5L,
636   5.312653119599957228630544772499197307195E4L,
637   6.434329437514083776052669599834938898255E3L,
638   3.385414416983114598582554037612347549220E2L,
639   4.907821957946273805080625052510832015792E0L
640 };
641 #define NRD1r25 8
642 static const long double RD1r25[NRD1r25 + 1] =
643 {
644   3.980939377333448005389084785896660309000E5L,
645   1.429634893085231519692365775184490465542E6L,
646   2.145438946455476062850151428438668234336E6L,
647   1.743786661358280837020848127465970357893E6L,
648   8.316364251289743923178092656080441655273E5L,
649   2.355732939106812496699621491135458324294E5L,
650   3.822267399625696880571810137601310855419E4L,
651   3.228463206479133236028576845538387620856E3L,
652   1.152133170470059555646301189220117965514E2L
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 long double RN1[NRN1 + 1] =
663 {
664   -9.987560186094800756471055681088744738818E3L,
665   -2.506039379419574361949680225279376329742E4L,
666   -1.386770737662176516403363873617457652991E4L,
667   1.439445846078103202928677244188837130744E4L,
668   2.159612048879650471489449668295139990693E4L,
669   1.047439813638144485276023138173676047079E4L,
670   2.250316398054332592560412486630769139961E3L,
671   1.958510425467720733041971651126443864041E2L,
672   4.516830313569454663374271993200291219855E0L
673 };
674 #define NRD1 7
675 static const long double RD1[NRD1 + 1] =
676 {
677   1.730299573175751778863269333703788214547E4L,
678   6.807080914851328611903744668028014678148E4L,
679   1.090071629101496938655806063184092302439E5L,
680   9.124354356415154289343303999616003884080E4L,
681   4.262071638655772404431164427024003253954E4L,
682   1.096981664067373953673982635805821283581E4L,
683   1.431229503796575892151252708527595787588E3L,
684   7.734110684303689320830401788262295992921E1L
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 long double RNr9[NRNr9 + 1] =
695 {
696   4.441379198241760069548832023257571176884E5L,
697   1.273072988367176540909122090089580368732E6L,
698   9.732422305818501557502584486510048387724E5L,
699   -5.040539994443998275271644292272870348684E5L,
700   -1.208719055525609446357448132109723786736E6L,
701   -7.434275365370936547146540554419058907156E5L,
702   -2.075642969983377738209203358199008185741E5L,
703   -2.565534860781128618589288075109372218042E4L,
704   -1.032901669542994124131223797515913955938E3L,
705 };
706 #define NRDr9 8
707 static const long double RDr9[NRDr9 + 1] =
708 {
709   -7.694488331323118759486182246005193998007E5L,
710   -3.301918855321234414232308938454112213751E6L,
711   -5.856830900232338906742924836032279404702E6L,
712   -5.540672519616151584486240871424021377540E6L,
713   -3.006530901041386626148342989181721176919E6L,
714   -9.350378280513062139466966374330795935163E5L,
715   -1.566179100031063346901755685375732739511E5L,
716   -1.205016539620260779274902967231510804992E4L,
717   -2.724583156305709733221564484006088794284E2L
718 /* 1.0E0 */
719 };
720 
721 
722 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
723 
724 static long double
725 neval (long double x, const long double *p, int n)
726 {
727   long double 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 long double
743 deval (long double x, const long double *p, int n)
744 {
745   long double 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 long double
759 lgammal_r(long double x, int *signgamp)
760 {
761   long double p, q, w, z, nx;
762   int i, nn;
763 
764   *signgamp = 1;
765 
766   if (! finite (x))
767     return x * x;
768 
769   if (x == 0.0L)
770     {
771       if (signbit (x))
772 	*signgamp = -1;
773       return one / fabsl (x);
774     }
775 
776   if (x < 0.0L)
777     {
778       q = -x;
779       p = floorl (q);
780       if (p == q)
781 	return (one / (p - p));
782       i = p;
783       if ((i & 1) == 0)
784 	*signgamp = -1;
785       else
786 	*signgamp = 1;
787       z = q - p;
788       if (z > 0.5L)
789 	{
790 	  p += 1.0L;
791 	  z = p - q;
792 	}
793       z = q * sinl (PIL * z);
794       if (z == 0.0L)
795 	return (*signgamp * huge * huge);
796       w = lgammal (q);
797       z = logl (PIL / z) - w;
798       return (z);
799     }
800 
801   if (x < 13.5L)
802     {
803       p = 0.0L;
804       nx = floorl (x + 0.5L);
805       nn = nx;
806       switch (nn)
807 	{
808 	case 0:
809 	  /* log gamma (x + 1) = log(x) + log gamma(x) */
810 	  if (x <= 0.125)
811 	    {
812 	      p = x * neval (x, RN1, NRN1) / deval (x, RD1, NRD1);
813 	    }
814 	  else if (x <= 0.375)
815 	    {
816 	      z = x - 0.25L;
817 	      p = z * neval (z, RN1r25, NRN1r25) / deval (z, RD1r25, NRD1r25);
818 	      p += lgam1r25b;
819 	      p += lgam1r25a;
820 	    }
821 	  else if (x <= 0.625)
822 	    {
823 	      z = x + (1.0L - x0a);
824 	      z = z - x0b;
825 	      p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
826 	      p = p * z * z;
827 	      p = p + y0b;
828 	      p = p + y0a;
829 	    }
830 	  else if (x <= 0.875)
831 	    {
832 	      z = x - 0.75L;
833 	      p = z * neval (z, RN1r75, NRN1r75) / deval (z, RD1r75, NRD1r75);
834 	      p += lgam1r75b;
835 	      p += lgam1r75a;
836 	    }
837 	  else
838 	    {
839 	      z = x - 1.0L;
840 	      p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
841 	    }
842 	  p = p - logl (x);
843 	  break;
844 
845 	case 1:
846 	  if (x < 0.875L)
847 	    {
848 	      if (x <= 0.625)
849 		{
850 		  z = x + (1.0L - x0a);
851 		  z = z - x0b;
852 		  p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
853 		  p = p * z * z;
854 		  p = p + y0b;
855 		  p = p + y0a;
856 		}
857 	      else if (x <= 0.875)
858 		{
859 		  z = x - 0.75L;
860 		  p = z * neval (z, RN1r75, NRN1r75)
861 			/ deval (z, RD1r75, NRD1r75);
862 		  p += lgam1r75b;
863 		  p += lgam1r75a;
864 		}
865 	      else
866 		{
867 		  z = x - 1.0L;
868 		  p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
869 		}
870 	      p = p - logl (x);
871 	    }
872 	  else if (x < 1.0L)
873 	    {
874 	      z = x - 1.0L;
875 	      p = z * neval (z, RNr9, NRNr9) / deval (z, RDr9, NRDr9);
876 	    }
877 	  else if (x == 1.0L)
878 	    p = 0.0L;
879 	  else if (x <= 1.125L)
880 	    {
881 	      z = x - 1.0L;
882 	      p = z * neval (z, RN1, NRN1) / deval (z, RD1, NRD1);
883 	    }
884 	  else if (x <= 1.375)
885 	    {
886 	      z = x - 1.25L;
887 	      p = z * neval (z, RN1r25, NRN1r25) / deval (z, RD1r25, NRD1r25);
888 	      p += lgam1r25b;
889 	      p += lgam1r25a;
890 	    }
891 	  else
892 	    {
893 	      /* 1.375 <= x+x0 <= 1.625 */
894 	      z = x - x0a;
895 	      z = z - x0b;
896 	      p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
897 	      p = p * z * z;
898 	      p = p + y0b;
899 	      p = p + y0a;
900 	    }
901 	  break;
902 
903 	case 2:
904 	  if (x < 1.625L)
905 	    {
906 	      z = x - x0a;
907 	      z = z - x0b;
908 	      p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
909 	      p = p * z * z;
910 	      p = p + y0b;
911 	      p = p + y0a;
912 	    }
913 	  else if (x < 1.875L)
914 	    {
915 	      z = x - 1.75L;
916 	      p = z * neval (z, RN1r75, NRN1r75) / deval (z, RD1r75, NRD1r75);
917 	      p += lgam1r75b;
918 	      p += lgam1r75a;
919 	    }
920 	  else if (x == 2.0L)
921 	    p = 0.0L;
922 	  else if (x < 2.375L)
923 	    {
924 	      z = x - 2.0L;
925 	      p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
926 	    }
927 	  else
928 	    {
929 	      z = x - 2.5L;
930 	      p = z * neval (z, RN2r5, NRN2r5) / deval (z, RD2r5, NRD2r5);
931 	      p += lgam2r5b;
932 	      p += lgam2r5a;
933 	    }
934 	  break;
935 
936 	case 3:
937 	  if (x < 2.75)
938 	    {
939 	      z = x - 2.5L;
940 	      p = z * neval (z, RN2r5, NRN2r5) / deval (z, RD2r5, NRD2r5);
941 	      p += lgam2r5b;
942 	      p += lgam2r5a;
943 	    }
944 	  else
945 	    {
946 	      z = x - 3.0L;
947 	      p = z * neval (z, RN3, NRN3) / deval (z, RD3, NRD3);
948 	      p += lgam3b;
949 	      p += lgam3a;
950 	    }
951 	  break;
952 
953 	case 4:
954 	  z = x - 4.0L;
955 	  p = z * neval (z, RN4, NRN4) / deval (z, RD4, NRD4);
956 	  p += lgam4b;
957 	  p += lgam4a;
958 	  break;
959 
960 	case 5:
961 	  z = x - 5.0L;
962 	  p = z * neval (z, RN5, NRN5) / deval (z, RD5, NRD5);
963 	  p += lgam5b;
964 	  p += lgam5a;
965 	  break;
966 
967 	case 6:
968 	  z = x - 6.0L;
969 	  p = z * neval (z, RN6, NRN6) / deval (z, RD6, NRD6);
970 	  p += lgam6b;
971 	  p += lgam6a;
972 	  break;
973 
974 	case 7:
975 	  z = x - 7.0L;
976 	  p = z * neval (z, RN7, NRN7) / deval (z, RD7, NRD7);
977 	  p += lgam7b;
978 	  p += lgam7a;
979 	  break;
980 
981 	case 8:
982 	  z = x - 8.0L;
983 	  p = z * neval (z, RN8, NRN8) / deval (z, RD8, NRD8);
984 	  p += lgam8b;
985 	  p += lgam8a;
986 	  break;
987 
988 	case 9:
989 	  z = x - 9.0L;
990 	  p = z * neval (z, RN9, NRN9) / deval (z, RD9, NRD9);
991 	  p += lgam9b;
992 	  p += lgam9a;
993 	  break;
994 
995 	case 10:
996 	  z = x - 10.0L;
997 	  p = z * neval (z, RN10, NRN10) / deval (z, RD10, NRD10);
998 	  p += lgam10b;
999 	  p += lgam10a;
1000 	  break;
1001 
1002 	case 11:
1003 	  z = x - 11.0L;
1004 	  p = z * neval (z, RN11, NRN11) / deval (z, RD11, NRD11);
1005 	  p += lgam11b;
1006 	  p += lgam11a;
1007 	  break;
1008 
1009 	case 12:
1010 	  z = x - 12.0L;
1011 	  p = z * neval (z, RN12, NRN12) / deval (z, RD12, NRD12);
1012 	  p += lgam12b;
1013 	  p += lgam12a;
1014 	  break;
1015 
1016 	case 13:
1017 	  z = x - 13.0L;
1018 	  p = z * neval (z, RN13, NRN13) / deval (z, RD13, NRD13);
1019 	  p += lgam13b;
1020 	  p += lgam13a;
1021 	  break;
1022 	}
1023       return p;
1024     }
1025 
1026   if (x > MAXLGM)
1027     return (*signgamp * huge * huge);
1028 
1029   q = ls2pi - x;
1030   q = (x - 0.5L) * logl (x) + q;
1031   if (x > 1.0e18L)
1032     return (q);
1033 
1034   p = 1.0L / (x * x);
1035   q += neval (p, RASY, NRASY) / x;
1036   return (q);
1037 }
1038