forked from yixuan/RcppNumerical
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathErf.c
More file actions
51 lines (47 loc) · 1.39 KB
/
Erf.c
File metadata and controls
51 lines (47 loc) · 1.39 KB
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
/*
Erf.c
Gaussian error function
= 2/Sqrt[Pi] Integrate[Exp[-t^2], {t, 0, x}]
Code from Takuya Ooura's gamerf2a.f
http://www.kurims.kyoto-u.ac.jp/~ooura/gamerf.html
last modified 12 Mar 15 th
*/
static real Erfc(creal x)
{
static creal c[] = {
2.96316885199227378e-01, 6.12158644495538758e-02,
1.81581125134637070e-01, 5.50942780056002085e-01,
6.81866451424939493e-02, 1.53039662058770397e+00,
1.56907543161966709e-02, 2.99957952311300634e+00,
2.21290116681517573e-03, 4.95867777128246701e+00,
1.91395813098742864e-04, 7.41471251099335407e+00,
9.71013284010551623e-06, 1.04765104356545238e+01,
1.66642447174307753e-07, 1.48455557345597957e+01,
6.10399733098688199e+00, 1.26974899965115684e+01 };
real y = x*x;
y = expx(-y)*x*(
c[0]/(y + c[1]) + c[2]/(y + c[3]) +
c[4]/(y + c[5]) + c[6]/(y + c[7]) +
c[8]/(y + c[9]) + c[10]/(y + c[11]) +
c[12]/(y + c[13]) + c[14]/(y + c[15]) );
if( x < c[16] ) y += 2/(expx(c[17]*x) + 1);
return y;
}
static real Erf(creal x)
{
static creal c[] = {
1.12837916709551257e+00,
-3.76126389031833602e-01,
1.12837916706621301e-01,
-2.68661698447642378e-02,
5.22387877685618101e-03,
-8.49202435186918470e-04 };
real y = fabsx(x);
if( y > .125 ) {
y = 1 - Erfc(y);
return (x > 0) ? y : -y;
}
y *= y;
return x*(c[0] + y*(c[1] + y*(c[2] +
y*(c[3] + y*(c[4] + y*c[5])))));
}