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