GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
as241.c
Go to the documentation of this file.
1#include <math.h>
2
3/*-
4 * algorithm as241 appl. statist. (1988) 37(3):477-484.
5 * produces the normal deviate z corresponding to a given lower tail
6 * area of p; z is accurate to about 1 part in 10**7.
7 *
8 * the hash sums below are the sums of the mantissas of the coefficients.
9 * they are included for use in checking transcription.
10 */
11double Cdhc_ppnd7(double p)
12{
13 static double zero = 0.0, one = 1.0, half = 0.5;
14 static double split1 = 0.425, split2 = 5.0;
15 static double const1 = 0.180625, const2 = 1.6;
16
17 /* coefficients for p close to 0.5 */
18 static double a[4] = {3.3871327179, 5.0434271938e+01, 1.5929113202e+02,
19 5.9109374720e+01};
20 static double b[4] = {0.0, 1.7895169469e+01, 7.8757757664e+01,
21 6.7187563600e+01};
22
23 /* hash sum ab 32.3184577772 */
24 /* coefficients for p not close to 0, 0.5 or 1. */
25 static double c[4] = {1.4234372777e+00, 2.7568153900e+00, 1.3067284816e+00,
26 1.7023821103e-01};
27 static double d[3] = {0.0, 7.3700164250e-01, 1.2021132975e-01};
28
29 /* hash sum cd 15.7614929821 */
30 /* coefficients for p near 0 or 1. */
31 static double e[4] = {6.6579051150e+00, 3.0812263860e+00, 4.2868294337e-01,
32 1.7337203997e-02};
33 static double f[3] = {0.0, 2.4197894225e-01, 1.2258202635e-02};
34
35 /* hash sum ef 19.4052910204 */
36 double q, r, ret;
37
38 q = p - half;
39 if (fabs(q) <= split1) {
40 r = const1 - q * q;
41 ret = q * (((a[3] * r + a[2]) * r + a[1]) * r + a[0]) /
42 (((b[3] * r + b[2]) * r + b[1]) * r + one);
43
44 return ret;
45 ;
46 }
47 /* else */
48
49 if (q < zero)
50 r = p;
51 else
52 r = one - p;
53
54 if (r <= zero)
55 return zero;
56
57 r = sqrt(-log(r));
58 if (r <= split2) {
59 r = r - const2;
60 ret = (((c[3] * r + c[2]) * r + c[1]) * r + c[0]) /
61 ((d[2] * r + d[1]) * r + one);
62 }
63 else {
64 r = r - split2;
65 ret = (((e[3] * r + e[2]) * r + e[1]) * r + e[0]) /
66 ((f[2] * r + f[1]) * r + one);
67 }
68
69 if (q < zero)
70 ret = -ret;
71
72 return ret;
73 ;
74}
75
76/*-
77 * algorithm as241 appl. statist. (1988) 37(3):
78 *
79 * produces the normal deviate z corresponding to a given lower
80 * tail area of p; z is accurate to about 1 part in 10**16.
81 *
82 * the hash sums below are the sums of the mantissas of the
83 * coefficients. they are included for use in checking
84 * transcription.
85 */
86double ppnd16(double p)
87{
88 static double zero = 0.0, one = 1.0, half = 0.5;
89 static double split1 = 0.425, split2 = 5.0;
90 static double const1 = 0.180625, const2 = 1.6;
91
92 /* coefficients for p close to 0.5 */
93 static double a[8] = {3.3871328727963666080e0, 1.3314166789178437745e+2,
94 1.9715909503065514427e+3, 1.3731693765509461125e+4,
95 4.5921953931549871457e+4, 6.7265770927008700853e+4,
96 3.3430575583588128105e+4, 2.5090809287301226727e+3};
97 static double b[8] = {0.0,
98 4.2313330701600911252e+1,
99 6.8718700749205790830e+2,
100 5.3941960214247511077e+3,
101 2.1213794301586595867e+4,
102 3.9307895800092710610e+4,
103 2.8729085735721942674e+4,
104 5.2264952788528545610e+3};
105
106 /* hash sum ab 55.8831928806149014439 */
107 /* coefficients for p not close to 0, 0.5 or 1. */
108 static double c[8] = {1.42343711074968357734e0, 4.63033784615654529590e0,
109 5.76949722146069140550e0, 3.64784832476320460504e0,
110 1.27045825245236838258e0, 2.41780725177450611770e-1,
111 2.27238449892691845833e-2, 7.74545014278341407640e-4};
112 static double d[8] = {0.0,
113 2.05319162663775882187e0,
114 1.67638483018380384940e0,
115 6.89767334985100004550e-1,
116 1.48103976427480074590e-1,
117 1.51986665636164571966e-2,
118 5.47593808499534494600e-4,
119 1.05075007164441684324e-9};
120
121 /* hash sum cd 49.33206503301610289036 */
122 /* coefficients for p near 0 or 1. */
123 static double e[8] = {6.65790464350110377720e0, 5.46378491116411436990e0,
124 1.78482653991729133580e0, 2.96560571828504891230e-1,
125 2.65321895265761230930e-2, 1.24266094738807843860e-3,
126 2.71155556874348757815e-5, 2.01033439929228813265e-7};
127 static double f[8] = {0.0,
128 5.99832206555887937690e-1,
129 1.36929880922735805310e-1,
130 1.48753612908506148525e-2,
131 7.86869131145613259100e-4,
132 1.84631831751005468180e-5,
133 1.42151175831644588870e-7,
134 2.04426310338993978564e-15};
135
136 /* hash sum ef 47.52583317549289671629 */
137 double q, r, ret;
138
139 q = p - half;
140 if (fabs(q) <= split1) {
141 r = const1 - q * q;
142 ret = q *
143 (((((((a[7] * r + a[6]) * r + a[5]) * r + a[4]) * r + a[3]) * r +
144 a[2]) *
145 r +
146 a[1]) *
147 r +
148 a[0]) /
149 (((((((b[7] * r + b[6]) * r + b[5]) * r + b[4]) * r + b[3]) * r +
150 b[2]) *
151 r +
152 b[1]) *
153 r +
154 one);
155
156 return ret;
157 }
158 /* else */
159
160 if (q < zero)
161 r = p;
162 else
163 r = one - p;
164
165 if (r <= zero)
166 return zero;
167
168 r = sqrt(-log(r));
169 if (r <= split2) {
170 r -= const2;
171 ret = (((((((c[7] * r + c[6]) * r + c[5]) * r + c[4]) * r + c[3]) * r +
172 c[2]) *
173 r +
174 c[1]) *
175 r +
176 c[0]) /
177 (((((((d[7] * r + d[6]) * r + d[5]) * r + d[4]) * r + d[3]) * r +
178 d[2]) *
179 r +
180 d[1]) *
181 r +
182 one);
183 }
184 else {
185 r -= split2;
186 ret = (((((((e[7] * r + e[6]) * r + e[5]) * r + e[4]) * r + e[3]) * r +
187 e[2]) *
188 r +
189 e[1]) *
190 r +
191 e[0]) /
192 (((((((f[7] * r + f[6]) * r + f[5]) * r + f[4]) * r + f[3]) * r +
193 f[2]) *
194 r +
195 f[1]) *
196 r +
197 one);
198 }
199
200 if (q < zero)
201 ret = -ret;
202
203 return ret;
204}
double ppnd16(double p)
Definition as241.c:86
double Cdhc_ppnd7(double p)
Definition as241.c:11
double b
double r