libquadmath/math/lgammaq.c

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