370 lines
8.2 KiB
C
370 lines
8.2 KiB
C
/* lgam()
|
|
*
|
|
* Natural logarithm of gamma function
|
|
*
|
|
*
|
|
*
|
|
* SYNOPSIS:
|
|
*
|
|
* double x, y, __lgamma_r();
|
|
* int* sgngam;
|
|
* y = __lgamma_r( x, sgngam );
|
|
*
|
|
* double x, y, lgamma();
|
|
* y = lgamma( x);
|
|
*
|
|
*
|
|
*
|
|
* DESCRIPTION:
|
|
*
|
|
* Returns the base e (2.718...) logarithm of the absolute
|
|
* value of the gamma function of the argument. In the reentrant
|
|
* version, the sign (+1 or -1) of the gamma function is returned
|
|
* in the variable referenced by sgngam.
|
|
*
|
|
* For arguments greater than 13, the logarithm of the gamma
|
|
* function is approximated by the logarithmic version of
|
|
* Stirling's formula using a polynomial approximation of
|
|
* degree 4. Arguments between -33 and +33 are reduced by
|
|
* recurrence to the interval [2,3] of a rational approximation.
|
|
* The cosecant reflection formula is employed for arguments
|
|
* less than -33.
|
|
*
|
|
* Arguments greater than MAXLGM return MAXNUM and an error
|
|
* message. MAXLGM = 2.035093e36 for DEC
|
|
* arithmetic or 2.556348e305 for IEEE arithmetic.
|
|
*
|
|
*
|
|
*
|
|
* ACCURACY:
|
|
*
|
|
*
|
|
* arithmetic domain # trials peak rms
|
|
* DEC 0, 3 7000 5.2e-17 1.3e-17
|
|
* DEC 2.718, 2.035e36 5000 3.9e-17 9.9e-18
|
|
* IEEE 0, 3 28000 5.4e-16 1.1e-16
|
|
* IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17
|
|
* The error criterion was relative when the function magnitude
|
|
* was greater than one but absolute when it was less than one.
|
|
*
|
|
* The following test used the relative error criterion, though
|
|
* at certain points the relative error could be much higher than
|
|
* indicated.
|
|
* IEEE -200, -4 10000 4.8e-16 1.3e-16
|
|
*
|
|
*/
|
|
|
|
/*
|
|
* Cephes Math Library Release 2.8: June, 2000
|
|
* Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
|
|
*/
|
|
|
|
/*
|
|
* 26-11-2002 Modified for mingw.
|
|
* Danny Smith <dannysmith@users.sourceforge.net>
|
|
*/
|
|
|
|
|
|
#ifndef __MINGW32__
|
|
#include "mconf.h"
|
|
#ifdef ANSIPROT
|
|
extern double pow ( double, double );
|
|
extern double log ( double );
|
|
extern double exp ( double );
|
|
extern double sin ( double );
|
|
extern double polevl ( double, void *, int );
|
|
extern double p1evl ( double, void *, int );
|
|
extern double floor ( double );
|
|
extern double fabs ( double );
|
|
extern int isnan ( double );
|
|
extern int isfinite ( double );
|
|
#else
|
|
double pow(), log(), exp(), sin(), polevl(), p1evl(), floor(), fabs();
|
|
int isnan(), isfinite();
|
|
#endif
|
|
#ifdef INFINITIES
|
|
extern double INFINITY;
|
|
#endif
|
|
#ifdef NANS
|
|
extern double NAN;
|
|
#endif
|
|
#else /* __MINGW32__ */
|
|
#include "cephes_mconf.h"
|
|
#endif /* __MINGW32__ */
|
|
|
|
|
|
/* A[]: Stirling's formula expansion of log gamma
|
|
* B[], C[]: log gamma function between 2 and 3
|
|
*/
|
|
#ifdef UNK
|
|
static double A[] = {
|
|
8.11614167470508450300E-4,
|
|
-5.95061904284301438324E-4,
|
|
7.93650340457716943945E-4,
|
|
-2.77777777730099687205E-3,
|
|
8.33333333333331927722E-2
|
|
};
|
|
static double B[] = {
|
|
-1.37825152569120859100E3,
|
|
-3.88016315134637840924E4,
|
|
-3.31612992738871184744E5,
|
|
-1.16237097492762307383E6,
|
|
-1.72173700820839662146E6,
|
|
-8.53555664245765465627E5
|
|
};
|
|
static double C[] = {
|
|
/* 1.00000000000000000000E0, */
|
|
-3.51815701436523470549E2,
|
|
-1.70642106651881159223E4,
|
|
-2.20528590553854454839E5,
|
|
-1.13933444367982507207E6,
|
|
-2.53252307177582951285E6,
|
|
-2.01889141433532773231E6
|
|
};
|
|
/* log( sqrt( 2*pi ) ) */
|
|
static double LS2PI = 0.91893853320467274178;
|
|
#define MAXLGM 2.556348e305
|
|
static double LOGPI = 1.14472988584940017414;
|
|
#endif
|
|
|
|
#ifdef DEC
|
|
static const unsigned short A[] = {
|
|
0035524,0141201,0034633,0031405,
|
|
0135433,0176755,0126007,0045030,
|
|
0035520,0006371,0003342,0172730,
|
|
0136066,0005540,0132605,0026407,
|
|
0037252,0125252,0125252,0125132
|
|
};
|
|
static const unsigned short B[] = {
|
|
0142654,0044014,0077633,0035410,
|
|
0144027,0110641,0125335,0144760,
|
|
0144641,0165637,0142204,0047447,
|
|
0145215,0162027,0146246,0155211,
|
|
0145322,0026110,0010317,0110130,
|
|
0145120,0061472,0120300,0025363
|
|
};
|
|
static const unsigned short C[] = {
|
|
/*0040200,0000000,0000000,0000000*/
|
|
0142257,0164150,0163630,0112622,
|
|
0143605,0050153,0156116,0135272,
|
|
0144527,0056045,0145642,0062332,
|
|
0145213,0012063,0106250,0001025,
|
|
0145432,0111254,0044577,0115142,
|
|
0145366,0071133,0050217,0005122
|
|
};
|
|
/* log( sqrt( 2*pi ) ) */
|
|
static const unsigned short LS2P[] = {040153,037616,041445,0172645,};
|
|
#define LS2PI *(double *)LS2P
|
|
#define MAXLGM 2.035093e36
|
|
static const unsigned short LPI[4] = {
|
|
0040222,0103202,0043475,0006750,
|
|
};
|
|
#define LOGPI *(double *)LPI
|
|
|
|
#endif
|
|
|
|
#ifdef IBMPC
|
|
static const unsigned short A[] = {
|
|
0x6661,0x2733,0x9850,0x3f4a,
|
|
0xe943,0xb580,0x7fbd,0xbf43,
|
|
0x5ebb,0x20dc,0x019f,0x3f4a,
|
|
0xa5a1,0x16b0,0xc16c,0xbf66,
|
|
0x554b,0x5555,0x5555,0x3fb5
|
|
};
|
|
static const unsigned short B[] = {
|
|
0x6761,0x8ff3,0x8901,0xc095,
|
|
0xb93e,0x355b,0xf234,0xc0e2,
|
|
0x89e5,0xf890,0x3d73,0xc114,
|
|
0xdb51,0xf994,0xbc82,0xc131,
|
|
0xf20b,0x0219,0x4589,0xc13a,
|
|
0x055e,0x5418,0x0c67,0xc12a
|
|
};
|
|
static const unsigned short C[] = {
|
|
/*0x0000,0x0000,0x0000,0x3ff0,*/
|
|
0x12b2,0x1cf3,0xfd0d,0xc075,
|
|
0xd757,0x7b89,0xaa0d,0xc0d0,
|
|
0x4c9b,0xb974,0xeb84,0xc10a,
|
|
0x0043,0x7195,0x6286,0xc131,
|
|
0xf34c,0x892f,0x5255,0xc143,
|
|
0xe14a,0x6a11,0xce4b,0xc13e
|
|
};
|
|
/* log( sqrt( 2*pi ) ) */
|
|
static const union
|
|
{
|
|
unsigned short s[4];
|
|
double d;
|
|
} ls2p = {{0xbeb5,0xc864,0x67f1,0x3fed}};
|
|
#define LS2PI (ls2p.d)
|
|
#define MAXLGM 2.556348e305
|
|
/* log (pi) */
|
|
static const union
|
|
{
|
|
unsigned short s[4];
|
|
double d;
|
|
} lpi = {{0xa1bd,0x48e7,0x50d0,0x3ff2}};
|
|
#define LOGPI (lpi.d)
|
|
#endif
|
|
|
|
#ifdef MIEEE
|
|
static const unsigned short A[] = {
|
|
0x3f4a,0x9850,0x2733,0x6661,
|
|
0xbf43,0x7fbd,0xb580,0xe943,
|
|
0x3f4a,0x019f,0x20dc,0x5ebb,
|
|
0xbf66,0xc16c,0x16b0,0xa5a1,
|
|
0x3fb5,0x5555,0x5555,0x554b
|
|
};
|
|
static const unsigned short B[] = {
|
|
0xc095,0x8901,0x8ff3,0x6761,
|
|
0xc0e2,0xf234,0x355b,0xb93e,
|
|
0xc114,0x3d73,0xf890,0x89e5,
|
|
0xc131,0xbc82,0xf994,0xdb51,
|
|
0xc13a,0x4589,0x0219,0xf20b,
|
|
0xc12a,0x0c67,0x5418,0x055e
|
|
};
|
|
static const unsigned short C[] = {
|
|
0xc075,0xfd0d,0x1cf3,0x12b2,
|
|
0xc0d0,0xaa0d,0x7b89,0xd757,
|
|
0xc10a,0xeb84,0xb974,0x4c9b,
|
|
0xc131,0x6286,0x7195,0x0043,
|
|
0xc143,0x5255,0x892f,0xf34c,
|
|
0xc13e,0xce4b,0x6a11,0xe14a
|
|
};
|
|
/* log( sqrt( 2*pi ) ) */
|
|
static const union
|
|
{
|
|
unsigned short s[4];
|
|
double d;
|
|
} ls2p = {{0x3fed,0x67f1,0xc864,0xbeb5}};
|
|
#define LS2PI ls2p.d
|
|
#define MAXLGM 2.556348e305
|
|
/* log (pi) */
|
|
static const union
|
|
{
|
|
unsigned short s[4];
|
|
double d;
|
|
} lpi = {{0x3ff2, 0x50d0, 0x48e7, 0xa1bd}};
|
|
#define LOGPI (lpi.d)
|
|
#endif
|
|
|
|
|
|
/* Logarithm of gamma function */
|
|
/* Reentrant version */
|
|
|
|
double __lgamma_r(double x, int* sgngam)
|
|
{
|
|
double p, q, u, w, z;
|
|
int i;
|
|
|
|
*sgngam = 1;
|
|
#ifdef NANS
|
|
if( isnan(x) )
|
|
return(x);
|
|
#endif
|
|
|
|
#ifdef INFINITIES
|
|
if( !isfinite(x) )
|
|
return(INFINITY);
|
|
#endif
|
|
|
|
if( x < -34.0 )
|
|
{
|
|
q = -x;
|
|
w = __lgamma_r(q, sgngam); /* note this modifies sgngam! */
|
|
p = floor(q);
|
|
if( p == q )
|
|
{
|
|
lgsing:
|
|
_SET_ERRNO(EDOM);
|
|
mtherr( "lgam", SING );
|
|
#ifdef INFINITIES
|
|
return (INFINITY);
|
|
#else
|
|
return (MAXNUM);
|
|
#endif
|
|
}
|
|
i = p;
|
|
if( (i & 1) == 0 )
|
|
*sgngam = -1;
|
|
else
|
|
*sgngam = 1;
|
|
z = q - p;
|
|
if( z > 0.5 )
|
|
{
|
|
p += 1.0;
|
|
z = p - q;
|
|
}
|
|
z = q * sin( PI * z );
|
|
if( z == 0.0 )
|
|
goto lgsing;
|
|
/* z = log(PI) - log( z ) - w;*/
|
|
z = LOGPI - log( z ) - w;
|
|
return( z );
|
|
}
|
|
|
|
if( x < 13.0 )
|
|
{
|
|
z = 1.0;
|
|
p = 0.0;
|
|
u = x;
|
|
while( u >= 3.0 )
|
|
{
|
|
p -= 1.0;
|
|
u = x + p;
|
|
z *= u;
|
|
}
|
|
while( u < 2.0 )
|
|
{
|
|
if( u == 0.0 )
|
|
goto lgsing;
|
|
z /= u;
|
|
p += 1.0;
|
|
u = x + p;
|
|
}
|
|
if( z < 0.0 )
|
|
{
|
|
*sgngam = -1;
|
|
z = -z;
|
|
}
|
|
else
|
|
*sgngam = 1;
|
|
if( u == 2.0 )
|
|
return( log(z) );
|
|
p -= 2.0;
|
|
x = x + p;
|
|
p = x * polevl( x, B, 5 ) / p1evl( x, C, 6);
|
|
return( log(z) + p );
|
|
}
|
|
|
|
if( x > MAXLGM )
|
|
{
|
|
_SET_ERRNO(ERANGE);
|
|
mtherr( "lgamma", OVERFLOW );
|
|
#ifdef INFINITIES
|
|
return( *sgngam * INFINITY );
|
|
#else
|
|
return( *sgngam * MAXNUM );
|
|
#endif
|
|
}
|
|
|
|
q = ( x - 0.5 ) * log(x) - x + LS2PI;
|
|
if( x > 1.0e8 )
|
|
return( q );
|
|
|
|
p = 1.0/(x*x);
|
|
if( x >= 1000.0 )
|
|
q += (( 7.9365079365079365079365e-4 * p
|
|
- 2.7777777777777777777778e-3) *p
|
|
+ 0.0833333333333333333333) / x;
|
|
else
|
|
q += polevl( p, A, 4 ) / x;
|
|
return( q );
|
|
}
|
|
|
|
/* This is the C99 version */
|
|
|
|
double lgamma(double x)
|
|
{
|
|
int local_sgngam=0;
|
|
return (__lgamma_r(x, &local_sgngam));
|
|
}
|