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