GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-d6dec75dd4
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  */
11 double 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  */
86 double 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 }
ppnd16
double ppnd16(double p)
Definition: as241.c:86
Cdhc_ppnd7
double Cdhc_ppnd7(double p)
Definition: as241.c:11
b
double b
Definition: r_raster.c:39
r
double r
Definition: r_raster.c:39