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
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
40
NormalUnit
() {}
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_
110
*
NormalUnit<T,IRNG,stateTag>::random
();
111
}
112
113
private
:
114
T
mean_
;
115
T
standardDeviation_
;
116
};
117
118
}
119
120
#endif
// BZ_RANDOM_NORMAL
ranlib::NormalUnit
Definition
normal.h:36
ranlib::NormalUnit::NormalUnit
NormalUnit()
Definition
normal.h:40
ranlib::NormalUnit::NormalUnit
NormalUnit(unsigned int i)
Definition
normal.h:42
ranlib::NormalUnit::random
T random()
Definition
normal.h:45
ranlib::NormalUnit::T_numtype
T T_numtype
Definition
normal.h:38
ranlib::Normal
Definition
normal.h:89
ranlib::Normal::T_numtype
T T_numtype
Definition
normal.h:92
ranlib::Normal::mean_
T mean_
Definition
normal.h:114
ranlib::Normal::Normal
Normal(T mean, T standardDeviation)
Definition
normal.h:94
ranlib::Normal::standardDeviation_
T standardDeviation_
Definition
normal.h:115
ranlib::Normal::Normal
Normal(T mean, T standardDeviation, unsigned int i)
Definition
normal.h:100
ranlib::Normal::random
T random()
Definition
normal.h:107
ranlib::UniformOpen
Definition
uniform.h:377
ranlib::UniformOpen::getUniform
T getUniform()
Definition
uniform.h:397
ranlib
Definition
beta.h:50
ranlib::defaultState
sharedState defaultState
Definition
default.h:55
ranlib::defaultIRNG
MersenneTwister defaultIRNG
Definition
default.h:120
uniform.h
random
normal.h
Generated by
1.10.0