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