GRASS GIS 8 Programmer's Manual  8.4.0dev(2024)-2653d7b383
as66.c
Go to the documentation of this file.
1 /*-Algorithm AS 66
2  * The Normal Integral, by I. D. Hill, 1973.
3  * Applied Statistics 22(3):424-427.
4  *
5  * Translation to C by James Darrell McCauley, mccauley@ecn.purdue.edu.
6  *
7  * Calculates the upper or lower tail area of the standardized normal
8  * curve corresponding to any given argument.
9  *
10  * x - the argument value
11  * upper: 1 -> the area from x to \infty
12  * 0 -> the area from -\infty to x
13  *
14  * Notes:
15  * The constant LTONE should be set to the value at which the
16  * lower tail area becomes 1.0 to the accuracy of the machine.
17  * LTONE=(n+9)/3 gives the required value accurately enough, for a
18  * machine that produces n decimal digits in its real numbers.
19  *
20  * The constant UTZERO should be set to the value at which the upper
21  * tail area becomes 0.0 to the accuracy of the machine. This may be
22  * taken as the value such that exp(-0.5 * UTZERO * UTZERO) /
23  * (UTZERO * sqrt(2*M_PI)) is just greater than the smallest allowable
24  * real numbers.
25  */
26 
27 #include <math.h>
28 
29 #define LTONE 7.0
30 #define UTZERO 18.66
31 
32 double Cdhc_alnorm(double x, int upper)
33 {
34  double ret, z, y;
35  int up;
36 
37  up = upper;
38  z = x;
39 
40  if (x < 0.0) {
41  up = up == 0 ? 1 : 0;
42  z = -x;
43  }
44 
45  if (!(z <= LTONE || (up == 1 && z <= UTZERO)))
46  ret = 0.0;
47  else {
48  y = 0.5 * z * z;
49  if (z <= 1.28)
50  ret = 0.5 - z * (0.398942280444 -
51  0.399903438504 * y /
52  (y + 5.75885480458 -
53  29.8213557808 /
54  (y + 2.62433121679 +
55  48.6959930692 / (y + 5.92885724438))));
56  else
57  ret = 0.398942280385 * exp(-y) /
58  (z - 3.8052e-8 +
59  1.00000615302 /
60  (z + 3.98064794e-4 +
61  1.98615381364 /
62  (z - 0.151679116635 +
63  5.29330324926 /
64  (z + 4.8385912808 -
65  15.1508972451 /
66  (z + 0.742380924027 +
67  30.789933034 / (z + 3.99019417011))))));
68  }
69 
70  if (up == 0)
71  ret = 1.0 - ret;
72 
73  return ret;
74 }
double Cdhc_alnorm(double x, int upper)
Definition: as66.c:32
#define UTZERO
Definition: as66.c:30
#define LTONE
Definition: as66.c:29
#define x