GRASS Programmer's Manual
6.4.2(2012)
|
00001 #include <math.h> 00002 00003 00004 /*- 00005 * algorithm as241 appl. statist. (1988) 37(3):477-484. 00006 * produces the normal deviate z corresponding to a given lower tail 00007 * area of p; z is accurate to about 1 part in 10**7. 00008 * 00009 * the hash sums below are the sums of the mantissas of the coefficients. 00010 * they are included for use in checking transcription. 00011 */ 00012 double ppnd7(double p) 00013 { 00014 static double zero = 0.0, one = 1.0, half = 0.5; 00015 static double split1 = 0.425, split2 = 5.0; 00016 static double const1 = 0.180625, const2 = 1.6; 00017 00018 /* coefficients for p close to 0.5 */ 00019 static double a[4] = { 3.3871327179, 5.0434271938e+01, 00020 1.5929113202e+02, 5.9109374720e+01 00021 }; 00022 static double b[4] = { 0.0, 1.7895169469e+01, 7.8757757664e+01, 00023 6.7187563600e+01 00024 }; 00025 00026 /* hash sum ab 32.3184577772 */ 00027 /* coefficients for p not close to 0, 0.5 or 1. */ 00028 static double c[4] = { 1.4234372777e+00, 2.7568153900e+00, 00029 1.3067284816e+00, 1.7023821103e-01 00030 }; 00031 static double d[3] = { 0.0, 7.3700164250e-01, 1.2021132975e-01 }; 00032 00033 /* hash sum cd 15.7614929821 */ 00034 /* coefficients for p near 0 or 1. */ 00035 static double e[4] = { 6.6579051150e+00, 3.0812263860e+00, 00036 4.2868294337e-01, 1.7337203997e-02 00037 }; 00038 static double f[3] = { 0.0, 2.4197894225e-01, 1.2258202635e-02 }; 00039 00040 /* hash sum ef 19.4052910204 */ 00041 double q, r, ret; 00042 00043 q = p - half; 00044 if (fabs(q) <= split1) { 00045 r = const1 - q * q; 00046 ret = q * (((a[3] * r + a[2]) * r + a[1]) * r + a[0]) / 00047 (((b[3] * r + b[2]) * r + b[1]) * r + one); 00048 00049 return ret;; 00050 } 00051 /* else */ 00052 00053 if (q < zero) 00054 r = p; 00055 else 00056 r = one - p; 00057 00058 if (r <= zero) 00059 return zero; 00060 00061 r = sqrt(-log(r)); 00062 if (r <= split2) { 00063 r = r - const2; 00064 ret = (((c[3] * r + c[2]) * r + c[1]) * r + c[0]) / 00065 ((d[2] * r + d[1]) * r + one); 00066 } 00067 else { 00068 r = r - split2; 00069 ret = (((e[3] * r + e[2]) * r + e[1]) * r + e[0]) / 00070 ((f[2] * r + f[1]) * r + one); 00071 } 00072 00073 if (q < zero) 00074 ret = -ret; 00075 00076 return ret;; 00077 } 00078 00079 00080 /*- 00081 * algorithm as241 appl. statist. (1988) 37(3): 00082 * 00083 * produces the normal deviate z corresponding to a given lower 00084 * tail area of p; z is accurate to about 1 part in 10**16. 00085 * 00086 * the hash sums below are the sums of the mantissas of the 00087 * coefficients. they are included for use in checking 00088 * transcription. 00089 */ 00090 double ppnd16(double p) 00091 { 00092 static double zero = 0.0, one = 1.0, half = 0.5; 00093 static double split1 = 0.425, split2 = 5.0; 00094 static double const1 = 0.180625, const2 = 1.6; 00095 00096 /* coefficients for p close to 0.5 */ 00097 static double a[8] = { 00098 3.3871328727963666080e0, 00099 1.3314166789178437745e+2, 00100 1.9715909503065514427e+3, 00101 1.3731693765509461125e+4, 00102 4.5921953931549871457e+4, 00103 6.7265770927008700853e+4, 00104 3.3430575583588128105e+4, 00105 2.5090809287301226727e+3 00106 }; 00107 static double b[8] = { 0.0, 00108 4.2313330701600911252e+1, 00109 6.8718700749205790830e+2, 00110 5.3941960214247511077e+3, 00111 2.1213794301586595867e+4, 00112 3.9307895800092710610e+4, 00113 2.8729085735721942674e+4, 00114 5.2264952788528545610e+3 00115 }; 00116 00117 /* hash sum ab 55.8831928806149014439 */ 00118 /* coefficients for p not close to 0, 0.5 or 1. */ 00119 static double c[8] = { 00120 1.42343711074968357734e0, 00121 4.63033784615654529590e0, 00122 5.76949722146069140550e0, 00123 3.64784832476320460504e0, 00124 1.27045825245236838258e0, 00125 2.41780725177450611770e-1, 00126 2.27238449892691845833e-2, 00127 7.74545014278341407640e-4 00128 }; 00129 static double d[8] = { 0.0, 00130 2.05319162663775882187e0, 00131 1.67638483018380384940e0, 00132 6.89767334985100004550e-1, 00133 1.48103976427480074590e-1, 00134 1.51986665636164571966e-2, 00135 5.47593808499534494600e-4, 00136 1.05075007164441684324e-9 00137 }; 00138 00139 /* hash sum cd 49.33206503301610289036 */ 00140 /* coefficients for p near 0 or 1. */ 00141 static double e[8] = { 00142 6.65790464350110377720e0, 00143 5.46378491116411436990e0, 00144 1.78482653991729133580e0, 00145 2.96560571828504891230e-1, 00146 2.65321895265761230930e-2, 00147 1.24266094738807843860e-3, 00148 2.71155556874348757815e-5, 00149 2.01033439929228813265e-7 00150 }; 00151 static double f[8] = { 0.0, 00152 5.99832206555887937690e-1, 00153 1.36929880922735805310e-1, 00154 1.48753612908506148525e-2, 00155 7.86869131145613259100e-4, 00156 1.84631831751005468180e-5, 00157 1.42151175831644588870e-7, 00158 2.04426310338993978564e-15 00159 }; 00160 00161 /* hash sum ef 47.52583317549289671629 */ 00162 double q, r, ret; 00163 00164 q = p - half; 00165 if (fabs(q) <= split1) { 00166 r = const1 - q * q; 00167 ret = q * (((((((a[7] * r + a[6]) * r + a[5]) * r + a[4]) * r + a[3]) 00168 * r + a[2]) * r + a[1]) * r + a[0]) / 00169 (((((((b[7] * r + b[6]) * r + b[5]) * r + b[4]) * r + b[3]) 00170 * r + b[2]) * r + b[1]) * r + one); 00171 00172 return ret; 00173 } 00174 /* else */ 00175 00176 if (q < zero) 00177 r = p; 00178 else 00179 r = one - p; 00180 00181 if (r <= zero) 00182 return zero; 00183 00184 r = sqrt(-log(r)); 00185 if (r <= split2) { 00186 r -= const2; 00187 ret = (((((((c[7] * r + c[6]) * r + c[5]) * r + c[4]) * r + c[3]) 00188 * r + c[2]) * r + c[1]) * r + c[0]) / 00189 (((((((d[7] * r + d[6]) * r + d[5]) * r + d[4]) * r + d[3]) 00190 * r + d[2]) * r + d[1]) * r + one); 00191 } 00192 else { 00193 r -= split2; 00194 ret = (((((((e[7] * r + e[6]) * r + e[5]) * r + e[4]) * r + e[3]) 00195 * r + e[2]) * r + e[1]) * r + e[0]) / 00196 (((((((f[7] * r + f[6]) * r + f[5]) * r + f[4]) * r + f[3]) 00197 * r + f[2]) * r + f[1]) * r + one); 00198 } 00199 00200 if (q < zero) 00201 ret = -ret; 00202 00203 return ret; 00204 }