1 /* Natural logarithm of gamma function. IBM Extended Precision version.
2 Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
4 This library is free software; you can redistribute it and/or
5 modify it under the terms of the GNU Lesser General Public
6 License as published by the Free Software Foundation; either
7 version 2.1 of the License, or (at your option) any later version.
9 This library is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 Lesser General Public License for more details.
14 You should have received a copy of the GNU Lesser General Public
15 License along with this library; if not, see
16 <https://www.gnu.org/licenses/>. */
18 /* This file was copied from sysdeps/ieee754/ldbl-128/e_lgammal_r.c. */
22 #include <math_private.h>
24 #include <libm-alias-finite.h>
26 static const long double PIL
= 3.1415926535897932384626433832795028841972E0L
;
27 static const long double MAXLGM
= 0x5.d53649e2d469dbc1f01e99fd66p
+1012L;
28 static const long double one
= 1;
29 static const long double huge
= LDBL_MAX
;
31 /* log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x P(1/x^2)
32 1/x <= 0.0741 (x >= 13.495...)
33 Peak relative error 1.5e-36 */
34 static const long double ls2pi
= 9.1893853320467274178032973640561763986140E-1L;
36 static const long double RASY
[NRASY
+ 1] =
38 8.333333333333333333333333333310437112111E-2L,
39 -2.777777777777777777777774789556228296902E-3L,
40 7.936507936507936507795933938448586499183E-4L,
41 -5.952380952380952041799269756378148574045E-4L,
42 8.417508417507928904209891117498524452523E-4L,
43 -1.917526917481263997778542329739806086290E-3L,
44 6.410256381217852504446848671499409919280E-3L,
45 -2.955064066900961649768101034477363301626E-2L,
46 1.796402955865634243663453415388336954675E-1L,
47 -1.391522089007758553455753477688592767741E0L
,
48 1.326130089598399157988112385013829305510E1L
,
49 -1.420412699593782497803472576479997819149E2L
,
50 1.218058922427762808938869872528846787020E3L
54 /* log gamma(x+13) = log gamma(13) + x P(x)/Q(x)
57 Peak relative error 1.1e-36 */
58 static const long double lgam13a
= 1.9987213134765625E1L
;
59 static const long double lgam13b
= 1.3608962611495173623870550785125024484248E-6L;
61 static const long double RN13
[NRN13
+ 1] =
63 8.591478354823578150238226576156275285700E11L
,
64 2.347931159756482741018258864137297157668E11L
,
65 2.555408396679352028680662433943000804616E10L
,
66 1.408581709264464345480765758902967123937E9L
,
67 4.126759849752613822953004114044451046321E7L
,
68 6.133298899622688505854211579222889943778E5L
,
69 3.929248056293651597987893340755876578072E3L
,
70 6.850783280018706668924952057996075215223E0L
73 static const long double RD13
[NRD13
+ 1] =
75 3.401225382297342302296607039352935541669E11L
,
76 8.756765276918037910363513243563234551784E10L
,
77 8.873913342866613213078554180987647243903E9L
,
78 4.483797255342763263361893016049310017973E8L
,
79 1.178186288833066430952276702931512870676E7L
,
80 1.519928623743264797939103740132278337476E5L
,
81 7.989298844938119228411117593338850892311E2L
86 /* log gamma(x+12) = log gamma(12) + x P(x)/Q(x)
89 Peak relative error 4.1e-36 */
90 static const long double lgam12a
= 1.75023040771484375E1L
;
91 static const long double lgam12b
= 3.7687254483392876529072161996717039575982E-6L;
93 static const long double RN12
[NRN12
+ 1] =
95 4.709859662695606986110997348630997559137E11L
,
96 1.398713878079497115037857470168777995230E11L
,
97 1.654654931821564315970930093932954900867E10L
,
98 9.916279414876676861193649489207282144036E8L
,
99 3.159604070526036074112008954113411389879E7L
,
100 5.109099197547205212294747623977502492861E5L
,
101 3.563054878276102790183396740969279826988E3L
,
102 6.769610657004672719224614163196946862747E0L
105 static const long double RD12
[NRD12
+ 1] =
107 1.928167007860968063912467318985802726613E11L
,
108 5.383198282277806237247492369072266389233E10L
,
109 5.915693215338294477444809323037871058363E9L
,
110 3.241438287570196713148310560147925781342E8L
,
111 9.236680081763754597872713592701048455890E6L
,
112 1.292246897881650919242713651166596478850E5L
,
113 7.366532445427159272584194816076600211171E2L
118 /* log gamma(x+11) = log gamma(11) + x P(x)/Q(x)
121 Peak relative error 1.8e-35 */
122 static const long double lgam11a
= 1.5104400634765625E1L
;
123 static const long double lgam11b
= 1.1938309890295225709329251070371882250744E-5L;
125 static const long double RN11
[NRN11
+ 1] =
127 2.446960438029415837384622675816736622795E11L
,
128 7.955444974446413315803799763901729640350E10L
,
129 1.030555327949159293591618473447420338444E10L
,
130 6.765022131195302709153994345470493334946E8L
,
131 2.361892792609204855279723576041468347494E7L
,
132 4.186623629779479136428005806072176490125E5L
,
133 3.202506022088912768601325534149383594049E3L
,
134 6.681356101133728289358838690666225691363E0L
137 static const long double RD11
[NRD11
+ 1] =
139 1.040483786179428590683912396379079477432E11L
,
140 3.172251138489229497223696648369823779729E10L
,
141 3.806961885984850433709295832245848084614E9L
,
142 2.278070344022934913730015420611609620171E8L
,
143 7.089478198662651683977290023829391596481E6L
,
144 1.083246385105903533237139380509590158658E5L
,
145 6.744420991491385145885727942219463243597E2L
150 /* log gamma(x+10) = log gamma(10) + x P(x)/Q(x)
153 Peak relative error 5.4e-37 */
154 static const long double lgam10a
= 1.280181884765625E1L
;
155 static const long double lgam10b
= 8.6324252196112077178745667061642811492557E-6L;
157 static const long double RN10
[NRN10
+ 1] =
159 -1.239059737177249934158597996648808363783E14L
,
160 -4.725899566371458992365624673357356908719E13L
,
161 -7.283906268647083312042059082837754850808E12L
,
162 -5.802855515464011422171165179767478794637E11L
,
163 -2.532349691157548788382820303182745897298E10L
,
164 -5.884260178023777312587193693477072061820E8L
,
165 -6.437774864512125749845840472131829114906E6L
,
166 -2.350975266781548931856017239843273049384E4L
169 static const long double RD10
[NRD10
+ 1] =
171 -5.502645997581822567468347817182347679552E13L
,
172 -1.970266640239849804162284805400136473801E13L
,
173 -2.819677689615038489384974042561531409392E12L
,
174 -2.056105863694742752589691183194061265094E11L
,
175 -8.053670086493258693186307810815819662078E9L
,
176 -1.632090155573373286153427982504851867131E8L
,
177 -1.483575879240631280658077826889223634921E6L
,
178 -4.002806669713232271615885826373550502510E3L
183 /* log gamma(x+9) = log gamma(9) + x P(x)/Q(x)
186 Peak relative error 3.6e-36 */
187 static const long double lgam9a
= 1.06045989990234375E1L
;
188 static const long double lgam9b
= 3.9037218127284172274007216547549861681400E-6L;
190 static const long double RN9
[NRN9
+ 1] =
192 -4.936332264202687973364500998984608306189E13L
,
193 -2.101372682623700967335206138517766274855E13L
,
194 -3.615893404644823888655732817505129444195E12L
,
195 -3.217104993800878891194322691860075472926E11L
,
196 -1.568465330337375725685439173603032921399E10L
,
197 -4.073317518162025744377629219101510217761E8L
,
198 -4.983232096406156139324846656819246974500E6L
,
199 -2.036280038903695980912289722995505277253E4L
202 static const long double RD9
[NRD9
+ 1] =
204 -2.306006080437656357167128541231915480393E13L
,
205 -9.183606842453274924895648863832233799950E12L
,
206 -1.461857965935942962087907301194381010380E12L
,
207 -1.185728254682789754150068652663124298303E11L
,
208 -5.166285094703468567389566085480783070037E9L
,
209 -1.164573656694603024184768200787835094317E8L
,
210 -1.177343939483908678474886454113163527909E6L
,
211 -3.529391059783109732159524500029157638736E3L
216 /* log gamma(x+8) = log gamma(8) + x P(x)/Q(x)
219 Peak relative error 2.4e-37 */
220 static const long double lgam8a
= 8.525146484375E0L
;
221 static const long double lgam8b
= 1.4876690414300165531036347125050759667737E-5L;
223 static const long double RN8
[NRN8
+ 1] =
225 6.600775438203423546565361176829139703289E11L
,
226 3.406361267593790705240802723914281025800E11L
,
227 7.222460928505293914746983300555538432830E10L
,
228 8.102984106025088123058747466840656458342E9L
,
229 5.157620015986282905232150979772409345927E8L
,
230 1.851445288272645829028129389609068641517E7L
,
231 3.489261702223124354745894067468953756656E5L
,
232 2.892095396706665774434217489775617756014E3L
,
233 6.596977510622195827183948478627058738034E0L
236 static const long double RD8
[NRD8
+ 1] =
238 3.274776546520735414638114828622673016920E11L
,
239 1.581811207929065544043963828487733970107E11L
,
240 3.108725655667825188135393076860104546416E10L
,
241 3.193055010502912617128480163681842165730E9L
,
242 1.830871482669835106357529710116211541839E8L
,
243 5.790862854275238129848491555068073485086E6L
,
244 9.305213264307921522842678835618803553589E4L
,
245 6.216974105861848386918949336819572333622E2L
250 /* log gamma(x+7) = log gamma(7) + x P(x)/Q(x)
253 Peak relative error 3.2e-36 */
254 static const long double lgam7a
= 6.5792388916015625E0L
;
255 static const long double lgam7b
= 1.2320408538495060178292903945321122583007E-5L;
257 static const long double RN7
[NRN7
+ 1] =
259 2.065019306969459407636744543358209942213E11L
,
260 1.226919919023736909889724951708796532847E11L
,
261 2.996157990374348596472241776917953749106E10L
,
262 3.873001919306801037344727168434909521030E9L
,
263 2.841575255593761593270885753992732145094E8L
,
264 1.176342515359431913664715324652399565551E7L
,
265 2.558097039684188723597519300356028511547E5L
,
266 2.448525238332609439023786244782810774702E3L
,
267 6.460280377802030953041566617300902020435E0L
270 static const long double RD7
[NRD7
+ 1] =
272 1.102646614598516998880874785339049304483E11L
,
273 6.099297512712715445879759589407189290040E10L
,
274 1.372898136289611312713283201112060238351E10L
,
275 1.615306270420293159907951633566635172343E9L
,
276 1.061114435798489135996614242842561967459E8L
,
277 3.845638971184305248268608902030718674691E6L
,
278 7.081730675423444975703917836972720495507E4L
,
279 5.423122582741398226693137276201344096370E2L
284 /* log gamma(x+6) = log gamma(6) + x P(x)/Q(x)
287 Peak relative error 6.2e-37 */
288 static const long double lgam6a
= 4.7874908447265625E0L
;
289 static const long double lgam6b
= 8.9805548349424770093452324304839959231517E-7L;
291 static const long double RN6
[NRN6
+ 1] =
293 -3.538412754670746879119162116819571823643E13L
,
294 -2.613432593406849155765698121483394257148E13L
,
295 -8.020670732770461579558867891923784753062E12L
,
296 -1.322227822931250045347591780332435433420E12L
,
297 -1.262809382777272476572558806855377129513E11L
,
298 -7.015006277027660872284922325741197022467E9L
,
299 -2.149320689089020841076532186783055727299E8L
,
300 -3.167210585700002703820077565539658995316E6L
,
301 -1.576834867378554185210279285358586385266E4L
304 static const long double RD6
[NRD6
+ 1] =
306 -2.073955870771283609792355579558899389085E13L
,
307 -1.421592856111673959642750863283919318175E13L
,
308 -4.012134994918353924219048850264207074949E12L
,
309 -6.013361045800992316498238470888523722431E11L
,
310 -5.145382510136622274784240527039643430628E10L
,
311 -2.510575820013409711678540476918249524123E9L
,
312 -6.564058379709759600836745035871373240904E7L
,
313 -7.861511116647120540275354855221373571536E5L
,
314 -2.821943442729620524365661338459579270561E3L
319 /* log gamma(x+5) = log gamma(5) + x P(x)/Q(x)
322 Peak relative error 3.4e-37 */
323 static const long double lgam5a
= 3.17803955078125E0L
;
324 static const long double lgam5b
= 1.4279566695619646941601297055408873990961E-5L;
326 static const long double RN5
[NRN5
+ 1] =
328 2.010952885441805899580403215533972172098E11L
,
329 1.916132681242540921354921906708215338584E11L
,
330 7.679102403710581712903937970163206882492E10L
,
331 1.680514903671382470108010973615268125169E10L
,
332 2.181011222911537259440775283277711588410E9L
,
333 1.705361119398837808244780667539728356096E8L
,
334 7.792391565652481864976147945997033946360E6L
,
335 1.910741381027985291688667214472560023819E5L
,
336 2.088138241893612679762260077783794329559E3L
,
337 6.330318119566998299106803922739066556550E0L
340 static const long double RD5
[NRD5
+ 1] =
342 1.335189758138651840605141370223112376176E11L
,
343 1.174130445739492885895466097516530211283E11L
,
344 4.308006619274572338118732154886328519910E10L
,
345 8.547402888692578655814445003283720677468E9L
,
346 9.934628078575618309542580800421370730906E8L
,
347 6.847107420092173812998096295422311820672E7L
,
348 2.698552646016599923609773122139463150403E6L
,
349 5.526516251532464176412113632726150253215E4L
,
350 4.772343321713697385780533022595450486932E2L
355 /* log gamma(x+4) = log gamma(4) + x P(x)/Q(x)
358 Peak relative error 6.7e-37 */
359 static const long double lgam4a
= 1.791748046875E0L
;
360 static const long double lgam4b
= 1.1422353055000812477358380702272722990692E-5L;
362 static const long double RN4
[NRN4
+ 1] =
364 -1.026583408246155508572442242188887829208E13L
,
365 -1.306476685384622809290193031208776258809E13L
,
366 -7.051088602207062164232806511992978915508E12L
,
367 -2.100849457735620004967624442027793656108E12L
,
368 -3.767473790774546963588549871673843260569E11L
,
369 -4.156387497364909963498394522336575984206E10L
,
370 -2.764021460668011732047778992419118757746E9L
,
371 -1.036617204107109779944986471142938641399E8L
,
372 -1.895730886640349026257780896972598305443E6L
,
373 -1.180509051468390914200720003907727988201E4L
376 static const long double RD4
[NRD4
+ 1] =
378 -8.172669122056002077809119378047536240889E12L
,
379 -9.477592426087986751343695251801814226960E12L
,
380 -4.629448850139318158743900253637212801682E12L
,
381 -1.237965465892012573255370078308035272942E12L
,
382 -1.971624313506929845158062177061297598956E11L
,
383 -1.905434843346570533229942397763361493610E10L
,
384 -1.089409357680461419743730978512856675984E9L
,
385 -3.416703082301143192939774401370222822430E7L
,
386 -4.981791914177103793218433195857635265295E5L
,
387 -2.192507743896742751483055798411231453733E3L
392 /* log gamma(x+3) = log gamma(3) + x P(x)/Q(x)
395 Peak relative error 6.0e-37 */
396 static const long double lgam3a
= 6.93145751953125E-1L;
397 static const long double lgam3b
= 1.4286068203094172321214581765680755001344E-6L;
400 static const long double RN3
[NRN3
+ 1] =
402 -4.813901815114776281494823863935820876670E11L
,
403 -8.425592975288250400493910291066881992620E11L
,
404 -6.228685507402467503655405482985516909157E11L
,
405 -2.531972054436786351403749276956707260499E11L
,
406 -6.170200796658926701311867484296426831687E10L
,
407 -9.211477458528156048231908798456365081135E9L
,
408 -8.251806236175037114064561038908691305583E8L
,
409 -4.147886355917831049939930101151160447495E7L
,
410 -1.010851868928346082547075956946476932162E6L
,
411 -8.333374463411801009783402800801201603736E3L
414 static const long double RD3
[NRD3
+ 1] =
416 -5.216713843111675050627304523368029262450E11L
,
417 -8.014292925418308759369583419234079164391E11L
,
418 -5.180106858220030014546267824392678611990E11L
,
419 -1.830406975497439003897734969120997840011E11L
,
420 -3.845274631904879621945745960119924118925E10L
,
421 -4.891033385370523863288908070309417710903E9L
,
422 -3.670172254411328640353855768698287474282E8L
,
423 -1.505316381525727713026364396635522516989E7L
,
424 -2.856327162923716881454613540575964890347E5L
,
425 -1.622140448015769906847567212766206894547E3L
430 /* log gamma(x+2.5) = log gamma(2.5) + x P(x)/Q(x)
432 2.375 <= x+2.5 <= 2.75 */
433 static const long double lgam2r5a
= 2.8466796875E-1L;
434 static const long double lgam2r5b
= 1.4901722919159632494669682701924320137696E-5L;
436 static const long double RN2r5
[NRN2r5
+ 1] =
438 -4.676454313888335499356699817678862233205E9L
,
439 -9.361888347911187924389905984624216340639E9L
,
440 -7.695353600835685037920815799526540237703E9L
,
441 -3.364370100981509060441853085968900734521E9L
,
442 -8.449902011848163568670361316804900559863E8L
,
443 -1.225249050950801905108001246436783022179E8L
,
444 -9.732972931077110161639900388121650470926E6L
,
445 -3.695711763932153505623248207576425983573E5L
,
446 -4.717341584067827676530426007495274711306E3L
449 static const long double RD2r5
[NRD2r5
+ 1] =
451 -6.650657966618993679456019224416926875619E9L
,
452 -1.099511409330635807899718829033488771623E10L
,
453 -7.482546968307837168164311101447116903148E9L
,
454 -2.702967190056506495988922973755870557217E9L
,
455 -5.570008176482922704972943389590409280950E8L
,
456 -6.536934032192792470926310043166993233231E7L
,
457 -4.101991193844953082400035444146067511725E6L
,
458 -1.174082735875715802334430481065526664020E5L
,
459 -9.932840389994157592102947657277692978511E2L
464 /* log gamma(x+2) = x P(x)/Q(x)
465 -0.125 <= x <= +0.375
466 1.875 <= x+2 <= 2.375
467 Peak relative error 4.6e-36 */
469 static const long double RN2
[NRN2
+ 1] =
471 -3.716661929737318153526921358113793421524E9L
,
472 -1.138816715030710406922819131397532331321E10L
,
473 -1.421017419363526524544402598734013569950E10L
,
474 -9.510432842542519665483662502132010331451E9L
,
475 -3.747528562099410197957514973274474767329E9L
,
476 -8.923565763363912474488712255317033616626E8L
,
477 -1.261396653700237624185350402781338231697E8L
,
478 -9.918402520255661797735331317081425749014E6L
,
479 -3.753996255897143855113273724233104768831E5L
,
480 -4.778761333044147141559311805999540765612E3L
483 static const long double RD2
[NRD2
+ 1] =
485 -8.790916836764308497770359421351673950111E9L
,
486 -2.023108608053212516399197678553737477486E10L
,
487 -1.958067901852022239294231785363504458367E10L
,
488 -1.035515043621003101254252481625188704529E10L
,
489 -3.253884432621336737640841276619272224476E9L
,
490 -6.186383531162456814954947669274235815544E8L
,
491 -6.932557847749518463038934953605969951466E7L
,
492 -4.240731768287359608773351626528479703758E6L
,
493 -1.197343995089189188078944689846348116630E5L
,
494 -1.004622911670588064824904487064114090920E3L
499 /* log gamma(x+1.75) = log gamma(1.75) + x P(x)/Q(x)
500 -0.125 <= x <= +0.125
501 1.625 <= x+1.75 <= 1.875
502 Peak relative error 9.2e-37 */
503 static const long double lgam1r75a
= -8.441162109375E-2L;
504 static const long double lgam1r75b
= 1.0500073264444042213965868602268256157604E-5L;
506 static const long double RN1r75
[NRN1r75
+ 1] =
508 -5.221061693929833937710891646275798251513E7L
,
509 -2.052466337474314812817883030472496436993E8L
,
510 -2.952718275974940270675670705084125640069E8L
,
511 -2.132294039648116684922965964126389017840E8L
,
512 -8.554103077186505960591321962207519908489E7L
,
513 -1.940250901348870867323943119132071960050E7L
,
514 -2.379394147112756860769336400290402208435E6L
,
515 -1.384060879999526222029386539622255797389E5L
,
516 -2.698453601378319296159355612094598695530E3L
519 static const long double RD1r75
[NRD1r75
+ 1] =
521 -2.109754689501705828789976311354395393605E8L
,
522 -5.036651829232895725959911504899241062286E8L
,
523 -4.954234699418689764943486770327295098084E8L
,
524 -2.589558042412676610775157783898195339410E8L
,
525 -7.731476117252958268044969614034776883031E7L
,
526 -1.316721702252481296030801191240867486965E7L
,
527 -1.201296501404876774861190604303728810836E6L
,
528 -5.007966406976106636109459072523610273928E4L
,
529 -6.155817990560743422008969155276229018209E2L
534 /* log gamma(x+x0) = y0 + x^2 P(x)/Q(x)
535 -0.0867 <= x <= +0.1634
536 1.374932... <= x+x0 <= 1.625032...
537 Peak relative error 4.0e-36 */
538 static const long double x0a
= 1.4616241455078125L;
539 static const long double x0b
= 7.9994605498412626595423257213002588621246E-6L;
540 static const long double y0a
= -1.21490478515625E-1L;
541 static const long double y0b
= 4.1879797753919044854428223084178486438269E-6L;
543 static const long double RN1r5
[NRN1r5
+ 1] =
545 6.827103657233705798067415468881313128066E5L
,
546 1.910041815932269464714909706705242148108E6L
,
547 2.194344176925978377083808566251427771951E6L
,
548 1.332921400100891472195055269688876427962E6L
,
549 4.589080973377307211815655093824787123508E5L
,
550 8.900334161263456942727083580232613796141E4L
,
551 9.053840838306019753209127312097612455236E3L
,
552 4.053367147553353374151852319743594873771E2L
,
553 5.040631576303952022968949605613514584950E0L
556 static const long double RD1r5
[NRD1r5
+ 1] =
558 1.411036368843183477558773688484699813355E6L
,
559 4.378121767236251950226362443134306184849E6L
,
560 5.682322855631723455425929877581697918168E6L
,
561 3.999065731556977782435009349967042222375E6L
,
562 1.653651390456781293163585493620758410333E6L
,
563 4.067774359067489605179546964969435858311E5L
,
564 5.741463295366557346748361781768833633256E4L
,
565 4.226404539738182992856094681115746692030E3L
,
566 1.316980975410327975566999780608618774469E2L
,
571 /* log gamma(x+1.25) = log gamma(1.25) + x P(x)/Q(x)
573 1.125 <= x+1.25 <= 1.375
574 Peak relative error = 4.9e-36 */
575 static const long double lgam1r25a
= -9.82818603515625E-2L;
576 static const long double lgam1r25b
= 1.0023929749338536146197303364159774377296E-5L;
578 static const long double RN1r25
[NRN1r25
+ 1] =
580 -9.054787275312026472896002240379580536760E4L
,
581 -8.685076892989927640126560802094680794471E4L
,
582 2.797898965448019916967849727279076547109E5L
,
583 6.175520827134342734546868356396008898299E5L
,
584 5.179626599589134831538516906517372619641E5L
,
585 2.253076616239043944538380039205558242161E5L
,
586 5.312653119599957228630544772499197307195E4L
,
587 6.434329437514083776052669599834938898255E3L
,
588 3.385414416983114598582554037612347549220E2L
,
589 4.907821957946273805080625052510832015792E0L
592 static const long double RD1r25
[NRD1r25
+ 1] =
594 3.980939377333448005389084785896660309000E5L
,
595 1.429634893085231519692365775184490465542E6L
,
596 2.145438946455476062850151428438668234336E6L
,
597 1.743786661358280837020848127465970357893E6L
,
598 8.316364251289743923178092656080441655273E5L
,
599 2.355732939106812496699621491135458324294E5L
,
600 3.822267399625696880571810137601310855419E4L
,
601 3.228463206479133236028576845538387620856E3L
,
602 1.152133170470059555646301189220117965514E2L
607 /* log gamma(x + 1) = x P(x)/Q(x)
610 Peak relative error 1.1e-35 */
612 static const long double RN1
[NRN1
+ 1] =
614 -9.987560186094800756471055681088744738818E3L
,
615 -2.506039379419574361949680225279376329742E4L
,
616 -1.386770737662176516403363873617457652991E4L
,
617 1.439445846078103202928677244188837130744E4L
,
618 2.159612048879650471489449668295139990693E4L
,
619 1.047439813638144485276023138173676047079E4L
,
620 2.250316398054332592560412486630769139961E3L
,
621 1.958510425467720733041971651126443864041E2L
,
622 4.516830313569454663374271993200291219855E0L
625 static const long double RD1
[NRD1
+ 1] =
627 1.730299573175751778863269333703788214547E4L
,
628 6.807080914851328611903744668028014678148E4L
,
629 1.090071629101496938655806063184092302439E5L
,
630 9.124354356415154289343303999616003884080E4L
,
631 4.262071638655772404431164427024003253954E4L
,
632 1.096981664067373953673982635805821283581E4L
,
633 1.431229503796575892151252708527595787588E3L
,
634 7.734110684303689320830401788262295992921E1L
639 /* log gamma(x + 1) = x P(x)/Q(x)
642 Peak relative error 7.0e-37 */
644 static const long double RNr9
[NRNr9
+ 1] =
646 4.441379198241760069548832023257571176884E5L
,
647 1.273072988367176540909122090089580368732E6L
,
648 9.732422305818501557502584486510048387724E5L
,
649 -5.040539994443998275271644292272870348684E5L
,
650 -1.208719055525609446357448132109723786736E6L
,
651 -7.434275365370936547146540554419058907156E5L
,
652 -2.075642969983377738209203358199008185741E5L
,
653 -2.565534860781128618589288075109372218042E4L
,
654 -1.032901669542994124131223797515913955938E3L
,
657 static const long double RDr9
[NRDr9
+ 1] =
659 -7.694488331323118759486182246005193998007E5L
,
660 -3.301918855321234414232308938454112213751E6L
,
661 -5.856830900232338906742924836032279404702E6L
,
662 -5.540672519616151584486240871424021377540E6L
,
663 -3.006530901041386626148342989181721176919E6L
,
664 -9.350378280513062139466966374330795935163E5L
,
665 -1.566179100031063346901755685375732739511E5L
,
666 -1.205016539620260779274902967231510804992E4L
,
667 -2.724583156305709733221564484006088794284E2L
672 /* Evaluate P[n] x^n + P[n-1] x^(n-1) + ... + P[0] */
675 neval (long double x
, const long double *p
, int n
)
690 /* Evaluate x^n+1 + P[n] x^(n) + P[n-1] x^(n-1) + ... + P[0] */
693 deval (long double x
, const long double *p
, int n
)
709 __ieee754_lgammal_r (long double x
, int *signgamp
)
711 long double p
, q
, w
, z
, nx
;
727 if (x
< -2 && x
> -48)
728 return __lgamma_negl (x
, signgamp
);
732 return (one
/ fabsl (p
- p
));
733 long double halfp
= p
* 0.5L;
734 if (halfp
== floorl (halfp
))
746 z
= q
* __sinl (PIL
* z
);
747 w
= __ieee754_lgammal_r (q
, &i
);
748 z
= __logl (PIL
/ z
) - w
;
755 nx
= floorl (x
+ 0.5L);
760 /* log gamma (x + 1) = log(x) + log gamma(x) */
765 p
= x
* neval (x
, RN1
, NRN1
) / deval (x
, RD1
, NRD1
);
770 p
= z
* neval (z
, RN1r25
, NRN1r25
) / deval (z
, RD1r25
, NRD1r25
);
778 p
= neval (z
, RN1r5
, NRN1r5
) / deval (z
, RD1r5
, NRD1r5
);
786 p
= z
* neval (z
, RN1r75
, NRN1r75
) / deval (z
, RD1r75
, NRD1r75
);
793 p
= z
* neval (z
, RN2
, NRN2
) / deval (z
, RD2
, NRD2
);
805 p
= neval (z
, RN1r5
, NRN1r5
) / deval (z
, RD1r5
, NRD1r5
);
813 p
= z
* neval (z
, RN1r75
, NRN1r75
)
814 / deval (z
, RD1r75
, NRD1r75
);
821 p
= z
* neval (z
, RN2
, NRN2
) / deval (z
, RD2
, NRD2
);
828 p
= z
* neval (z
, RNr9
, NRNr9
) / deval (z
, RDr9
, NRDr9
);
832 else if (x
<= 1.125L)
835 p
= z
* neval (z
, RN1
, NRN1
) / deval (z
, RD1
, NRD1
);
840 p
= z
* neval (z
, RN1r25
, NRN1r25
) / deval (z
, RD1r25
, NRD1r25
);
846 /* 1.375 <= x+x0 <= 1.625 */
849 p
= neval (z
, RN1r5
, NRN1r5
) / deval (z
, RD1r5
, NRD1r5
);
861 p
= neval (z
, RN1r5
, NRN1r5
) / deval (z
, RD1r5
, NRD1r5
);
869 p
= z
* neval (z
, RN1r75
, NRN1r75
) / deval (z
, RD1r75
, NRD1r75
);
878 p
= z
* neval (z
, RN2
, NRN2
) / deval (z
, RD2
, NRD2
);
883 p
= z
* neval (z
, RN2r5
, NRN2r5
) / deval (z
, RD2r5
, NRD2r5
);
893 p
= z
* neval (z
, RN2r5
, NRN2r5
) / deval (z
, RD2r5
, NRD2r5
);
900 p
= z
* neval (z
, RN3
, NRN3
) / deval (z
, RD3
, NRD3
);
908 p
= z
* neval (z
, RN4
, NRN4
) / deval (z
, RD4
, NRD4
);
915 p
= z
* neval (z
, RN5
, NRN5
) / deval (z
, RD5
, NRD5
);
922 p
= z
* neval (z
, RN6
, NRN6
) / deval (z
, RD6
, NRD6
);
929 p
= z
* neval (z
, RN7
, NRN7
) / deval (z
, RD7
, NRD7
);
936 p
= z
* neval (z
, RN8
, NRN8
) / deval (z
, RD8
, NRD8
);
943 p
= z
* neval (z
, RN9
, NRN9
) / deval (z
, RD9
, NRD9
);
950 p
= z
* neval (z
, RN10
, NRN10
) / deval (z
, RD10
, NRD10
);
957 p
= z
* neval (z
, RN11
, NRN11
) / deval (z
, RD11
, NRD11
);
964 p
= z
* neval (z
, RN12
, NRN12
) / deval (z
, RD12
, NRD12
);
971 p
= z
* neval (z
, RN13
, NRN13
) / deval (z
, RD13
, NRD13
);
980 return (*signgamp
* huge
* huge
);
983 return x
* (__logl (x
) - 1);
985 q
= (x
- 0.5L) * __logl (x
) + q
;
990 q
+= neval (p
, RASY
, NRASY
) / x
;
993 libm_alias_finite (__ieee754_lgammal_r
, __lgammal_r
)