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