libquadmath/math/lgammaq.c

   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
  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
 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
 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
 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 }