MuseScore/aeolus/rngen.cpp
2013-04-22 17:18:03 +02:00

150 lines
2.5 KiB
C++

/*
Copyright (C) 2003-2008 Fons Adriaensen <fons@kokkinizita.net>
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 <time.h>
#include "rngen.h"
const double Rngen::_p31 = 2147483648.0;
const double Rngen::_p32 = 4294967296.0;
const float Rngen::_p31f = 2147483648.0f;
const float Rngen::_p32f = 4294967296.0f;
Rngen::Rngen (void)
{
init (0);
}
Rngen::~Rngen (void)
{
}
void Rngen::init (uint32_t s)
{
int i, j;
Prbsgen G;
if (s == 0) s = time (0);
G.set_poly (Prbsgen::G32);
G.set_stat (s);
for (i = 0; i < 55; i++)
{
G.step ();
j = G.stat () & 4095;
while (j--) G.step ();
_a [i] = G.stat ();
}
_i = 0;
_md = false;
_mf = false;
}
double Rngen::grand (void)
{
double a, b, r;
if (_md)
{
_md = false;
return _vd;
}
do
{
a = irand () / _p31 - 1.0;
b = irand () / _p31 - 1.0;
r = a * a + b * b;
}
while ((r > 1.0) || (r < 1e-20));
r = sqrt (-2.0 * log (r) / r);
_md = true;
_vd = r * b;
return r * a;
}
void Rngen::grand (double *x, double *y)
{
double a, b, r;
do
{
a = irand () / _p31 - 1.0;
b = irand () / _p31 - 1.0;
r = a * a + b * b;
}
while ((r > 1.0) || (r < 1e-20));
r = sqrt (-log (r) / r);
*x = r * a;
*y = r * b;
}
float Rngen::grandf (void)
{
float a, b, r;
if (_mf)
{
_mf = false;
return _vf;
}
do
{
a = irand () / _p31f - 1.0f;
b = irand () / _p31f - 1.0f;
r = a * a + b * b;
}
while ((r > 1.0f) || (r < 1e-20f));
r = sqrtf (-2.0f * logf (r) / r);
_mf = true;
_vf = r * b;
return r * a;
}
void Rngen::grandf (float *x, float *y)
{
float a, b, r;
do
{
a = irand () / _p31f - 1.0f;
b = irand () / _p31f - 1.0f;
r = a * a + b * b;
}
while ((r > 1.0f) || (r < 1e-20f));
r = sqrtf (-logf (r) / r);
*x = r * a;
*y = r * b;
}