blitz Version 1.0.2
Loading...
Searching...
No Matches
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
31namespace ranlib {
32
33template<typename T = double, typename IRNG = defaultIRNG,
34 typename stateTag = defaultState>
35class NormalUnit : public UniformOpen<T,IRNG,stateTag>
36{
37public:
38 typedef T T_numtype;
39
41
42 explicit NormalUnit(unsigned int i) :
43 UniformOpen<T,IRNG,stateTag>(i) {};
44
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
87template<typename T = double, typename IRNG = defaultIRNG,
88 typename stateTag = defaultState>
89class Normal : public NormalUnit<T,IRNG,stateTag> {
90
91public:
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
112
113private:
116};
117
118}
119
120#endif // BZ_RANDOM_NORMAL
Definition normal.h:36
NormalUnit()
Definition normal.h:40
NormalUnit(unsigned int i)
Definition normal.h:42
T random()
Definition normal.h:45
T T_numtype
Definition normal.h:38
Definition normal.h:89
T mean_
Definition normal.h:114
Normal(T mean, T standardDeviation)
Definition normal.h:94
T standardDeviation_
Definition normal.h:115
T T_numtype
Definition normal.h:92
Normal(T mean, T standardDeviation, unsigned int i)
Definition normal.h:100
T random()
Definition normal.h:107
Definition uniform.h:377
Definition beta.h:50
sharedState defaultState
Definition default.h:55
MersenneTwister defaultIRNG
Definition default.h:120