blitz  Version 0.9
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
normal.h
Go to the documentation of this file.
1 /*
2  * This is a modification of the Kinderman + Monahan algorithm for
3  * generating normal random numbers, due to Leva:
4  *
5  * J.L. Leva, Algorithm 712. A normal random number generator, ACM Trans.
6  * Math. Softw. 18 (1992) 454--455.
7  *
8  * http://www.acm.org/pubs/citations/journals/toms/1992-18-4/p449-leva/
9  *
10  * Note: Some of the constants used below look like they have dubious
11  * precision. These constants are used for an approximate bounding
12  * region test (see the paper). If the approximate test fails,
13  * then an exact region test is performed.
14  *
15  * Only 0.012 logarithm evaluations are required per random number
16  * generated, making this method comparatively fast.
17  *
18  * Adapted to C++ by T. Veldhuizen.
19  */
20 
21 #ifndef BZ_RANDOM_NORMAL
22 #define BZ_RANDOM_NORMAL
23 
24 #ifndef BZ_RANDOM_UNIFORM
25  #include <random/uniform.h>
26 #endif
27 
28 BZ_NAMESPACE(ranlib)
29 
30 template<typename T = double, typename IRNG = defaultIRNG,
31  typename stateTag = defaultState>
32 class NormalUnit : public UniformOpen<T,IRNG,stateTag>
33 {
34 public:
35  typedef T T_numtype;
36 
37  T random()
38  {
39  const T s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472;
40  const T r1 = 0.27597, r2 = 0.27846;
41 
42  T u, v;
43 
44  for (;;) {
45  // Generate P = (u,v) uniform in rectangle enclosing
46  // acceptance region:
47  // 0 < u < 1
48  // - sqrt(2/e) < v < sqrt(2/e)
49  // The constant below is 2*sqrt(2/e).
50 
51  u = this->getUniform();
52  v = 1.715527769921413592960379282557544956242L
53  * (this->getUniform() - 0.5);
54 
55  // Evaluate the quadratic form
56  T x = u - s;
57  T y = fabs(v) - t;
58  T q = x*x + y*(a*y - b*x);
59 
60  // Accept P if inside inner ellipse
61  if (q < r1)
62  break;
63 
64  // Reject P if outside outer ellipse
65  if (q > r2)
66  continue;
67 
68  // Between ellipses: perform exact test
69  if (v*v <= -4.0 * log(u)*u*u)
70  break;
71  }
72 
73  return v/u;
74  }
75 
76 };
77 
78 
79 template<typename T = double, typename IRNG = defaultIRNG,
80  typename stateTag = defaultState>
81 class Normal : public NormalUnit<T,IRNG,stateTag> {
82 
83 public:
84  typedef T T_numtype;
85 
86  Normal(T mean, T standardDeviation)
87  {
88  mean_ = mean;
89  standardDeviation_ = standardDeviation;
90  }
91 
92  T random()
93  {
94  return mean_ + standardDeviation_
96  }
97 
98 private:
99  T mean_;
101 };
102 
104 
105 #endif // BZ_RANDOM_NORMAL