blitz  Version 0.9
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
rand-normal.h
Go to the documentation of this file.
1 /***************************************************************************
2  * blitz/rand-normal.h Random Gaussian (Normal) generator
3  *
4  * $Id: rand-normal.h,v 1.4 2003/12/11 03:44:22 julianc Exp $
5  *
6  * Copyright (C) 1997-2001 Todd Veldhuizen <tveldhui@oonumerics.org>
7  *
8  * This program is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU General Public License
10  * as published by the Free Software Foundation; either version 2
11  * of the License, or (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * Suggestions: blitz-dev@oonumerics.org
19  * Bugs: blitz-bugs@oonumerics.org
20  *
21  * For more information, please see the Blitz++ Home Page:
22  * http://oonumerics.org/blitz/
23  *
24  ***************************************************************************
25  *
26  * This generator transforms a (0,1] uniform distribution into
27  * a Normal distribution. Let u,v be (0,1] random variables. Then
28  *
29  * x = sqrt(-2 ln v) cos(pi (2u-1))
30  *
31  * is N(0,1) distributed.
32  *
33  * Reference: Athanasios Papoulis, "Probability, random variables,
34  * and stochastic processes," McGraw-Hill : Toronto, 1991.
35  *
36  ***************************************************************************/
37 
38 #ifndef BZ_RAND_NORMAL_H
39 #define BZ_RAND_NORMAL_H
40 
41 #ifndef BZ_RANDOM_H
42  #include <blitz/random.h>
43 #endif
44 
45 #ifndef BZ_RAND_UNIFORM_H
46  #include <blitz/rand-uniform.h>
47 #endif
48 
49 #include <math.h>
50 
51 BZ_NAMESPACE(blitz)
52 
53 template<typename P_uniform BZ_TEMPLATE_DEFAULT(Uniform)>
54 class Normal {
55 
56 public:
57  typedef double T_numtype;
58 
59  Normal(double mean = 0.0, double variance = 1.0, double = 0.0)
60  : mean_(mean), sigma_(::sqrt(variance))
61  {
62  }
63 
64  void randomize()
65  {
66  uniform_.randomize();
67  }
68 
69  double random()
70  {
71  double u, v;
72 
73  do {
74  u = uniform_.random();
75  v = uniform_.random();
76  } while (v == 0);
77 
78  return mean_ + sigma_ * ::sqrt(-2*::log(v)) * ::cos(M_PI * (2*u - 1));
79  }
80 
81 private:
82  double mean_, sigma_;
83  P_uniform uniform_;
84 };
85 
87 
88 #endif // BZ_RAND_NORMAL_H
89