GRASS Programmer's Manual  6.4.2(2012)
as66.c
Go to the documentation of this file.
00001 
00002 /*-Algorithm AS 66
00003  * The Normal Integral, by I. D. Hill, 1973.
00004  * Applied Statistics 22(3):424-427.
00005  *
00006  * Translation to C by James Darrell McCauley, mccauley@ecn.purdue.edu.
00007  *
00008  * Calculates the upper or lower tail area of the standardized normal
00009  * curve corresponding to any given argument.
00010  *
00011  * x - the argument value
00012  * upper:  1 -> the area from x to \infty
00013  *         0 -> the area from -\infty to x
00014  *
00015  * Notes:
00016  * The constant LTONE should be set to the value at which the
00017  * lower tail area becomes 1.0 to the accuracy of the machine.
00018  * LTONE=(n+9)/3 gives the required value accurately enough, for a
00019  * machine that produces n decimal digits in its real numbers.
00020  *
00021  * The constant UTZERO should be set to the value at which the upper
00022  * tail area becomes 0.0 to the accuracy of the machine. This may be
00023  * taken as the value such that exp(-0.5 * UTZERO * UTZERO) /
00024  * (UTZERO * sqrt(2*M_PI)) is just greater than the smallest allowable
00025  * real numbers.
00026  */
00027 
00028 #include <math.h>
00029 
00030 
00031 #define LTONE 7.0
00032 #define UTZERO 18.66
00033 
00034 
00035 double alnorm(double x, int upper)
00036 {
00037     double ret, z, y;
00038     int up;
00039 
00040     up = upper;
00041     z = x;
00042 
00043     if (x < 0.0) {
00044         up = up == 0 ? 1 : 0;
00045         z = -x;
00046     }
00047 
00048     if (!(z <= LTONE || (up == 1 && z <= UTZERO)))
00049         ret = 0.0;
00050     else {
00051         y = 0.5 * z * z;
00052         if (z <= 1.28)
00053             ret = 0.5 - z * (0.398942280444 - 0.399903438504 * y /
00054                              (y + 5.75885480458 - 29.8213557808 /
00055                               (y + 2.62433121679 + 48.6959930692 /
00056                                (y + 5.92885724438))));
00057         else
00058             ret = 0.398942280385 * exp(-y) /
00059                 (z - 3.8052e-8 + 1.00000615302 /
00060                  (z + 3.98064794e-4 + 1.98615381364 /
00061                   (z - 0.151679116635 + 5.29330324926 /
00062                    (z + 4.8385912808 - 15.1508972451 /
00063                     (z + 0.742380924027 + 30.789933034 /
00064                      (z + 3.99019417011))))));
00065     }
00066 
00067     if (up == 0)
00068         ret = 1.0 - ret;
00069 
00070     return ret;
00071 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines