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
|
/* erf.c
reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten
(New Algorithm handbook in C language) (Gijyutsu hyouron
sha, Tokyo, 1991) p.227 [in Japanese] */
#include <stdio.h>
#include <math.h>
#ifdef _WIN32
# include <float.h>
# if !defined __MINGW32__ || defined __NO_ISOCEXT
# ifndef isnan
# define isnan(x) _isnan(x)
# endif
# ifndef isinf
# define isinf(x) (!_finite(x) && !_isnan(x))
# endif
# ifndef finite
# define finite(x) _finite(x)
# endif
# endif
#endif
static double q_gamma(double, double, double);
/* Incomplete gamma function
1 / Gamma(a) * Int_0^x exp(-t) t^(a-1) dt */
static double p_gamma(a, x, loggamma_a)
double a, x, loggamma_a;
{
int k;
double result, term, previous;
if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a);
if (x == 0) return 0;
result = term = exp(a * log(x) - x - loggamma_a) / a;
for (k = 1; k < 1000; k++) {
term *= x / (a + k);
previous = result; result += term;
if (result == previous) return result;
}
fprintf(stderr, "erf.c:%d:p_gamma() could not converge.", __LINE__);
return result;
}
/* Incomplete gamma function
1 / Gamma(a) * Int_x^inf exp(-t) t^(a-1) dt */
static double q_gamma(a, x, loggamma_a)
double a, x, loggamma_a;
{
int k;
double result, w, temp, previous;
double la = 1, lb = 1 + x - a; /* Laguerre polynomial */
if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a);
w = exp(a * log(x) - x - loggamma_a);
result = w / lb;
for (k = 2; k < 1000; k++) {
temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;
la = lb; lb = temp;
w *= (k - 1 - a) / k;
temp = w / (la * lb);
previous = result; result += temp;
if (result == previous) return result;
}
fprintf(stderr, "erf.c:%d:q_gamma() could not converge.", __LINE__);
return result;
}
#define LOG_PI_OVER_2 0.572364942924700087071713675675 /* log_e(PI)/2 */
double erf(x)
double x;
{
if (!finite(x)) {
if (isnan(x)) return x; /* erf(NaN) = NaN */
return (x>0 ? 1.0 : -1.0); /* erf(+-inf) = +-1.0 */
}
if (x >= 0) return p_gamma(0.5, x * x, LOG_PI_OVER_2);
else return - p_gamma(0.5, x * x, LOG_PI_OVER_2);
}
double erfc(x)
double x;
{
if (!finite(x)) {
if (isnan(x)) return x; /* erfc(NaN) = NaN */
return (x>0 ? 0.0 : 2.0); /* erfc(+-inf) = 0.0, 2.0 */
}
if (x >= 0) return q_gamma(0.5, x * x, LOG_PI_OVER_2);
else return 1 + p_gamma(0.5, x * x, LOG_PI_OVER_2);
}
|
