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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
|
/*
* Copyright (c) 1985 Regents of the University of California.
* All rights reserved. The Berkeley software License Agreement
* specifies the terms and conditions for redistribution.
*/
#if 0
#ifndef lint
static char sccsid[] = "@(#)lgamma.c 5.3 (Berkeley) 9/22/88";
#endif /* not lint */
#endif
/*
C program for floating point log Gamma function
lgamma(x) computes the log of the absolute
value of the Gamma function.
The sign of the Gamma function is returned in the
external quantity signgam.
The coefficients for expansion around zero
are #5243 from Hart & Cheney; for expansion
around infinity they are #5404.
Calls log, floor and sin.
*/
#if 0
#include "mathimpl.h"
#else
#include <math.h>
#endif
#if defined(vax)||defined(tahoe)
#include <errno.h>
#endif /* defined(vax)||defined(tahoe) */
int signgam;
static const double goobie = 0.9189385332046727417803297; /* log(2*pi)/2 */
static const double pi = 3.1415926535897932384626434;
#define M 6
#define N 8
static const double p1[] = {
0.83333333333333101837e-1,
-.277777777735865004e-2,
0.793650576493454e-3,
-.5951896861197e-3,
0.83645878922e-3,
-.1633436431e-2,
};
static const double p2[] = {
-.42353689509744089647e5,
-.20886861789269887364e5,
-.87627102978521489560e4,
-.20085274013072791214e4,
-.43933044406002567613e3,
-.50108693752970953015e2,
-.67449507245925289918e1,
0.0,
};
static const double q2[] = {
-.42353689509744090010e5,
-.29803853309256649932e4,
0.99403074150827709015e4,
-.15286072737795220248e4,
-.49902852662143904834e3,
0.18949823415702801641e3,
-.23081551524580124562e2,
0.10000000000000000000e1,
};
static double pos(double), neg(double), asym(double);
double
lgamma(double arg)
{
signgam = 1;
if(arg <= 0.) return(neg(arg));
if(arg > 8.) return(asym(arg));
return(log(pos(arg)));
}
static double
asym(double arg)
{
double n, argsq;
int i;
/* Avoid overflow of arg*arg. */
n = 0.0;
if(arg > 1.0e8)
goto noasym;
argsq = 1./(arg*arg);
for(i=M-1; i>=0; i--){
n = n*argsq + p1[i];
}
noasym:
return((arg-.5)*log(arg) - arg + goobie + n/arg);
}
static double
neg(double arg)
{
double t;
arg = -arg;
/*
* to see if arg were a true integer, the old code used the
* mathematically correct observation:
* sin(n*pi) = 0 <=> n is an integer.
* but in finite precision arithmetic, sin(n*PI) will NEVER
* be zero simply because n*PI is a rational number. hence
* it failed to work with our newer, more accurate sin()
* which uses true pi to do the argument reduction...
* temp = sin(pi*arg);
*/
t = floor(arg);
if (arg - t > 0.5e0)
t += 1.e0; /* t := integer nearest arg */
#if defined(vax)||defined(tahoe)
if (arg == t) {
return(infnan(ERANGE)); /* +INF */
}
#endif /* defined(vax)||defined(tahoe) */
signgam = (int) (t - 2*floor(t/2)); /* signgam = 1 if t was odd, */
/* 0 if t was even */
signgam = signgam - 1 + signgam; /* signgam = 1 if t was odd, */
/* -1 if t was even */
t = arg - t; /* -0.5 <= t <= 0.5 */
if (t < 0.e0) {
t = -t;
signgam = -signgam;
}
return(-log(arg*pos(arg)*sin(pi*t)/pi));
}
static double
pos(double arg)
{
double n, d, s;
register i;
if(arg < 2.) return(pos(arg+1.)/arg);
if(arg > 3.) return((arg-1.)*pos(arg-1.));
s = arg - 2.;
for(n=0,d=0,i=N-1; i>=0; i--){
n = n*s + p2[i];
d = d*s + q2[i];
}
return(n/d);
}
|
