-
-
Notifications
You must be signed in to change notification settings - Fork 288
/
enormp.c
62 lines (57 loc) · 1.86 KB
/
enormp.c
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
#include <math.h>
double Cdhc_enormp(double x)
{
double x1, x2, x3, x4, ret_val;
static double xp[5] = { 7.7105849500132e-5, -0.00133733772997339,
0.0323076579225834, 0.0479137145607681, 0.128379167095513
};
static double xq[3] = { 0.00301048631703895, 0.0538971687740286,
0.375795757275549
};
static double xr[8] = { -1.36864857382717e-7, 0.564195517478974,
7.21175825088309, 43.1622272220567, 152.98928504694,
339.320816734344, 451.918953711873, 300.459261020162
};
static double xs[8] = { 1.0, 12.7827273196294, 77.0001529352295,
277.585444743988, 638.980264465631, 931.35409485061,
790.950925327898, 300.459260956983
};
static double xt[5] = { 2.10144126479064, 26.2370141675169,
21.3688200555087, 4.6580782871847, 0.282094791773523
};
static double xu[4] = { 94.153775055546, 187.11481179959,
99.0191814623914, 18.0124575948747
};
double yy1, yy2;
x3 = (double)0.564189583547756;
x1 = fabs(x);
if (x1 <= .5) {
x4 = x * x;
yy1 =
(((xp[0] * x4 + xp[1]) * x4 + xp[2]) * x4 + xp[3]) * x4 + xp[4] +
1.;
yy2 = ((xq[0] * x4 + xq[1]) * x4 + xq[2]) * x4 + 1.;
ret_val = x * (yy1 / yy2);
}
else if (x1 <= 4.) {
yy1 = ((((((xr[0] * x1 + xr[1]) * x1 + xr[2]) * x1 + xr[3])
* x1 + xr[4]) * x1 + xr[5]) * x1 + xr[6]) * x1 + xr[7];
yy2 = ((((((xs[0] * x1 + xs[1]) * x1 + xs[2]) * x1 + xs[3])
* x1 + xs[4]) * x1 + xs[5]) * x1 + xs[6]) * x1 + xs[7];
ret_val = 1. - exp(-x * x) * yy1 / yy2;
if (x < 0.)
ret_val *= -1.0;
}
else {
x2 = x * x;
x4 *= 1.; /* huh, what? See original FORTRAN */
x4 = 0.5;
yy1 = (((xt[0] * x4 + xt[1]) * x4 + xt[2]) * x4 + xt[3]) * x4 + xt[4];
yy2 = (((xu[0] * x4 + xu[1]) * x4 + xu[2]) * x4 + xu[3]) * x4 + 1.;
ret_val = x3 / x1 - yy1 * x1 / (x2 * yy2);
ret_val = 1. - exp(-x2) * ret_val;
if (x < 0.)
ret_val = -ret_val;
}
return ret_val;
}