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