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