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
lgammaq(__float128 x)76 lgammaq (__float128 x)
77 {
78 #ifndef HAVE_MATH_H_SIGNGAM
79 int signgam;
80 #endif
81 return __quadmath_lgammaq_r (x, &signgam);
82 }
83
84 static const __float128 PIL = 3.1415926535897932384626433832795028841972E0Q;
85 static const __float128 MAXLGM = 1.0485738685148938358098967157129705071571E4928Q;
86 static const __float128 one = 1;
87 static 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 */
92 static const __float128 ls2pi = 9.1893853320467274178032973640561763986140E-1Q;
93 #define NRASY 12
94 static 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 */
116 static const __float128 lgam13a = 1.9987213134765625E1Q;
117 static const __float128 lgam13b = 1.3608962611495173623870550785125024484248E-6Q;
118 #define NRN13 7
119 static 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
131 static 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 */
148 static const __float128 lgam12a = 1.75023040771484375E1Q;
149 static const __float128 lgam12b = 3.7687254483392876529072161996717039575982E-6Q;
150 #define NRN12 7
151 static 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
163 static 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 */
180 static const __float128 lgam11a = 1.5104400634765625E1Q;
181 static const __float128 lgam11b = 1.1938309890295225709329251070371882250744E-5Q;
182 #define NRN11 7
183 static 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
195 static 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 */
212 static const __float128 lgam10a = 1.280181884765625E1Q;
213 static const __float128 lgam10b = 8.6324252196112077178745667061642811492557E-6Q;
214 #define NRN10 7
215 static 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
227 static 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 */
245 static const __float128 lgam9a = 1.06045989990234375E1Q;
246 static const __float128 lgam9b = 3.9037218127284172274007216547549861681400E-6Q;
247 #define NRN9 7
248 static 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
260 static 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 */
278 static const __float128 lgam8a = 8.525146484375E0Q;
279 static const __float128 lgam8b = 1.4876690414300165531036347125050759667737E-5Q;
280 #define NRN8 8
281 static 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
294 static 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 */
312 static const __float128 lgam7a = 6.5792388916015625E0Q;
313 static const __float128 lgam7b = 1.2320408538495060178292903945321122583007E-5Q;
314 #define NRN7 8
315 static 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
328 static 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 */
346 static const __float128 lgam6a = 4.7874908447265625E0Q;
347 static const __float128 lgam6b = 8.9805548349424770093452324304839959231517E-7Q;
348 #define NRN6 8
349 static 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
362 static 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 */
381 static const __float128 lgam5a = 3.17803955078125E0Q;
382 static const __float128 lgam5b = 1.4279566695619646941601297055408873990961E-5Q;
383 #define NRN5 9
384 static 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
398 static 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 */
417 static const __float128 lgam4a = 1.791748046875E0Q;
418 static const __float128 lgam4b = 1.1422353055000812477358380702272722990692E-5Q;
419 #define NRN4 9
420 static 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
434 static 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 */
454 static const __float128 lgam3a = 6.93145751953125E-1Q;
455 static const __float128 lgam3b = 1.4286068203094172321214581765680755001344E-6Q;
456
457 #define NRN3 9
458 static 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
472 static 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 */
491 static const __float128 lgam2r5a = 2.8466796875E-1Q;
492 static const __float128 lgam2r5b = 1.4901722919159632494669682701924320137696E-5Q;
493 #define NRN2r5 8
494 static 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
507 static 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
527 static 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
541 static 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 */
561 static const __float128 lgam1r75a = -8.441162109375E-2Q;
562 static const __float128 lgam1r75b = 1.0500073264444042213965868602268256157604E-5Q;
563 #define NRN1r75 8
564 static 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
577 static 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 */
596 static const __float128 x0a = 1.4616241455078125Q;
597 static const __float128 x0b = 7.9994605498412626595423257213002588621246E-6Q;
598 static const __float128 y0a = -1.21490478515625E-1Q;
599 static const __float128 y0b = 4.1879797753919044854428223084178486438269E-6Q;
600 #define NRN1r5 8
601 static 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
614 static 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 */
633 static const __float128 lgam1r25a = -9.82818603515625E-2Q;
634 static const __float128 lgam1r25b = 1.0023929749338536146197303364159774377296E-5Q;
635 #define NRN1r25 9
636 static 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
650 static 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
670 static 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
683 static 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
702 static 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
715 static 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
732 static __float128
neval(__float128 x,const __float128 * p,int n)733 neval (__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
750 static __float128
deval(__float128 x,const __float128 * p,int n)751 deval (__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
__quadmath_lgammaq_r(__float128 x,int * signgamp)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