1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
|
#include "dietlibm.h"
/*--------------------------------------------------------------------------*
Name gamma, lgamma - gamma function
Usage double gamma (double x);
double lgamma(double x);
extern int signgam;
Prototype in math.h
Description gamma returns the logarithm of the absolute value of the
gamma function. So it is possible â(x) for very large x.
The sign is stored in signgam, a extern variable
overwritten during every call to gamma(). lgamma() is
a synonym for gamma().
You can calculate â(x) by the following sequence:
double gammafunction(double x)
{ double y=exp(gamma(x));
return signgam ? -y : +y;
}
Return value gamma returns a value in range (-0.1208, +oo). For a input
value of zero, it returns +oo and errno is set to:
ERANGE Result out of range
*---------------------------------------------------------------------------*/
#include <stdlib.h>
#include <math.h>
#define B0 + 1.0l/ 6/ 1/ 2
#define B1 - 1.0l/ 30/ 3/ 4
#define B2 + 1.0l/ 42/ 5/ 6
#define B3 - 1.0l/ 30/ 7/ 8
#define B4 + 5.0l/ 66/ 9/10
#define B5 - 691.0l/2730/11/12
#define B6 + 7.0l/ 6/13/14
#define B7 - 3617.0l/ 510/15/16
#define B8 + 43867.0l/ 798/17/18
#define B9 - 174611.0l/ 330/19/20
#define B10 + 854513.0l/ 138/21/22
#define B11 - 236364091.0l/2730/23/24
#define B12 + 8553103.0l/ 6/25/26
static const double coeff[] = { B0, B1, B2, B3, B4, B5, B6, B7, B8, B9, B10 };
int signgam;
#define EXPL(x) (((short *)&x)[4] & 0x7FFF)
static double logfact ( long double x )
{
long double z = 2. * M_PI * x;
register int e = EXPL (x);
static unsigned char list [] = { 6, 4, 3, 3, 2, 2 };
return (log(x) - 1) * x + 0.5*log(z) + __poly (1./(x*x), e<0x4003 ? 10 : (e>0x4008 ? 1 : list [e-0x4003] ), coeff) / x;
}
double lgamma ( double x )
{
register int k = floor (x);
long double w;
long double y;
long double z;
signgam = 0;
if ( k >= 7 )
return logfact (x-1);
if ( k == x )
switch (k) {
case 1 :
case 2 : return 0.000000000000000000000000000l;
case 3 : return 0.693147180559945309432805516l;
case 4 : return 1.791759469228055000858148560l;
case 5 : return 3.178053830347945619723759592l;
case 6 : return 4.787491742782045994244981560l;
default: return 1./0.; /* ignore the gcc warning, this is intentional */
}
z = logfact (y = x - k + 7.0 - 1);
w = 1;
for ( k = 7 - k; k--; )
w *= y, y -= 1.;
signgam = k >= 0 ? 0 : k & 1;
return z - log (w);
}
double gamma ( double val ) __attribute__ ((weak,alias("lgamma")));
|
