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