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, write to the Free Software
  69     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
  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 #ifdef __STDC__
 758 long double
 759 __ieee754_lgammal_r (long double x, int *signgamp)
 760 #else
 761 long double
 762 __ieee754_lgammal_r (x, signgamp)
 763      long double x;
 764      int *signgamp;
 765 #endif
 766 {
 767   long double p, q, w, z, nx;
 768   int i, nn;
 769
 770   *signgamp = 1;
 771
 772   if (! __finitel (x))
 773     return x * x;
 774
 775   if (x == 0.0L)
 776     {
 777       if (__signbitl (x))
 778         *signgamp = -1;
 779     }
 780
 781   if (x < 0.0L)
 782     {
 783       q = -x;
 784       p = __floorl (q);
 785       if (p == q)
 786         return (one / (p - p));
 787       i = p;
 788       if ((i & 1) == 0)
 789         *signgamp = -1;
 790       else
 791         *signgamp = 1;
 792       z = q - p;
 793       if (z > 0.5L)
 794         {
 795           p += 1.0L;
 796           z = p - q;
 797         }
 798       z = q * __sinl (PIL * z);
 799       if (z == 0.0L)
 800         return (*signgamp * huge * huge);
 801       w = __ieee754_lgammal_r (q, &i);
 802       z = __logl (PIL / z) - w;
 803       return (z);
 804     }
 805
 806   if (x < 13.5L)
 807     {
 808       p = 0.0L;
 809       nx = __floorl (x + 0.5L);
 810       nn = nx;
 811       switch (nn)
 812         {
 813         case 0:
 814           /* log gamma (x + 1) = log(x) + log gamma(x) */
 815           if (x <= 0.125)
 816             {
 817               p = x * neval (x, RN1, NRN1) / deval (x, RD1, NRD1);
 818             }
 819           else if (x <= 0.375)
 820             {
 821               z = x - 0.25L;
 822               p = z * neval (z, RN1r25, NRN1r25) / deval (z, RD1r25, NRD1r25);
 823               p += lgam1r25b;
 824               p += lgam1r25a;
 825             }
 826           else if (x <= 0.625)
 827             {
 828               z = x + (1.0L - x0a);
 829               z = z - x0b;
 830               p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
 831               p = p * z * z;
 832               p = p + y0b;
 833               p = p + y0a;
 834             }
 835           else if (x <= 0.875)
 836             {
 837               z = x - 0.75L;
 838               p = z * neval (z, RN1r75, NRN1r75) / deval (z, RD1r75, NRD1r75);
 839               p += lgam1r75b;
 840               p += lgam1r75a;
 841             }
 842           else
 843             {
 844               z = x - 1.0L;
 845               p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
 846             }
 847           p = p - __logl (x);
 848           break;
 849
 850         case 1:
 851           if (x < 0.875L)
 852             {
 853               if (x <= 0.625)
 854                 {
 855                   z = x + (1.0L - x0a);
 856                   z = z - x0b;
 857                   p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
 858                   p = p * z * z;
 859                   p = p + y0b;
 860                   p = p + y0a;
 861                 }
 862               else if (x <= 0.875)
 863                 {
 864                   z = x - 0.75L;
 865                   p = z * neval (z, RN1r75, NRN1r75)
 866                         / deval (z, RD1r75, NRD1r75);
 867                   p += lgam1r75b;
 868                   p += lgam1r75a;
 869                 }
 870               else
 871                 {
 872                   z = x - 1.0L;
 873                   p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
 874                 }
 875               p = p - __logl (x);
 876             }
 877           else if (x < 1.0L)
 878             {
 879               z = x - 1.0L;
 880               p = z * neval (z, RNr9, NRNr9) / deval (z, RDr9, NRDr9);
 881             }
 882           else if (x == 1.0L)
 883             p = 0.0L;
 884           else if (x <= 1.125L)
 885             {
 886               z = x - 1.0L;
 887               p = z * neval (z, RN1, NRN1) / deval (z, RD1, NRD1);
 888             }
 889           else if (x <= 1.375)
 890             {
 891               z = x - 1.25L;
 892               p = z * neval (z, RN1r25, NRN1r25) / deval (z, RD1r25, NRD1r25);
 893               p += lgam1r25b;
 894               p += lgam1r25a;
 895             }
 896           else
 897             {
 898               /* 1.375 <= x+x0 <= 1.625 */
 899               z = x - x0a;
 900               z = z - x0b;
 901               p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
 902               p = p * z * z;
 903               p = p + y0b;
 904               p = p + y0a;
 905             }
 906           break;
 907
 908         case 2:
 909           if (x < 1.625L)
 910             {
 911               z = x - x0a;
 912               z = z - x0b;
 913               p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
 914               p = p * z * z;
 915               p = p + y0b;
 916               p = p + y0a;
 917             }
 918           else if (x < 1.875L)
 919             {
 920               z = x - 1.75L;
 921               p = z * neval (z, RN1r75, NRN1r75) / deval (z, RD1r75, NRD1r75);
 922               p += lgam1r75b;
 923               p += lgam1r75a;
 924             }
 925           else if (x == 2.0L)
 926             p = 0.0L;
 927           else if (x < 2.375L)
 928             {
 929               z = x - 2.0L;
 930               p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
 931             }
 932           else
 933             {
 934               z = x - 2.5L;
 935               p = z * neval (z, RN2r5, NRN2r5) / deval (z, RD2r5, NRD2r5);
 936               p += lgam2r5b;
 937               p += lgam2r5a;
 938             }
 939           break;
 940
 941         case 3:
 942           if (x < 2.75)
 943             {
 944               z = x - 2.5L;
 945               p = z * neval (z, RN2r5, NRN2r5) / deval (z, RD2r5, NRD2r5);
 946               p += lgam2r5b;
 947               p += lgam2r5a;
 948             }
 949           else
 950             {
 951               z = x - 3.0L;
 952               p = z * neval (z, RN3, NRN3) / deval (z, RD3, NRD3);
 953               p += lgam3b;
 954               p += lgam3a;
 955             }
 956           break;
 957
 958         case 4:
 959           z = x - 4.0L;
 960           p = z * neval (z, RN4, NRN4) / deval (z, RD4, NRD4);
 961           p += lgam4b;
 962           p += lgam4a;
 963           break;
 964
 965         case 5:
 966           z = x - 5.0L;
 967           p = z * neval (z, RN5, NRN5) / deval (z, RD5, NRD5);
 968           p += lgam5b;
 969           p += lgam5a;
 970           break;
 971
 972         case 6:
 973           z = x - 6.0L;
 974           p = z * neval (z, RN6, NRN6) / deval (z, RD6, NRD6);
 975           p += lgam6b;
 976           p += lgam6a;
 977           break;
 978
 979         case 7:
 980           z = x - 7.0L;
 981           p = z * neval (z, RN7, NRN7) / deval (z, RD7, NRD7);
 982           p += lgam7b;
 983           p += lgam7a;
 984           break;
 985
 986         case 8:
 987           z = x - 8.0L;
 988           p = z * neval (z, RN8, NRN8) / deval (z, RD8, NRD8);
 989           p += lgam8b;
 990           p += lgam8a;
 991           break;
 992
 993         case 9:
 994           z = x - 9.0L;
 995           p = z * neval (z, RN9, NRN9) / deval (z, RD9, NRD9);
 996           p += lgam9b;
 997           p += lgam9a;
 998           break;
 999
1000         case 10:
1001           z = x - 10.0L;
1002           p = z * neval (z, RN10, NRN10) / deval (z, RD10, NRD10);
1003           p += lgam10b;
1004           p += lgam10a;
1005           break;
1006
1007         case 11:
1008           z = x - 11.0L;
1009           p = z * neval (z, RN11, NRN11) / deval (z, RD11, NRD11);
1010           p += lgam11b;
1011           p += lgam11a;
1012           break;
1013
1014         case 12:
1015           z = x - 12.0L;
1016           p = z * neval (z, RN12, NRN12) / deval (z, RD12, NRD12);
1017           p += lgam12b;
1018           p += lgam12a;
1019           break;
1020
1021         case 13:
1022           z = x - 13.0L;
1023           p = z * neval (z, RN13, NRN13) / deval (z, RD13, NRD13);
1024           p += lgam13b;
1025           p += lgam13a;
1026           break;
1027         }
1028       return p;
1029     }
1030
1031   if (x > MAXLGM)
1032     return (*signgamp * huge * huge);
1033
1034   q = ls2pi - x;
1035   q = (x - 0.5L) * __logl (x) + q;
1036   if (x > 1.0e18L)
1037     return (q);
1038
1039   p = 1.0L / (x * x);
1040   q += neval (p, RASY, NRASY) / x;
1041   return (q);
1042 }