GRASS Programmer's Manual  6.4.2(2012)
as241.c
Go to the documentation of this file.
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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines