blitz Version 0.9
|
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