home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / diverses / jpl_c / randnorm.c < prev    next >
Encoding:
C/C++ Source or Header  |  1986-05-18  |  1.6 KB  |  56 lines
  1. /* 1.0  11-12-84 */
  2. /************************************************************************
  3.  *            Robert C. Tausworthe                *
  4.  *            Jet Propulsion Laboratory            *
  5.  *            Pasadena, CA 91009        1984        *
  6.  ************************************************************************
  7.  *    Normal-Distribution Random Number Generator.
  8.  *
  9.  * Computes a sample value from a normal distribution characterized by
  10.  * zero mean and unit standard deviation.  Uses "Box-Muller polar method"
  11.  * for calculation of normal deviates,
  12.  *
  13.  * Box, G. E. P., and Muller, M. E., "A Note on the Generation of Normal
  14.  * Deviates," ANNALS OF MATH. STAT., Vol. 29, pp. 610-611, 1958.
  15.  *
  16.  * Refined in 
  17.  *
  18.  * Knuth, D. E., FUNDAMENTAL ALGORITHMS, Vol II, "Seminumerical Algorithms,"
  19.  * Addison-Wesley Pub. Co., page 117.
  20.  *
  21.  */
  22.  
  23. #include "defs.h"
  24. #include "stdtyp.h"
  25. #include "mathtyp.h"
  26.  
  27. /************************************************************************/
  28.     double
  29. randnorm()        /* Return a normally distributed random value
  30.                with zero mean and unit standard deviation.    */
  31. /*----------------------------------------------------------------------*/
  32. {
  33.     LOCAL BOOL first = TRUE;
  34.     double d, v1, v2;
  35.  
  36.     if (first--)
  37.     {    do
  38.         {    v1 = 2.0 * random() - 1.0;
  39.             v2 = 2.0 * random() - 1.0;
  40.             /* v1 and v2 are each uniform on [-1, 1] */
  41.             d = v1 * v1 + v2 * v2;
  42.         } while ((d <= 0.0) OR (d >= 1.0));
  43.         /* now, the point (v1, v2) is uniform inside the unit circle,
  44.          * excluding the origin.
  45.          */
  46.         d = sqrt(-2.0 * log(d) / d);
  47.         v1 *= d;
  48.         v2 *= d;
  49.         return (v1);
  50.     }
  51.     else
  52.     {    first = TRUE;
  53.         return (v2);
  54.     }
  55. }
  56.