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