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