Back to index

glibc  2.9
e_lgammal_r.c
Go to the documentation of this file.
00001 /*                                                      lgammal
00002  *
00003  *      Natural logarithm of gamma function
00004  *
00005  *
00006  *
00007  * SYNOPSIS:
00008  *
00009  * long double x, y, lgammal();
00010  * extern int sgngam;
00011  *
00012  * y = lgammal(x);
00013  *
00014  *
00015  *
00016  * DESCRIPTION:
00017  *
00018  * Returns the base e (2.718...) logarithm of the absolute
00019  * value of the gamma function of the argument.
00020  * The sign (+1 or -1) of the gamma function is returned in a
00021  * global (extern) variable named sgngam.
00022  *
00023  * The positive domain is partitioned into numerous segments for approximation.
00024  * For x > 10,
00025  *   log gamma(x) = (x - 0.5) log(x) - x + log sqrt(2 pi) + 1/x R(1/x^2)
00026  * Near the minimum at x = x0 = 1.46... the approximation is
00027  *   log gamma(x0 + z) = log gamma(x0) + z^2 P(z)/Q(z)
00028  * for small z.
00029  * Elsewhere between 0 and 10,
00030  *   log gamma(n + z) = log gamma(n) + z P(z)/Q(z)
00031  * for various selected n and small z.
00032  *
00033  * The cosecant reflection formula is employed for negative arguments.
00034  *
00035  *
00036  *
00037  * ACCURACY:
00038  *
00039  *
00040  * arithmetic      domain        # trials     peak         rms
00041  *                                            Relative error:
00042  *    IEEE         10, 30         100000     3.9e-34     9.8e-35
00043  *    IEEE          0, 10         100000     3.8e-34     5.3e-35
00044  *                                            Absolute error:
00045  *    IEEE         -10, 0         100000     8.0e-34     8.0e-35
00046  *    IEEE         -30, -10       100000     4.4e-34     1.0e-34
00047  *    IEEE        -100, 100       100000                 1.0e-34
00048  *
00049  * The absolute error criterion is the same as relative error
00050  * when the function magnitude is greater than one but it is absolute
00051  * when the magnitude is less than one.
00052  *
00053  */
00054 
00055 /* Copyright 2001 by Stephen L. Moshier <moshier@na-net.ornl.gov>
00056 
00057     This library is free software; you can redistribute it and/or
00058     modify it under the terms of the GNU Lesser General Public
00059     License as published by the Free Software Foundation; either
00060     version 2.1 of the License, or (at your option) any later version.
00061 
00062     This library is distributed in the hope that it will be useful,
00063     but WITHOUT ANY WARRANTY; without even the implied warranty of
00064     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00065     Lesser General Public License for more details.
00066 
00067     You should have received a copy of the GNU Lesser General Public
00068     License along with this library; if not, write to the Free Software
00069     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00070 
00071 #include "math.h"
00072 #include "math_private.h"
00073 
00074 static const long double PIL = 3.1415926535897932384626433832795028841972E0L;
00075 static const long double MAXLGM = 1.0485738685148938358098967157129705071571E4928L;
00076 static const long double one = 1.0L;
00077 static const long double zero = 0.0L;
00078 static const long double huge = 1.0e4000L;
00079 
00080 /* log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x P(1/x^2)
00081    1/x <= 0.0741 (x >= 13.495...)
00082    Peak relative error 1.5e-36  */
00083 static const long double ls2pi = 9.1893853320467274178032973640561763986140E-1L;
00084 #define NRASY 12
00085 static const long double RASY[NRASY + 1] =
00086 {
00087   8.333333333333333333333333333310437112111E-2L,
00088  -2.777777777777777777777774789556228296902E-3L,
00089   7.936507936507936507795933938448586499183E-4L,
00090  -5.952380952380952041799269756378148574045E-4L,
00091   8.417508417507928904209891117498524452523E-4L,
00092  -1.917526917481263997778542329739806086290E-3L,
00093   6.410256381217852504446848671499409919280E-3L,
00094  -2.955064066900961649768101034477363301626E-2L,
00095   1.796402955865634243663453415388336954675E-1L,
00096  -1.391522089007758553455753477688592767741E0L,
00097   1.326130089598399157988112385013829305510E1L,
00098  -1.420412699593782497803472576479997819149E2L,
00099   1.218058922427762808938869872528846787020E3L
00100 };
00101 
00102 
00103 /* log gamma(x+13) = log gamma(13) +  x P(x)/Q(x)
00104    -0.5 <= x <= 0.5
00105    12.5 <= x+13 <= 13.5
00106    Peak relative error 1.1e-36  */
00107 static const long double lgam13a = 1.9987213134765625E1L;
00108 static const long double lgam13b = 1.3608962611495173623870550785125024484248E-6L;
00109 #define NRN13 7
00110 static const long double RN13[NRN13 + 1] =
00111 {
00112   8.591478354823578150238226576156275285700E11L,
00113   2.347931159756482741018258864137297157668E11L,
00114   2.555408396679352028680662433943000804616E10L,
00115   1.408581709264464345480765758902967123937E9L,
00116   4.126759849752613822953004114044451046321E7L,
00117   6.133298899622688505854211579222889943778E5L,
00118   3.929248056293651597987893340755876578072E3L,
00119   6.850783280018706668924952057996075215223E0L
00120 };
00121 #define NRD13 6
00122 static const long double RD13[NRD13 + 1] =
00123 {
00124   3.401225382297342302296607039352935541669E11L,
00125   8.756765276918037910363513243563234551784E10L,
00126   8.873913342866613213078554180987647243903E9L,
00127   4.483797255342763263361893016049310017973E8L,
00128   1.178186288833066430952276702931512870676E7L,
00129   1.519928623743264797939103740132278337476E5L,
00130   7.989298844938119228411117593338850892311E2L
00131  /* 1.0E0L */
00132 };
00133 
00134 
00135 /* log gamma(x+12) = log gamma(12) +  x P(x)/Q(x)
00136    -0.5 <= x <= 0.5
00137    11.5 <= x+12 <= 12.5
00138    Peak relative error 4.1e-36  */
00139 static const long double lgam12a = 1.75023040771484375E1L;
00140 static const long double lgam12b = 3.7687254483392876529072161996717039575982E-6L;
00141 #define NRN12 7
00142 static const long double RN12[NRN12 + 1] =
00143 {
00144   4.709859662695606986110997348630997559137E11L,
00145   1.398713878079497115037857470168777995230E11L,
00146   1.654654931821564315970930093932954900867E10L,
00147   9.916279414876676861193649489207282144036E8L,
00148   3.159604070526036074112008954113411389879E7L,
00149   5.109099197547205212294747623977502492861E5L,
00150   3.563054878276102790183396740969279826988E3L,
00151   6.769610657004672719224614163196946862747E0L
00152 };
00153 #define NRD12 6
00154 static const long double RD12[NRD12 + 1] =
00155 {
00156   1.928167007860968063912467318985802726613E11L,
00157   5.383198282277806237247492369072266389233E10L,
00158   5.915693215338294477444809323037871058363E9L,
00159   3.241438287570196713148310560147925781342E8L,
00160   9.236680081763754597872713592701048455890E6L,
00161   1.292246897881650919242713651166596478850E5L,
00162   7.366532445427159272584194816076600211171E2L
00163  /* 1.0E0L */
00164 };
00165 
00166 
00167 /* log gamma(x+11) = log gamma(11) +  x P(x)/Q(x)
00168    -0.5 <= x <= 0.5
00169    10.5 <= x+11 <= 11.5
00170    Peak relative error 1.8e-35  */
00171 static const long double lgam11a = 1.5104400634765625E1L;
00172 static const long double lgam11b = 1.1938309890295225709329251070371882250744E-5L;
00173 #define NRN11 7
00174 static const long double RN11[NRN11 + 1] =
00175 {
00176   2.446960438029415837384622675816736622795E11L,
00177   7.955444974446413315803799763901729640350E10L,
00178   1.030555327949159293591618473447420338444E10L,
00179   6.765022131195302709153994345470493334946E8L,
00180   2.361892792609204855279723576041468347494E7L,
00181   4.186623629779479136428005806072176490125E5L,
00182   3.202506022088912768601325534149383594049E3L,
00183   6.681356101133728289358838690666225691363E0L
00184 };
00185 #define NRD11 6
00186 static const long double RD11[NRD11 + 1] =
00187 {
00188   1.040483786179428590683912396379079477432E11L,
00189   3.172251138489229497223696648369823779729E10L,
00190   3.806961885984850433709295832245848084614E9L,
00191   2.278070344022934913730015420611609620171E8L,
00192   7.089478198662651683977290023829391596481E6L,
00193   1.083246385105903533237139380509590158658E5L,
00194   6.744420991491385145885727942219463243597E2L
00195  /* 1.0E0L */
00196 };
00197 
00198 
00199 /* log gamma(x+10) = log gamma(10) +  x P(x)/Q(x)
00200    -0.5 <= x <= 0.5
00201    9.5 <= x+10 <= 10.5
00202    Peak relative error 5.4e-37  */
00203 static const long double lgam10a = 1.280181884765625E1L;
00204 static const long double lgam10b = 8.6324252196112077178745667061642811492557E-6L;
00205 #define NRN10 7
00206 static const long double RN10[NRN10 + 1] =
00207 {
00208   -1.239059737177249934158597996648808363783E14L,
00209   -4.725899566371458992365624673357356908719E13L,
00210   -7.283906268647083312042059082837754850808E12L,
00211   -5.802855515464011422171165179767478794637E11L,
00212   -2.532349691157548788382820303182745897298E10L,
00213   -5.884260178023777312587193693477072061820E8L,
00214   -6.437774864512125749845840472131829114906E6L,
00215   -2.350975266781548931856017239843273049384E4L
00216 };
00217 #define NRD10 7
00218 static const long double RD10[NRD10 + 1] =
00219 {
00220   -5.502645997581822567468347817182347679552E13L,
00221   -1.970266640239849804162284805400136473801E13L,
00222   -2.819677689615038489384974042561531409392E12L,
00223   -2.056105863694742752589691183194061265094E11L,
00224   -8.053670086493258693186307810815819662078E9L,
00225   -1.632090155573373286153427982504851867131E8L,
00226   -1.483575879240631280658077826889223634921E6L,
00227   -4.002806669713232271615885826373550502510E3L
00228  /* 1.0E0L */
00229 };
00230 
00231 
00232 /* log gamma(x+9) = log gamma(9) +  x P(x)/Q(x)
00233    -0.5 <= x <= 0.5
00234    8.5 <= x+9 <= 9.5
00235    Peak relative error 3.6e-36  */
00236 static const long double lgam9a = 1.06045989990234375E1L;
00237 static const long double lgam9b = 3.9037218127284172274007216547549861681400E-6L;
00238 #define NRN9 7
00239 static const long double RN9[NRN9 + 1] =
00240 {
00241   -4.936332264202687973364500998984608306189E13L,
00242   -2.101372682623700967335206138517766274855E13L,
00243   -3.615893404644823888655732817505129444195E12L,
00244   -3.217104993800878891194322691860075472926E11L,
00245   -1.568465330337375725685439173603032921399E10L,
00246   -4.073317518162025744377629219101510217761E8L,
00247   -4.983232096406156139324846656819246974500E6L,
00248   -2.036280038903695980912289722995505277253E4L
00249 };
00250 #define NRD9 7
00251 static const long double RD9[NRD9 + 1] =
00252 {
00253   -2.306006080437656357167128541231915480393E13L,
00254   -9.183606842453274924895648863832233799950E12L,
00255   -1.461857965935942962087907301194381010380E12L,
00256   -1.185728254682789754150068652663124298303E11L,
00257   -5.166285094703468567389566085480783070037E9L,
00258   -1.164573656694603024184768200787835094317E8L,
00259   -1.177343939483908678474886454113163527909E6L,
00260   -3.529391059783109732159524500029157638736E3L
00261   /* 1.0E0L */
00262 };
00263 
00264 
00265 /* log gamma(x+8) = log gamma(8) +  x P(x)/Q(x)
00266    -0.5 <= x <= 0.5
00267    7.5 <= x+8 <= 8.5
00268    Peak relative error 2.4e-37  */
00269 static const long double lgam8a = 8.525146484375E0L;
00270 static const long double lgam8b = 1.4876690414300165531036347125050759667737E-5L;
00271 #define NRN8 8
00272 static const long double RN8[NRN8 + 1] =
00273 {
00274   6.600775438203423546565361176829139703289E11L,
00275   3.406361267593790705240802723914281025800E11L,
00276   7.222460928505293914746983300555538432830E10L,
00277   8.102984106025088123058747466840656458342E9L,
00278   5.157620015986282905232150979772409345927E8L,
00279   1.851445288272645829028129389609068641517E7L,
00280   3.489261702223124354745894067468953756656E5L,
00281   2.892095396706665774434217489775617756014E3L,
00282   6.596977510622195827183948478627058738034E0L
00283 };
00284 #define NRD8 7
00285 static const long double RD8[NRD8 + 1] =
00286 {
00287   3.274776546520735414638114828622673016920E11L,
00288   1.581811207929065544043963828487733970107E11L,
00289   3.108725655667825188135393076860104546416E10L,
00290   3.193055010502912617128480163681842165730E9L,
00291   1.830871482669835106357529710116211541839E8L,
00292   5.790862854275238129848491555068073485086E6L,
00293   9.305213264307921522842678835618803553589E4L,
00294   6.216974105861848386918949336819572333622E2L
00295   /* 1.0E0L */
00296 };
00297 
00298 
00299 /* log gamma(x+7) = log gamma(7) +  x P(x)/Q(x)
00300    -0.5 <= x <= 0.5
00301    6.5 <= x+7 <= 7.5
00302    Peak relative error 3.2e-36  */
00303 static const long double lgam7a = 6.5792388916015625E0L;
00304 static const long double lgam7b = 1.2320408538495060178292903945321122583007E-5L;
00305 #define NRN7 8
00306 static const long double RN7[NRN7 + 1] =
00307 {
00308   2.065019306969459407636744543358209942213E11L,
00309   1.226919919023736909889724951708796532847E11L,
00310   2.996157990374348596472241776917953749106E10L,
00311   3.873001919306801037344727168434909521030E9L,
00312   2.841575255593761593270885753992732145094E8L,
00313   1.176342515359431913664715324652399565551E7L,
00314   2.558097039684188723597519300356028511547E5L,
00315   2.448525238332609439023786244782810774702E3L,
00316   6.460280377802030953041566617300902020435E0L
00317 };
00318 #define NRD7 7
00319 static const long double RD7[NRD7 + 1] =
00320 {
00321   1.102646614598516998880874785339049304483E11L,
00322   6.099297512712715445879759589407189290040E10L,
00323   1.372898136289611312713283201112060238351E10L,
00324   1.615306270420293159907951633566635172343E9L,
00325   1.061114435798489135996614242842561967459E8L,
00326   3.845638971184305248268608902030718674691E6L,
00327   7.081730675423444975703917836972720495507E4L,
00328   5.423122582741398226693137276201344096370E2L
00329   /* 1.0E0L */
00330 };
00331 
00332 
00333 /* log gamma(x+6) = log gamma(6) +  x P(x)/Q(x)
00334    -0.5 <= x <= 0.5
00335    5.5 <= x+6 <= 6.5
00336    Peak relative error 6.2e-37  */
00337 static const long double lgam6a = 4.7874908447265625E0L;
00338 static const long double lgam6b = 8.9805548349424770093452324304839959231517E-7L;
00339 #define NRN6 8
00340 static const long double RN6[NRN6 + 1] =
00341 {
00342   -3.538412754670746879119162116819571823643E13L,
00343   -2.613432593406849155765698121483394257148E13L,
00344   -8.020670732770461579558867891923784753062E12L,
00345   -1.322227822931250045347591780332435433420E12L,
00346   -1.262809382777272476572558806855377129513E11L,
00347   -7.015006277027660872284922325741197022467E9L,
00348   -2.149320689089020841076532186783055727299E8L,
00349   -3.167210585700002703820077565539658995316E6L,
00350   -1.576834867378554185210279285358586385266E4L
00351 };
00352 #define NRD6 8
00353 static const long double RD6[NRD6 + 1] =
00354 {
00355   -2.073955870771283609792355579558899389085E13L,
00356   -1.421592856111673959642750863283919318175E13L,
00357   -4.012134994918353924219048850264207074949E12L,
00358   -6.013361045800992316498238470888523722431E11L,
00359   -5.145382510136622274784240527039643430628E10L,
00360   -2.510575820013409711678540476918249524123E9L,
00361   -6.564058379709759600836745035871373240904E7L,
00362   -7.861511116647120540275354855221373571536E5L,
00363   -2.821943442729620524365661338459579270561E3L
00364   /* 1.0E0L */
00365 };
00366 
00367 
00368 /* log gamma(x+5) = log gamma(5) +  x P(x)/Q(x)
00369    -0.5 <= x <= 0.5
00370    4.5 <= x+5 <= 5.5
00371    Peak relative error 3.4e-37  */
00372 static const long double lgam5a = 3.17803955078125E0L;
00373 static const long double lgam5b = 1.4279566695619646941601297055408873990961E-5L;
00374 #define NRN5 9
00375 static const long double RN5[NRN5 + 1] =
00376 {
00377   2.010952885441805899580403215533972172098E11L,
00378   1.916132681242540921354921906708215338584E11L,
00379   7.679102403710581712903937970163206882492E10L,
00380   1.680514903671382470108010973615268125169E10L,
00381   2.181011222911537259440775283277711588410E9L,
00382   1.705361119398837808244780667539728356096E8L,
00383   7.792391565652481864976147945997033946360E6L,
00384   1.910741381027985291688667214472560023819E5L,
00385   2.088138241893612679762260077783794329559E3L,
00386   6.330318119566998299106803922739066556550E0L
00387 };
00388 #define NRD5 8
00389 static const long double RD5[NRD5 + 1] =
00390 {
00391   1.335189758138651840605141370223112376176E11L,
00392   1.174130445739492885895466097516530211283E11L,
00393   4.308006619274572338118732154886328519910E10L,
00394   8.547402888692578655814445003283720677468E9L,
00395   9.934628078575618309542580800421370730906E8L,
00396   6.847107420092173812998096295422311820672E7L,
00397   2.698552646016599923609773122139463150403E6L,
00398   5.526516251532464176412113632726150253215E4L,
00399   4.772343321713697385780533022595450486932E2L
00400   /* 1.0E0L */
00401 };
00402 
00403 
00404 /* log gamma(x+4) = log gamma(4) +  x P(x)/Q(x)
00405    -0.5 <= x <= 0.5
00406    3.5 <= x+4 <= 4.5
00407    Peak relative error 6.7e-37  */
00408 static const long double lgam4a = 1.791748046875E0L;
00409 static const long double lgam4b = 1.1422353055000812477358380702272722990692E-5L;
00410 #define NRN4 9
00411 static const long double RN4[NRN4 + 1] =
00412 {
00413   -1.026583408246155508572442242188887829208E13L,
00414   -1.306476685384622809290193031208776258809E13L,
00415   -7.051088602207062164232806511992978915508E12L,
00416   -2.100849457735620004967624442027793656108E12L,
00417   -3.767473790774546963588549871673843260569E11L,
00418   -4.156387497364909963498394522336575984206E10L,
00419   -2.764021460668011732047778992419118757746E9L,
00420   -1.036617204107109779944986471142938641399E8L,
00421   -1.895730886640349026257780896972598305443E6L,
00422   -1.180509051468390914200720003907727988201E4L
00423 };
00424 #define NRD4 9
00425 static const long double RD4[NRD4 + 1] =
00426 {
00427   -8.172669122056002077809119378047536240889E12L,
00428   -9.477592426087986751343695251801814226960E12L,
00429   -4.629448850139318158743900253637212801682E12L,
00430   -1.237965465892012573255370078308035272942E12L,
00431   -1.971624313506929845158062177061297598956E11L,
00432   -1.905434843346570533229942397763361493610E10L,
00433   -1.089409357680461419743730978512856675984E9L,
00434   -3.416703082301143192939774401370222822430E7L,
00435   -4.981791914177103793218433195857635265295E5L,
00436   -2.192507743896742751483055798411231453733E3L
00437   /* 1.0E0L */
00438 };
00439 
00440 
00441 /* log gamma(x+3) = log gamma(3) +  x P(x)/Q(x)
00442    -0.25 <= x <= 0.5
00443    2.75 <= x+3 <= 3.5
00444    Peak relative error 6.0e-37  */
00445 static const long double lgam3a = 6.93145751953125E-1L;
00446 static const long double lgam3b = 1.4286068203094172321214581765680755001344E-6L;
00447 
00448 #define NRN3 9
00449 static const long double RN3[NRN3 + 1] =
00450 {
00451   -4.813901815114776281494823863935820876670E11L,
00452   -8.425592975288250400493910291066881992620E11L,
00453   -6.228685507402467503655405482985516909157E11L,
00454   -2.531972054436786351403749276956707260499E11L,
00455   -6.170200796658926701311867484296426831687E10L,
00456   -9.211477458528156048231908798456365081135E9L,
00457   -8.251806236175037114064561038908691305583E8L,
00458   -4.147886355917831049939930101151160447495E7L,
00459   -1.010851868928346082547075956946476932162E6L,
00460   -8.333374463411801009783402800801201603736E3L
00461 };
00462 #define NRD3 9
00463 static const long double RD3[NRD3 + 1] =
00464 {
00465   -5.216713843111675050627304523368029262450E11L,
00466   -8.014292925418308759369583419234079164391E11L,
00467   -5.180106858220030014546267824392678611990E11L,
00468   -1.830406975497439003897734969120997840011E11L,
00469   -3.845274631904879621945745960119924118925E10L,
00470   -4.891033385370523863288908070309417710903E9L,
00471   -3.670172254411328640353855768698287474282E8L,
00472   -1.505316381525727713026364396635522516989E7L,
00473   -2.856327162923716881454613540575964890347E5L,
00474   -1.622140448015769906847567212766206894547E3L
00475   /* 1.0E0L */
00476 };
00477 
00478 
00479 /* log gamma(x+2.5) = log gamma(2.5) +  x P(x)/Q(x)
00480    -0.125 <= x <= 0.25
00481    2.375 <= x+2.5 <= 2.75  */
00482 static const long double lgam2r5a = 2.8466796875E-1L;
00483 static const long double lgam2r5b = 1.4901722919159632494669682701924320137696E-5L;
00484 #define NRN2r5 8
00485 static const long double RN2r5[NRN2r5 + 1] =
00486 {
00487   -4.676454313888335499356699817678862233205E9L,
00488   -9.361888347911187924389905984624216340639E9L,
00489   -7.695353600835685037920815799526540237703E9L,
00490   -3.364370100981509060441853085968900734521E9L,
00491   -8.449902011848163568670361316804900559863E8L,
00492   -1.225249050950801905108001246436783022179E8L,
00493   -9.732972931077110161639900388121650470926E6L,
00494   -3.695711763932153505623248207576425983573E5L,
00495   -4.717341584067827676530426007495274711306E3L
00496 };
00497 #define NRD2r5 8
00498 static const long double RD2r5[NRD2r5 + 1] =
00499 {
00500   -6.650657966618993679456019224416926875619E9L,
00501   -1.099511409330635807899718829033488771623E10L,
00502   -7.482546968307837168164311101447116903148E9L,
00503   -2.702967190056506495988922973755870557217E9L,
00504   -5.570008176482922704972943389590409280950E8L,
00505   -6.536934032192792470926310043166993233231E7L,
00506   -4.101991193844953082400035444146067511725E6L,
00507   -1.174082735875715802334430481065526664020E5L,
00508   -9.932840389994157592102947657277692978511E2L
00509   /* 1.0E0L */
00510 };
00511 
00512 
00513 /* log gamma(x+2) = x P(x)/Q(x)
00514    -0.125 <= x <= +0.375
00515    1.875 <= x+2 <= 2.375
00516    Peak relative error 4.6e-36  */
00517 #define NRN2 9
00518 static const long double RN2[NRN2 + 1] =
00519 {
00520   -3.716661929737318153526921358113793421524E9L,
00521   -1.138816715030710406922819131397532331321E10L,
00522   -1.421017419363526524544402598734013569950E10L,
00523   -9.510432842542519665483662502132010331451E9L,
00524   -3.747528562099410197957514973274474767329E9L,
00525   -8.923565763363912474488712255317033616626E8L,
00526   -1.261396653700237624185350402781338231697E8L,
00527   -9.918402520255661797735331317081425749014E6L,
00528   -3.753996255897143855113273724233104768831E5L,
00529   -4.778761333044147141559311805999540765612E3L
00530 };
00531 #define NRD2 9
00532 static const long double RD2[NRD2 + 1] =
00533 {
00534   -8.790916836764308497770359421351673950111E9L,
00535   -2.023108608053212516399197678553737477486E10L,
00536   -1.958067901852022239294231785363504458367E10L,
00537   -1.035515043621003101254252481625188704529E10L,
00538   -3.253884432621336737640841276619272224476E9L,
00539   -6.186383531162456814954947669274235815544E8L,
00540   -6.932557847749518463038934953605969951466E7L,
00541   -4.240731768287359608773351626528479703758E6L,
00542   -1.197343995089189188078944689846348116630E5L,
00543   -1.004622911670588064824904487064114090920E3L
00544 /* 1.0E0 */
00545 };
00546 
00547 
00548 /* log gamma(x+1.75) = log gamma(1.75) +  x P(x)/Q(x)
00549    -0.125 <= x <= +0.125
00550    1.625 <= x+1.75 <= 1.875
00551    Peak relative error 9.2e-37 */
00552 static const long double lgam1r75a = -8.441162109375E-2L;
00553 static const long double lgam1r75b = 1.0500073264444042213965868602268256157604E-5L;
00554 #define NRN1r75 8
00555 static const long double RN1r75[NRN1r75 + 1] =
00556 {
00557   -5.221061693929833937710891646275798251513E7L,
00558   -2.052466337474314812817883030472496436993E8L,
00559   -2.952718275974940270675670705084125640069E8L,
00560   -2.132294039648116684922965964126389017840E8L,
00561   -8.554103077186505960591321962207519908489E7L,
00562   -1.940250901348870867323943119132071960050E7L,
00563   -2.379394147112756860769336400290402208435E6L,
00564   -1.384060879999526222029386539622255797389E5L,
00565   -2.698453601378319296159355612094598695530E3L
00566 };
00567 #define NRD1r75 8
00568 static const long double RD1r75[NRD1r75 + 1] =
00569 {
00570   -2.109754689501705828789976311354395393605E8L,
00571   -5.036651829232895725959911504899241062286E8L,
00572   -4.954234699418689764943486770327295098084E8L,
00573   -2.589558042412676610775157783898195339410E8L,
00574   -7.731476117252958268044969614034776883031E7L,
00575   -1.316721702252481296030801191240867486965E7L,
00576   -1.201296501404876774861190604303728810836E6L,
00577   -5.007966406976106636109459072523610273928E4L,
00578   -6.155817990560743422008969155276229018209E2L
00579   /* 1.0E0L */
00580 };
00581 
00582 
00583 /* log gamma(x+x0) = y0 +  x^2 P(x)/Q(x)
00584    -0.0867 <= x <= +0.1634
00585    1.374932... <= x+x0 <= 1.625032...
00586    Peak relative error 4.0e-36  */
00587 static const long double x0a = 1.4616241455078125L;
00588 static const long double x0b = 7.9994605498412626595423257213002588621246E-6L;
00589 static const long double y0a = -1.21490478515625E-1L;
00590 static const long double y0b = 4.1879797753919044854428223084178486438269E-6L;
00591 #define NRN1r5 8
00592 static const long double RN1r5[NRN1r5 + 1] =
00593 {
00594   6.827103657233705798067415468881313128066E5L,
00595   1.910041815932269464714909706705242148108E6L,
00596   2.194344176925978377083808566251427771951E6L,
00597   1.332921400100891472195055269688876427962E6L,
00598   4.589080973377307211815655093824787123508E5L,
00599   8.900334161263456942727083580232613796141E4L,
00600   9.053840838306019753209127312097612455236E3L,
00601   4.053367147553353374151852319743594873771E2L,
00602   5.040631576303952022968949605613514584950E0L
00603 };
00604 #define NRD1r5 8
00605 static const long double RD1r5[NRD1r5 + 1] =
00606 {
00607   1.411036368843183477558773688484699813355E6L,
00608   4.378121767236251950226362443134306184849E6L,
00609   5.682322855631723455425929877581697918168E6L,
00610   3.999065731556977782435009349967042222375E6L,
00611   1.653651390456781293163585493620758410333E6L,
00612   4.067774359067489605179546964969435858311E5L,
00613   5.741463295366557346748361781768833633256E4L,
00614   4.226404539738182992856094681115746692030E3L,
00615   1.316980975410327975566999780608618774469E2L,
00616   /* 1.0E0L */
00617 };
00618 
00619 
00620 /* log gamma(x+1.25) = log gamma(1.25) +  x P(x)/Q(x)
00621    -.125 <= x <= +.125
00622    1.125 <= x+1.25 <= 1.375
00623    Peak relative error = 4.9e-36 */
00624 static const long double lgam1r25a = -9.82818603515625E-2L;
00625 static const long double lgam1r25b = 1.0023929749338536146197303364159774377296E-5L;
00626 #define NRN1r25 9
00627 static const long double RN1r25[NRN1r25 + 1] =
00628 {
00629   -9.054787275312026472896002240379580536760E4L,
00630   -8.685076892989927640126560802094680794471E4L,
00631   2.797898965448019916967849727279076547109E5L,
00632   6.175520827134342734546868356396008898299E5L,
00633   5.179626599589134831538516906517372619641E5L,
00634   2.253076616239043944538380039205558242161E5L,
00635   5.312653119599957228630544772499197307195E4L,
00636   6.434329437514083776052669599834938898255E3L,
00637   3.385414416983114598582554037612347549220E2L,
00638   4.907821957946273805080625052510832015792E0L
00639 };
00640 #define NRD1r25 8
00641 static const long double RD1r25[NRD1r25 + 1] =
00642 {
00643   3.980939377333448005389084785896660309000E5L,
00644   1.429634893085231519692365775184490465542E6L,
00645   2.145438946455476062850151428438668234336E6L,
00646   1.743786661358280837020848127465970357893E6L,
00647   8.316364251289743923178092656080441655273E5L,
00648   2.355732939106812496699621491135458324294E5L,
00649   3.822267399625696880571810137601310855419E4L,
00650   3.228463206479133236028576845538387620856E3L,
00651   1.152133170470059555646301189220117965514E2L
00652   /* 1.0E0L */
00653 };
00654 
00655 
00656 /* log gamma(x + 1) = x P(x)/Q(x)
00657    0.0 <= x <= +0.125
00658    1.0 <= x+1 <= 1.125
00659    Peak relative error 1.1e-35  */
00660 #define NRN1 8
00661 static const long double RN1[NRN1 + 1] =
00662 {
00663   -9.987560186094800756471055681088744738818E3L,
00664   -2.506039379419574361949680225279376329742E4L,
00665   -1.386770737662176516403363873617457652991E4L,
00666   1.439445846078103202928677244188837130744E4L,
00667   2.159612048879650471489449668295139990693E4L,
00668   1.047439813638144485276023138173676047079E4L,
00669   2.250316398054332592560412486630769139961E3L,
00670   1.958510425467720733041971651126443864041E2L,
00671   4.516830313569454663374271993200291219855E0L
00672 };
00673 #define NRD1 7
00674 static const long double RD1[NRD1 + 1] =
00675 {
00676   1.730299573175751778863269333703788214547E4L,
00677   6.807080914851328611903744668028014678148E4L,
00678   1.090071629101496938655806063184092302439E5L,
00679   9.124354356415154289343303999616003884080E4L,
00680   4.262071638655772404431164427024003253954E4L,
00681   1.096981664067373953673982635805821283581E4L,
00682   1.431229503796575892151252708527595787588E3L,
00683   7.734110684303689320830401788262295992921E1L
00684  /* 1.0E0 */
00685 };
00686 
00687 
00688 /* log gamma(x + 1) = x P(x)/Q(x)
00689    -0.125 <= x <= 0
00690    0.875 <= x+1 <= 1.0
00691    Peak relative error 7.0e-37  */
00692 #define NRNr9 8
00693 static const long double RNr9[NRNr9 + 1] =
00694 {
00695   4.441379198241760069548832023257571176884E5L,
00696   1.273072988367176540909122090089580368732E6L,
00697   9.732422305818501557502584486510048387724E5L,
00698   -5.040539994443998275271644292272870348684E5L,
00699   -1.208719055525609446357448132109723786736E6L,
00700   -7.434275365370936547146540554419058907156E5L,
00701   -2.075642969983377738209203358199008185741E5L,
00702   -2.565534860781128618589288075109372218042E4L,
00703   -1.032901669542994124131223797515913955938E3L,
00704 };
00705 #define NRDr9 8
00706 static const long double RDr9[NRDr9 + 1] =
00707 {
00708   -7.694488331323118759486182246005193998007E5L,
00709   -3.301918855321234414232308938454112213751E6L,
00710   -5.856830900232338906742924836032279404702E6L,
00711   -5.540672519616151584486240871424021377540E6L,
00712   -3.006530901041386626148342989181721176919E6L,
00713   -9.350378280513062139466966374330795935163E5L,
00714   -1.566179100031063346901755685375732739511E5L,
00715   -1.205016539620260779274902967231510804992E4L,
00716   -2.724583156305709733221564484006088794284E2L
00717 /* 1.0E0 */
00718 };
00719 
00720 
00721 /* Evaluate P[n] x^n  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
00722 
00723 static long double
00724 neval (long double x, const long double *p, int n)
00725 {
00726   long double y;
00727 
00728   p += n;
00729   y = *p--;
00730   do
00731     {
00732       y = y * x + *p--;
00733     }
00734   while (--n > 0);
00735   return y;
00736 }
00737 
00738 
00739 /* Evaluate x^n+1  +  P[n] x^(n)  +  P[n-1] x^(n-1)  +  ...  +  P[0] */
00740 
00741 static long double
00742 deval (long double x, const long double *p, int n)
00743 {
00744   long double y;
00745 
00746   p += n;
00747   y = x + *p--;
00748   do
00749     {
00750       y = y * x + *p--;
00751     }
00752   while (--n > 0);
00753   return y;
00754 }
00755 
00756 
00757 #ifdef __STDC__
00758 long double
00759 __ieee754_lgammal_r (long double x, int *signgamp)
00760 #else
00761 long double
00762 __ieee754_lgammal_r (x, signgamp)
00763      long double x;
00764      int *signgamp;
00765 #endif
00766 {
00767   long double p, q, w, z, nx;
00768   int i, nn;
00769 
00770   *signgamp = 1;
00771 
00772   if (! __finitel (x))
00773     return x * x;
00774 
00775   if (x == 0.0L)
00776     {
00777       if (__signbitl (x))
00778         *signgamp = -1;
00779     }
00780 
00781   if (x < 0.0L)
00782     {
00783       q = -x;
00784       p = __floorl (q);
00785       if (p == q)
00786        return (one / (p - p));
00787       i = p;
00788       if ((i & 1) == 0)
00789        *signgamp = -1;
00790       else
00791        *signgamp = 1;
00792       z = q - p;
00793       if (z > 0.5L)
00794        {
00795          p += 1.0L;
00796          z = p - q;
00797        }
00798       z = q * __sinl (PIL * z);
00799       if (z == 0.0L)
00800        return (*signgamp * huge * huge);
00801       w = __ieee754_lgammal_r (q, &i);
00802       z = __logl (PIL / z) - w;
00803       return (z);
00804     }
00805 
00806   if (x < 13.5L)
00807     {
00808       p = 0.0L;
00809       nx = __floorl (x + 0.5L);
00810       nn = nx;
00811       switch (nn)
00812        {
00813        case 0:
00814          /* log gamma (x + 1) = log(x) + log gamma(x) */
00815          if (x <= 0.125)
00816            {
00817              p = x * neval (x, RN1, NRN1) / deval (x, RD1, NRD1);
00818            }
00819          else if (x <= 0.375)
00820            {
00821              z = x - 0.25L;
00822              p = z * neval (z, RN1r25, NRN1r25) / deval (z, RD1r25, NRD1r25);
00823              p += lgam1r25b;
00824              p += lgam1r25a;
00825            }
00826          else if (x <= 0.625)
00827            {
00828              z = x + (1.0L - x0a);
00829              z = z - x0b;
00830              p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
00831              p = p * z * z;
00832              p = p + y0b;
00833              p = p + y0a;
00834            }
00835          else if (x <= 0.875)
00836            {
00837              z = x - 0.75L;
00838              p = z * neval (z, RN1r75, NRN1r75) / deval (z, RD1r75, NRD1r75);
00839              p += lgam1r75b;
00840              p += lgam1r75a;
00841            }
00842          else
00843            {
00844              z = x - 1.0L;
00845              p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
00846            }
00847          p = p - __logl (x);
00848          break;
00849 
00850        case 1:
00851          if (x < 0.875L)
00852            {
00853              if (x <= 0.625)
00854               {
00855                 z = x + (1.0L - x0a);
00856                 z = z - x0b;
00857                 p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
00858                 p = p * z * z;
00859                 p = p + y0b;
00860                 p = p + y0a;
00861               }
00862              else if (x <= 0.875)
00863               {
00864                 z = x - 0.75L;
00865                 p = z * neval (z, RN1r75, NRN1r75)
00866                       / deval (z, RD1r75, NRD1r75);
00867                 p += lgam1r75b;
00868                 p += lgam1r75a;
00869               }
00870              else
00871               {
00872                 z = x - 1.0L;
00873                 p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
00874               }
00875              p = p - __logl (x);
00876            }
00877          else if (x < 1.0L)
00878            {
00879              z = x - 1.0L;
00880              p = z * neval (z, RNr9, NRNr9) / deval (z, RDr9, NRDr9);
00881            }
00882          else if (x == 1.0L)
00883            p = 0.0L;
00884          else if (x <= 1.125L)
00885            {
00886              z = x - 1.0L;
00887              p = z * neval (z, RN1, NRN1) / deval (z, RD1, NRD1);
00888            }
00889          else if (x <= 1.375)
00890            {
00891              z = x - 1.25L;
00892              p = z * neval (z, RN1r25, NRN1r25) / deval (z, RD1r25, NRD1r25);
00893              p += lgam1r25b;
00894              p += lgam1r25a;
00895            }
00896          else
00897            {
00898              /* 1.375 <= x+x0 <= 1.625 */
00899              z = x - x0a;
00900              z = z - x0b;
00901              p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
00902              p = p * z * z;
00903              p = p + y0b;
00904              p = p + y0a;
00905            }
00906          break;
00907 
00908        case 2:
00909          if (x < 1.625L)
00910            {
00911              z = x - x0a;
00912              z = z - x0b;
00913              p = neval (z, RN1r5, NRN1r5) / deval (z, RD1r5, NRD1r5);
00914              p = p * z * z;
00915              p = p + y0b;
00916              p = p + y0a;
00917            }
00918          else if (x < 1.875L)
00919            {
00920              z = x - 1.75L;
00921              p = z * neval (z, RN1r75, NRN1r75) / deval (z, RD1r75, NRD1r75);
00922              p += lgam1r75b;
00923              p += lgam1r75a;
00924            }
00925          else if (x == 2.0L)
00926            p = 0.0L;
00927          else if (x < 2.375L)
00928            {
00929              z = x - 2.0L;
00930              p = z * neval (z, RN2, NRN2) / deval (z, RD2, NRD2);
00931            }
00932          else
00933            {
00934              z = x - 2.5L;
00935              p = z * neval (z, RN2r5, NRN2r5) / deval (z, RD2r5, NRD2r5);
00936              p += lgam2r5b;
00937              p += lgam2r5a;
00938            }
00939          break;
00940 
00941        case 3:
00942          if (x < 2.75)
00943            {
00944              z = x - 2.5L;
00945              p = z * neval (z, RN2r5, NRN2r5) / deval (z, RD2r5, NRD2r5);
00946              p += lgam2r5b;
00947              p += lgam2r5a;
00948            }
00949          else
00950            {
00951              z = x - 3.0L;
00952              p = z * neval (z, RN3, NRN3) / deval (z, RD3, NRD3);
00953              p += lgam3b;
00954              p += lgam3a;
00955            }
00956          break;
00957 
00958        case 4:
00959          z = x - 4.0L;
00960          p = z * neval (z, RN4, NRN4) / deval (z, RD4, NRD4);
00961          p += lgam4b;
00962          p += lgam4a;
00963          break;
00964 
00965        case 5:
00966          z = x - 5.0L;
00967          p = z * neval (z, RN5, NRN5) / deval (z, RD5, NRD5);
00968          p += lgam5b;
00969          p += lgam5a;
00970          break;
00971 
00972        case 6:
00973          z = x - 6.0L;
00974          p = z * neval (z, RN6, NRN6) / deval (z, RD6, NRD6);
00975          p += lgam6b;
00976          p += lgam6a;
00977          break;
00978 
00979        case 7:
00980          z = x - 7.0L;
00981          p = z * neval (z, RN7, NRN7) / deval (z, RD7, NRD7);
00982          p += lgam7b;
00983          p += lgam7a;
00984          break;
00985 
00986        case 8:
00987          z = x - 8.0L;
00988          p = z * neval (z, RN8, NRN8) / deval (z, RD8, NRD8);
00989          p += lgam8b;
00990          p += lgam8a;
00991          break;
00992 
00993        case 9:
00994          z = x - 9.0L;
00995          p = z * neval (z, RN9, NRN9) / deval (z, RD9, NRD9);
00996          p += lgam9b;
00997          p += lgam9a;
00998          break;
00999 
01000        case 10:
01001          z = x - 10.0L;
01002          p = z * neval (z, RN10, NRN10) / deval (z, RD10, NRD10);
01003          p += lgam10b;
01004          p += lgam10a;
01005          break;
01006 
01007        case 11:
01008          z = x - 11.0L;
01009          p = z * neval (z, RN11, NRN11) / deval (z, RD11, NRD11);
01010          p += lgam11b;
01011          p += lgam11a;
01012          break;
01013 
01014        case 12:
01015          z = x - 12.0L;
01016          p = z * neval (z, RN12, NRN12) / deval (z, RD12, NRD12);
01017          p += lgam12b;
01018          p += lgam12a;
01019          break;
01020 
01021        case 13:
01022          z = x - 13.0L;
01023          p = z * neval (z, RN13, NRN13) / deval (z, RD13, NRD13);
01024          p += lgam13b;
01025          p += lgam13a;
01026          break;
01027        }
01028       return p;
01029     }
01030 
01031   if (x > MAXLGM)
01032     return (*signgamp * huge * huge);
01033 
01034   q = ls2pi - x;
01035   q = (x - 0.5L) * __logl (x) + q;
01036   if (x > 1.0e18L)
01037     return (q);
01038 
01039   p = 1.0L / (x * x);
01040   q += neval (p, RASY, NRASY) / x;
01041   return (q);
01042 }