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