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