VERB_code_2.2
2
Main Page
Related Pages
Namespaces
Classes
Files
File List
File Members
All
Classes
Namespaces
Files
Functions
Variables
Typedefs
Macros
Pages
erf.cpp
Go to the documentation of this file.
1
/**
2
* \file erf.cpp
3
* Error function.
4
*
5
* \author Numerical Recepies.
6
*/
7
8
#include <string>
9
#include "
erf.h
"
10
#include <iostream>
11
12
using namespace
std;
13
14
bool
str2bool
(
string
str);
15
string
bool2str
(
bool
b);
16
17
void
nrerror
(
const
char
* msg) {
18
cout << msg;
19
return
;
20
}
21
22
///Returns the value ln[?(xx)] for xx > 0.
23
double
gammln
(
double
xx)
24
{
25
//Internal arithmetic will be done in double precision, a nicety that you can omit if five-figure accuracy is good enough.
26
double
x,y,tmp,ser;
27
static
double
cof[6]={76.18009172947146,-86.50532032941677, 24.01409824083091,-1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5};
28
int
j;
29
y=x=xx;
30
tmp=x+5.5;
31
tmp -= (x+0.5)*log(tmp);
32
ser=1.000000000190015;
33
for
(j=0;j<=5;j++) ser += cof[j]/++y;
34
return
-tmp+log(2.5066282746310005*ser/x);
35
}
36
37
38
///Returns the incomplete gamma function P(a, x).
39
double
gammp
(
double
a,
double
x)
40
{
41
void
gcf
(
double
*gammcf,
double
a,
double
x,
double
*gln);
42
void
gser
(
double
*gamser,
double
a,
double
x,
double
*gln);
43
// void nrerror(char error_text[]);
44
double
gamser,gammcf,gln;
45
if
(x < 0.0 || a <= 0.0)
nrerror
(
"Invalid arguments in routine gammp"
);
46
if
(x < (a+1.0)) {
//Use the series representation.
47
gser
(&gamser,a,x,&gln);
48
return
gamser;
49
}
else
{
//Use the continued fraction representation
50
gcf
(&gammcf,a,x,&gln);
51
return
1.0-gammcf;
//and take its complement.
52
}
53
}
54
55
///Returns the incomplete gamma function Q(a, x) ? 1 ? P(a, x).
56
double
gammq
(
double
a,
double
x)
57
{
58
void
gcf
(
double
*gammcf,
double
a,
double
x,
double
*gln);
59
void
gser
(
double
*gamser,
double
a,
double
x,
double
*gln);
60
// void nrerror(char error_text[]);
61
double
gamser,gammcf,gln;
62
if
(x < 0.0 || a <= 0.0)
nrerror
(
"Invalid arguments in routine gammq"
);
63
if
(x < (a+1.0)) {
//Use the series representation
64
gser
(&gamser,a,x,&gln);
65
return
1.0-gamser;
//and take its complement.
66
}
else
{
//Use the continued fraction representation.
67
gcf
(&gammcf,a,x,&gln);
68
return
gammcf;
69
}
70
}
71
72
///Returns the incomplete gamma function P(a, x) evaluated by its series representation as gamser.
73
///Also returns ln ?(a) as gln.
74
void
gser
(
double
*gamser,
double
a,
double
x,
double
*gln)
75
{
76
double
gammln
(
double
xx);
77
// void nrerror(char error_text[]);
78
int
n;
79
double
sum,del,ap;
80
*gln=
gammln
(a);
81
if
(x <= 0.0) {
82
if
(x < 0.0)
nrerror
(
"x less than 0 in routine gser"
);
83
*gamser=0.0;
84
return
;
85
}
else
{
86
ap=a;
87
del=sum=1.0/a;
88
for
(n=1;n<=
ITMAX
;n++) {
89
++ap;
90
del *= x/ap;
91
sum += del;
92
if
(fabs(del) < fabs(sum)*
EPS
) {
93
*gamser=sum*
exp
(-x+a*log(x)-(*gln));
94
return
;
95
}
96
}
97
nrerror
(
"a too large, ITMAX too small in routine gser"
);
98
return
;
99
}
100
}
101
102
///Returns the incomplete gamma function Q(a, x) evaluated by its continued fraction representation as gammcf. Also returns ln?(a) as gln.
103
void
gcf
(
double
*gammcf,
double
a,
double
x,
double
*gln)
104
{
105
double
gammln
(
double
xx);
106
double
an;
107
double
b;
108
double
c1;
109
double
d;
110
double
del;
111
double
h;
112
// void nrerror(char error_text[]);
113
int
i;
114
*gln=
gammln
(a);
115
b=x+1.0-a;
//Set up for evaluating continued fraction by modified Lentz's method (5.2) with b0 = 0.
116
c1=1.0/
FPMIN
;
117
d=1.0/b;
118
h=d;
119
for
(i=1;i<=
ITMAX
;i++) {
//Iterate to convergence.
120
an = -i*(i-a);
121
b += 2.0;
122
d=an*d+b;
123
if
(fabs(d) <
FPMIN
) d=
FPMIN
;
124
c1=b+an/c1;
125
if
(fabs(c1) <
FPMIN
) c1=
FPMIN
;
126
d=1.0/d;
127
del=d*c1;
128
h *= del;
129
if
(fabs(del-1.0) <
EPS
)
break
;
130
}
131
if
(i >
ITMAX
)
nrerror
(
"a too large, ITMAX too small in gcf"
);
132
*gammcf=
exp
(-x+a*log(x)-(*gln))*h;
//Put factors in front.
133
}
134
135
136
///Returns the error function erf(x).
137
double
erf
(
double
x)
138
{
139
double
gammp
(
double
a,
double
x);
140
return
x < 0.0 ? -
gammp
(0.5,x*x) :
gammp
(0.5,x*x);
141
}
VariousFunctions
erf.cpp
Generated on Sat Nov 16 2013 09:31:46 for VERB_code_2.2 by
1.8.4