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
|
/*
* Copyright (C) 2002 by Phil Karn KA9Q
* Copyright (C) 2003 by Jonathan Naylor G4KLX
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
#include "Gaussian.h"
#include <wx/datetime.h>
#include <cmath>
#include <limits.h>
using namespace std;
CGaussian::CGaussian() :
m_gset(0.0),
m_iset(false)
{
wxDateTime now = wxDateTime::UNow();
::srandom(now.GetMillisecond());
}
CGaussian::~CGaussian()
{
}
/*
* Generate gaussian random double with specified mean and std_dev
* using the rejection method
*/
double CGaussian::rand(double mean, double standardDeviation)
{
if (m_iset) {
/* Already got one */
m_iset = false;
return mean + standardDeviation * m_gset;
}
/*
* Generate two evenly distributed numbers between -1 and +1
* that are inside the unit circle
*/
double rsq, v1, v2;
do {
v1 = 2.0 * double(::random()) / double(RAND_MAX - 1);
v2 = 2.0 * double(::random()) / double(RAND_MAX - 1);
rsq = v1 * v1 + v2 * v2;
} while (rsq >= 1.0 || rsq == 0.0);
double fac = ::sqrt(-2.0 * ::log(rsq) / rsq);
m_gset = v1 * fac;
m_iset = true;
return mean + standardDeviation * v2 * fac;
}
|
