blitz Version 0.10
random/normal.h
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 // $Id: normal.h,v 1.5 2011/03/25 22:41:17 julianc Exp $
00003 
00004 /*
00005  * This is a modification of the Kinderman + Monahan algorithm for
00006  * generating normal random numbers, due to Leva:
00007  *
00008  * J.L. Leva, Algorithm 712. A normal random number generator, ACM Trans. 
00009  * Math. Softw.  18 (1992) 454--455. 
00010  *
00011  * http://www.acm.org/pubs/citations/journals/toms/1992-18-4/p449-leva/
00012  *
00013  * Note: Some of the constants used below look like they have dubious
00014  * precision.  These constants are used for an approximate bounding 
00015  * region test (see the paper).  If the approximate test fails,
00016  * then an exact region test is performed.
00017  *
00018  * Only 0.012 logarithm evaluations are required per random number
00019  * generated, making this method comparatively fast.
00020  *
00021  * Adapted to C++ by T. Veldhuizen.
00022  */
00023 
00024 #ifndef BZ_RANDOM_NORMAL
00025 #define BZ_RANDOM_NORMAL
00026 
00027 #ifndef BZ_RANDOM_UNIFORM
00028  #include <random/uniform.h>
00029 #endif
00030 
00031 BZ_NAMESPACE(ranlib)
00032 
00033 template<typename T = double, typename IRNG = defaultIRNG, 
00034     typename stateTag = defaultState>
00035 class NormalUnit : public UniformOpen<T,IRNG,stateTag>
00036 {
00037 public:
00038     typedef T T_numtype;
00039 
00040   NormalUnit() {}
00041 
00042   explicit NormalUnit(unsigned int i) :
00043     UniformOpen<T,IRNG,stateTag>(i) {};
00044 
00045     T random()
00046     {
00047         const T s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472;
00048         const T r1 = 0.27597, r2 = 0.27846;
00049 
00050         T u, v;
00051 
00052         for (;;) {
00053           // Generate P = (u,v) uniform in rectangle enclosing
00054           // acceptance region:
00055           //   0 < u < 1
00056           // - sqrt(2/e) < v < sqrt(2/e)
00057           // The constant below is 2*sqrt(2/e).
00058 
00059           u = this->getUniform();
00060           v = 1.715527769921413592960379282557544956242L 
00061               * (this->getUniform() - 0.5);
00062 
00063           // Evaluate the quadratic form
00064           T x = u - s;
00065           T y = fabs(v) - t;
00066           T q = x*x + y*(a*y - b*x);
00067        
00068           // Accept P if inside inner ellipse
00069           if (q < r1)
00070             break;
00071 
00072           // Reject P if outside outer ellipse
00073           if (q > r2)
00074             continue;
00075 
00076           // Between ellipses: perform exact test
00077           if (v*v <= -4.0 * log(u)*u*u)
00078             break;
00079         }
00080 
00081         return v/u;
00082     }
00083 
00084 };
00085 
00086 
00087 template<typename T = double, typename IRNG = defaultIRNG, 
00088     typename stateTag = defaultState>
00089 class Normal : public NormalUnit<T,IRNG,stateTag> {
00090 
00091 public:
00092     typedef T T_numtype;
00093 
00094     Normal(T mean, T standardDeviation)
00095     {
00096         mean_ = mean;
00097         standardDeviation_ = standardDeviation;
00098     }
00099 
00100   Normal(T mean, T standardDeviation, unsigned int i) :
00101     NormalUnit<T,IRNG,stateTag>(i) 
00102   {
00103         mean_ = mean;
00104         standardDeviation_ = standardDeviation;
00105   };
00106   
00107     T random()
00108     {
00109         return mean_ + standardDeviation_ 
00110            * NormalUnit<T,IRNG,stateTag>::random();
00111     }
00112 
00113 private:
00114     T mean_;
00115     T standardDeviation_;
00116 };
00117 
00118 BZ_NAMESPACE_END
00119 
00120 #endif // BZ_RANDOM_NORMAL
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines