e_lgammal_r.c

Glibc 2.3.2源代码(解压后有100多M)

/*                                                      lgammal
 *
 *      Natural logarithm of gamma function
 *
 *
 *
 * SYNOPSIS:
 *
 * long double x, y, lgammal();
 * extern int sgngam;
 *
 * y = lgammal(x);
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns the base e (2.718...) logarithm of the absolute
 * value of the gamma function of the argument.
 * The sign (+1 or -1) of the gamma function is returned in a
 * global (extern) variable named sgngam.
 *
 * The positive domain is partitioned into numerous segments for approximation.
 * For x > 10,
 *   log gamma(x) = (x - 0.5) log(x) - x + log sqrt(2 pi) + 1/x R(1/x^2)
 * Near the minimum at x = x0 = 1.46... the approximation is
 *   log gamma(x0 + z) = log gamma(x0) + z^2 P(z)/Q(z)
 * for small z.
 * Elsewhere between 0 and 10,
 *   log gamma(n + z) = log gamma(n) + z P(z)/Q(z)
 * for various selected n and small z.
 *
 * The cosecant reflection formula is employed for negative arguments.
 *
 *
 *
 * ACCURACY:
 *
 *
 * arithmetic      domain        # trials     peak         rms
 *                                            Relative error:
 *    IEEE         10, 30         100000     3.9e-34     9.8e-35
 *    IEEE          0, 10         100000     3.8e-34     5.3e-35
 *                                            Absolute error:
 *    IEEE         -10, 0         100000     8.0e-34     8.0e-35
 *    IEEE         -30, -10       100000     4.4e-34     1.0e-34
 *    IEEE        -100, 100       100000                 1.0e-34
 *
 * The absolute error criterion is the same as relative error
 * when the function magnitude is greater than one but it is absolute
 * when the magnitude is less than one.
 *
 */

/* Copyright 2001 by Stephen L. Moshier <[e-mail]>

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this library; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */

#include "math.h"
#include "math_private.h"

static const long double PIL = 3.1415926535897932384626433832795028841972E0L;
static const long double MAXLGM = 1.0485738685148938358098967157129705071571E4928L;
static const long double one = 1.0L;
static const long double zero = 0.0L;
static const long double huge = 1.0e4000L;

/* log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x P(1/x^2)
   1/x <= 0.0741 (x >= 13.495...)
   Peak relative error 1.5e-36  */
static const long double ls2pi = 9.1893853320467274178032973640561763986140E-1L;
#define NRASY 12
static const long double RASY[NRASY + 1] =
{
  8.333333333333333333333333333310437112111E-2L,
 -2.777777777777777777777774789556228296902E-3L,
  7.936507936507936507795933938448586499183E-4L,
 -5.952380952380952041799269756378148574045E-4L,
  8.417508417507928904209891117498524452523E-4L,
 -1.917526917481263997778542329739806086290E-3L,
  6.410256381217852504446848671499409919280E-3L,
 -2.955064066900961649768101034477363301626E-2L,
  1.796402955865634243663453415388336954675E-1L,
 -1.391522089007758553455753477688592767741E0L,
  1.326130089598399157988112385013829305510E1L,
 -1.420412699593782497803472576479997819149E2L,
  1.218058922427762808938869872528846787020E3L
};


/* log gamma(x+13) = log gamma(13) +  x P(x)/Q(x)
   -0.5 <= x <= 0.5
   12.5 <= x+13 <= 13.5
   Peak relative error 1.1e-36  */
static const long double lgam13a = 1.9987213134765625E1L;
static const long double lgam13b = 1.3608962611495173623870550785125024484248E-6L;
#define NRN13 7
static const long double RN13[NRN13 + 1] =
{
  8.591478354823578150238226576156275285700E11L,
  2.347931159756482741018258864137297157668E11L,
  2.555408396679352028680662433943000804616E10L,
  1.408581709264464345480765758902967123937E9L,
  4.126759849752613822953004114044451046321E7L,
  6.133298899622688505854211579222889943778E5L,
  3.929248056293651597987893340755876578072E3L,
  6.850783280018706668924952057996075215223E0L
};
#define NRD13 6
static const long double RD13[NRD13 + 1] =
{
  3.401225382297342302296607039352935541669E11L,
  8.756765276918037910363513243563234551784E10L,
  8.873913342866613213078554180987647243903E9L,
  4.483797255342763263361893016049310017973E8L,
  1.178186288833066430952276702931512870676E7L,
  1.519928623743264797939103740132278337476E5L,
  7.989298844938119228411117593338850892311E2L
 /* 1.0E0L */
};


/* log gamma(x+12) = log gamma(12) +  x P(x)/Q(x)
   -0.5 <= x <= 0.5
   11.5 <= x+12 <= 12.5
   Peak relative error 4.1e-36  */
static const long double lgam12a = 1.75023040771484375E1L;
static const long double lgam12b = 3.7687254483392876529072161996717039575982E-6L;
#define NRN12 7
static const long double RN12[NRN12 + 1] =
{
  4.709859662695606986110997348630997559137E11L,
  1.398713878079497115037857470168777995230E11L,
  1.654654931821564315970930093932954900867E10L,
  9.916279414876676861193649489207282144036E8L,
  3.159604070526036074112008954113411389879E7L,
  5.109099197547205212294747623977502492861E5L,
  3.563054878276102790183396740969279826988E3L,
  6.769610657004672719224614163196946862747E0L
};
#define NRD12 6
static const long double RD12[NRD12 + 1] =
{
  1.928167007860968063912467318985802726613E11L,
  5.383198282277806237247492369072266389233E10L,
  5.915693215338294477444809323037871058363E9L,
  3.241438287570196713148310560147925781342E8L,
  9.236680081763754597872713592701048455890E6L,
  1.292246897881650919242713651166596478850E5L,
  7.366532445427159272584194816076600211171E2L
 /* 1.0E0L */
};


/* log gamma(x+11) = log gamma(11) +  x P(x)/Q(x)
   -0.5 <= x <= 0.5
   10.5 <= x+11 <= 11.5
   Peak relative error 1.8e-35  */
static const long double lgam11a = 1.5104400634765625E1L;
static const long double lgam11b = 1.1938309890295225709329251070371882250744E-5L;
#define NRN11 7
static const long double RN11[NRN11 + 1] =
{
  2.446960438029415837384622675816736622795E11L,
  7.955444974446413315803799763901729640350E10L,
  1.030555327949159293591618473447420338444E10L,
  6.765022131195302709153994345470493334946E8L,
  2.361892792609204855279723576041468347494E7L,
  4.186623629779479136428005806072176490125E5L,
  3.202506022088912768601325534149383594049E3L,
  6.681356101133728289358838690666225691363E0L
};
#define NRD11 6
static const long double RD11[NRD11 + 1] =
{
  1.040483786179428590683912396379079477432E11L,
  3.172251138489229497223696648369823779729E10L,
  3.806961885984850433709295832245848084614E9L,
  2.278070344022934913730015420611609620171E8L,
  7.089478198662651683977290023829391596481E6L,
  1.083246385105903533237139380509590158658E5L,
  6.744420991491385145885727942219463243597E2L
 /* 1.0E0L */
};


/* log gamma(x+10) = log gamma(10) +  x P(x)/Q(x)
   -0.5 <= x <= 0.5
   9.5 <= x+10 <= 10.5
   Peak relative error 5.4e-37  */
static const long double lgam10a = 1.280181884765625E1L;
static const long double lgam10b = 8.6324252196112077178745667061642811492557E-6L;
#define NRN10 7
static const long double RN10[NRN10 + 1] =
{
  -1.239059737177249934158597996648808363783E14L,
  -4.725899566371458992365624673357356908719E13L,
  -7.283906268647083312042059082837754850808E12L,
  -5.802855515464011422171165179767478794637E11L,
  -2.532349691157548788382820303182745897298E10L,
  -5.884260178023777312587193693477072061820E8L,
  -6.437774864512125749845840472131829114906E6L,
  -2.350975266781548931856017239843273049384E4L
};
#define NRD10 7
static const long double RD10[NRD10 + 1] =
{
  -5.502645997581822567468347817182347679552E13L,
  -1.970266640239849804162284805400136473801E13L,
  -2.819677689615038489384974042561531409392E12L,
  -2.056105863694742752589691183194061265094E11L,
  -8.053670086493258693186307810815819662078E9L,
  -1.632090155573373286153427982504851867131E8L,
  -1.483575879240631280658077826889223634921E6L,
  -4.002806669713232271615885826373550502510E3L
 /* 1.0E0L */
};


/* log gamma(x+9) = log gamma(9) +  x P(x)/Q(x)
   -0.5 <= x <= 0.5
   8.5 <= x+9 <= 9.5
   Peak relative error 3.6e-36  */
static const long double lgam9a = 1.06045989990234375E1L;
static const long double lgam9b = 3.9037218127284172274007216547549861681400E-6L;
#define NRN9 7
static const long double RN9[NRN9 + 1] =
{
  -4.936332264202687973364500998984608306189E13L,
  -2.101372682623700967335206138517766274855E13L,
  -3.615893404644823888655732817505129444195E12L,
  -3.217104993800878891194322691860075472926E11L,
  -1.568465330337375725685439173603032921399E10L,
  -4.073317518162025744377629219101510217761E8L,
  -4.983232096406156139324846656819246974500E6L,
  -2.036280038903695980912289722995505277253E4L
};
#define NRD9 7
static const long double RD9[NRD9 + 1] =
{
  -2.306006080437656357167128541231915480393E13L,
  -9.183606842453274924895648863832233799950E12L,
  -1.461857965935942962087907301194381010380E12L,
  -1.185728254682789754150068652663124298303E11L,
  -5.166285094703468567389566085480783070037E9L,
  -1.164573656694603024184768200787835094317E8L,
  -1.177343939483908678474886454113163527909E6L,
  -3.529391059783109732159524500029157638736E3L
  /* 1.0E0L */
};


/* log gamma(x+8) = log gamma(8) +  x P(x)/Q(x)
   -0.5 <= x <= 0.5
   7.5 <= x+8 <= 8.5
   Peak relative error 2.4e-37  */
static const long double lgam8a = 8.525146484375E0L;
static const long double lgam8b = 1.4876690414300165531036347125050759667737E-5L;
#define NRN8 8
static const long double RN8[NRN8 + 1] =
{
  6.600775438203423546565361176829139703289E11L,
  3.406361267593790705240802723914281025800E11L,
  7.222460928505293914746983300555538432830E10L,
  8.102984106025088123058747466840656458342E9L,
  5.157620015986282905232150979772409345927E8L,
  1.851445288272645829028129389609068641517E7L,
  3.489261702223124354745894067468953756656E5L,
  2.892095396706665774434217489775617756014E3L,
  6.596977510622195827183948478627058738034E0L
};
#define NRD8 7
static const long double RD8[NRD8 + 1] =
{
  3.274776546520735414638114828622673016920E11L,
  1.581811207929065544043963828487733970107E11L,
  3.108725655667825188135393076860104546416E10L,
  3.193055010502912617128480163681842165730E9L,
  1.830871482669835106357529710116211541839E8L,
  5.790862854275238129848491555068073485086E6L,
  9.305213264307921522842678835618803553589E4L,
  6.216974105861848386918949336819572333622E2L
  /* 1.0E0L */
};


/* log gamma(x+7) = log gamma(7) +  x P(x)/Q(x)
   -0.5 <= x <= 0.5
   6.5 <= x+7 <= 7.5
   Peak relative error 3.2e-36  */
static const long double lgam7a = 6.5792388916015625E0L;
static const long double lgam7b = 1.2320408538495060178292903945321122583007E-5L;
#define NRN7 8
static const long double RN7[NRN7 + 1] =
{
  2.065019306969459407636744543358209942213E11L,
  1.226919919023736909889724951708796532847E11L,
  2.996157990374348596472241776917953749106E10L,
  3.873001919306801037344727168434909521030E9L,
  2.841575255593761593270885753992732145094E8L,
  1.176342515359431913664715324652399565551E7L,
  2.558097039684188723597519300356028511547E5L,
  2.448525238332609439023786244782810774702E3L,
  6.460280377802030953041566617300902020435E0L
};
#define NRD7 7
static const long double RD7[NRD7 + 1] =
{
  1.102646614598516998880874785339049304483E11L,
  6.099297512712715445879759589407189290040E10L,
  1.372898136289611312713283201112060238351E10L,
  1.615306270420293159907951633566635172343E9L,
  1.061114435798489135996614242842561967459E8L,
  3.845638971184305248268608902030718674691E6L,
  7.081730675423444975703917836972720495507E4L,
  5.423122582741398226693137276201344096370E2L
  /* 1.0E0L */
};


/* log gamma(x+6) = log gamma(6) +  x P(x)/Q(x)
   -0.5 <= x <= 0.5
   5.5 <= x+6 <= 6.5
   Peak relative error 6.2e-37  */
static const long double lgam6a = 4.7874908447265625E0L;
static const long double lgam6b = 8.9805548349424770093452324304839959231517E-7L;
#define NRN6 8
static const long double RN6[NRN6 + 1] =
{
  -3.538412754670746879119162116819571823643E13L,
  -2.613432593406849155765698121483394257148E13L,
  -8.020670732770461579558867891923784753062E12L,
  -1.322227822931250045347591780332435433420E12L,
  -1.262809382777272476572558806855377129513E11L,