← Back to branch summary

~ubuntu-branches/ubuntu/vivid/aeolus/vivid

« back to all changes in this revision

Viewing changes to source/rngen.cc

  • Committer: Bazaar Package Importer
  • Author(s): Alessio Treglia
  • Date: 2010-04-19 19:12:51 UTC
  • mfrom: (1.1.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20100419191251-hgarjfcdfl7c0ryl
Tags: 0.8.4-3
http://bugs.debian.org/578278
debian/patches/01-makefile.patch: Drop -march=native flag, it isn't valid
for Debian packages as the results are unpredictable, thanks to
Bastian Blank for reporting this (Closes: #578278).

Show diffs side-by-side

added added

removed removed

Lines of Context:
source/rngen.cc
 
1
/*
 
2
  Copyright (C) 2003-2008 Fons Adriaensen <fons@kokkinizita.net>
 
3
    
 
4
  This program is free software; you can redistribute it and/or modify
 
5
  it under the terms of the GNU General Public License as published by
 
6
  the Free Software Foundation; either version 2 of the License, or
 
7
  (at your option) any later version.
 
8
 
 
9
  This program is distributed in the hope that it will be useful,
 
10
  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
  GNU General Public License for more details.
 
13
 
 
14
  You should have received a copy of the GNU General Public License
 
15
  along with this program; if not, write to the Free Software
 
16
  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 
17
*/
 
18
 
 
19
 
 
20
#include <math.h>
 
21
#include <time.h>
 
22
#include "rngen.h"
 
23
 
 
24
 
 
25
const double Rngen::_p31 = 2147483648.0;
 
26
const double Rngen::_p32 = 4294967296.0;
 
27
const float Rngen::_p31f = 2147483648.0f;
 
28
const float Rngen::_p32f = 4294967296.0f;
 
29
 
 
30
 
 
31
Rngen::Rngen (void)
 
32
{
 
33
    init (0);
 
34
}
 
35
 
 
36
 
 
37
Rngen::~Rngen (void)
 
38
{
 
39
}
 
40
 
 
41
 
 
42
void Rngen::init (uint32_t s)
 
43
{
 
44
    int i, j;
 
45
    Prbsgen G;
 
46
 
 
47
    if (s == 0) s = time (0);
 
48
 
 
49
    G.set_poly (Prbsgen::G32);
 
50
    G.set_stat (s);
 
51
 
 
52
    for (i = 0; i < 55; i++)
 
53
    {
 
54
        G.step ();
 
55
        j = G.stat () & 4095;
 
56
        while (j--) G.step ();
 
57
        _a [i] = G.stat ();
 
58
    }
 
59
 
 
60
    _i = 0;
 
61
    _md = false;
 
62
    _mf = false;
 
63
}
 
64
 
 
65
 
 
66
double Rngen::grand (void)
 
67
{
 
68
    double a, b, r;
 
69
 
 
70
    if (_md)
 
71
    {
 
72
        _md = false;
 
73
        return _vd;
 
74
    }
 
75
 
 
76
    do
 
77
    {
 
78
        a = irand () / _p31 - 1.0;
 
79
        b = irand () / _p31 - 1.0;
 
80
        r = a * a + b * b;
 
81
    }
 
82
    while ((r > 1.0) || (r < 1e-20));
 
83
 
 
84
    r = sqrt (-2.0 * log (r) / r);
 
85
    _md = true;
 
86
    _vd = r * b;
 
87
 
 
88
    return r * a;
 
89
}    
 
90
 
 
91
 
 
92
void Rngen::grand (double *x, double *y)
 
93
{
 
94
    double a, b, r;
 
95
 
 
96
    do
 
97
    {
 
98
        a = irand () / _p31 - 1.0;
 
99
        b = irand () / _p31 - 1.0;
 
100
        r = a * a + b * b;
 
101
    }
 
102
    while ((r > 1.0) || (r < 1e-20));
 
103
 
 
104
    r = sqrt (-log (r) / r);
 
105
    *x = r * a;
 
106
    *y = r * b;
 
107
}    
 
108
 
 
109
 
 
110
float Rngen::grandf (void)
 
111
{
 
112
    float a, b, r;
 
113
 
 
114
    if (_mf)
 
115
    {
 
116
        _mf = false;
 
117
        return _vf;
 
118
    }
 
119
 
 
120
    do
 
121
    {
 
122
        a = irand () / _p31f - 1.0f;
 
123
        b = irand () / _p31f - 1.0f;
 
124
        r = a * a + b * b;
 
125
    }
 
126
    while ((r > 1.0f) || (r < 1e-20f));
 
127
 
 
128
    r = sqrtf (-2.0f * logf (r) / r);
 
129
    _mf = true;
 
130
    _vf = r * b;
 
131
 
 
132
    return r * a;
 
133
}    
 
134
 
 
135
 
 
136
void Rngen::grandf (float *x, float *y)
 
137
{
 
138
    float a, b, r;
 
139
 
 
140
    do
 
141
    {
 
142
        a = irand () / _p31f - 1.0f;
 
143
        b = irand () / _p31f - 1.0f;
 
144
        r = a * a + b * b;
 
145
    }
 
146
    while ((r > 1.0f) || (r < 1e-20f));
 
147
 
 
148
    r = sqrtf (-logf (r) / r);
 
149
    *x = r * a;
 
150
    *y = r * b;
 
151
}