GRASS Programmer's Manual  6.4.2(2012)
normp.c
Go to the documentation of this file.
00001 #include <math.h>
00002 
00003 
00004 /*-
00005  *     SUBROUTINE NORMP(Z, P, Q, PDF)
00006  *
00007  *     Normal distribution probabilities accurate to 1.e-15.
00008  *     Z = no. of standard deviations from the mean.
00009  *     P, Q = probabilities to the left & right of Z.   P + Q = 1.
00010  *     PDF = the probability density.
00011  *
00012  *     Based upon algorithm 5666 for the error function, from:
00013  *     Hart, J.F. et al, 'Computer Approximations', Wiley 1968
00014  *
00015  *     Programmer: Alan Miller
00016  *
00017  *     Latest revision - 30 March 1986
00018  *
00019  */
00020 
00021 /* Conversion to C by James Darrell McCauley, 24 September 1994 */
00022 
00023 double normp(double z)
00024 {
00025     static double p[7] = { 220.2068679123761, 221.2135961699311,
00026         112.079291497870, 33.91286607838300, 6.37396220353165,
00027         0.7003830644436881, 0.352624965998910e-1
00028     };
00029     static double q[8] = { 440.4137358247522, 793.8265125199484,
00030         637.3336333788311, 296.5642487796737, 86.78073220294608,
00031         16.06417757920695, 1.755667163182642, 0.8838834764831844e-1
00032     };
00033     static double cutoff = 7.071, root2pi = 2.506628274631001;
00034     double zabs, expntl;
00035     double pee, queue, pdf;
00036 
00037     zabs = fabs(z);
00038 
00039     if (zabs > 37.0) {
00040         pdf = 0.0;
00041         if (z > 0.0) {
00042             pee = 1.0;
00043             queue = 0.0;
00044         }
00045         else {
00046             pee = 0.0;
00047             queue = 1.0;
00048         }
00049 
00050         return pee;
00051     }
00052 
00053     expntl = exp(-0.5 * zabs * zabs);
00054     pdf = expntl / root2pi;
00055 
00056     if (zabs < cutoff)
00057         pee = expntl * ((((((p[6] * zabs + p[5]) * zabs + p[4])
00058                            * zabs + p[3]) * zabs + p[2]) * zabs +
00059                          p[1]) * zabs + p[0])
00060             / (((((((q[7] * zabs + q[6]) * zabs + q[5]) * zabs + q[4])
00061                   * zabs + q[3]) * zabs + q[2]) * zabs + q[1]) * zabs + q[0]);
00062     else
00063         pee = pdf / (zabs + 1.0 / (zabs + 2.0 / (zabs + 3.0 / (zabs + 4.0 /
00064                                                                (zabs +
00065                                                                 0.65)))));
00066 
00067     if (z < 0.0)
00068         queue = 1.0 - pee;
00069     else {
00070         queue = pee;
00071         pee = 1.0 - queue;
00072     }
00073 
00074     return pee;
00075 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines