normal.h

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

Generated on Mon Oct 22 10:04:52 2007 for blitz by  doxygen 1.3.9.1