2010-08-08 10:06:34 +02:00
|
|
|
/*
|
|
|
|
Copyright (c) 1998-2010, Enno Rehling <enno@eressea.de>
|
|
|
|
Katja Zedel <katze@felidae.kn-bremen.de
|
|
|
|
Christian Schlittchen <corwin@amber.kn-bremen.de>
|
|
|
|
|
|
|
|
Permission to use, copy, modify, and/or distribute this software for any
|
|
|
|
purpose with or without fee is hereby granted, provided that the above
|
|
|
|
copyright notice and this permission notice appear in all copies.
|
|
|
|
|
|
|
|
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
|
|
|
|
WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
|
|
|
|
MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
|
|
|
|
ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
|
|
|
|
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
|
|
|
|
ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
|
|
|
|
OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
|
|
|
|
**/
|
|
|
|
|
|
|
|
#include <platform.h>
|
|
|
|
#include "rand.h"
|
|
|
|
#include "rng.h"
|
|
|
|
|
|
|
|
#include <assert.h>
|
|
|
|
#include <string.h>
|
|
|
|
#include <math.h>
|
|
|
|
#include <float.h>
|
|
|
|
#include <ctype.h>
|
|
|
|
|
2011-03-07 08:02:35 +01:00
|
|
|
#define M_PIl 3.1415926535897932384626433832795029L /* pi */
|
2010-08-08 10:06:34 +02:00
|
|
|
|
|
|
|
/* NormalRand aus python, random.py geklaut, dort ist Referenz auf
|
|
|
|
* den Algorithmus. mu = Mittelwert, sigma = Standardabweichung.
|
|
|
|
* http://de.wikipedia.org/wiki/Standardabweichung#Diskrete_Gleichverteilung.2C_W.C3.BCrfel
|
|
|
|
*/
|
2011-03-07 08:02:35 +01:00
|
|
|
double normalvariate(double mu, double sigma)
|
2010-08-08 10:06:34 +02:00
|
|
|
{
|
2011-03-07 08:02:35 +01:00
|
|
|
static const double NV_MAGICCONST = 1.7155277699214135; /* STATIC_CONST: a constant */
|
2010-08-08 10:06:34 +02:00
|
|
|
double z;
|
|
|
|
for (;;) {
|
|
|
|
double u1 = rng_double();
|
|
|
|
double u2 = 1.0 - rng_double();
|
2011-03-07 08:02:35 +01:00
|
|
|
z = NV_MAGICCONST * (u1 - 0.5) / u2;
|
|
|
|
if (z * z / 4.0 <= -log(u2)) {
|
2010-08-08 10:06:34 +02:00
|
|
|
break;
|
|
|
|
}
|
|
|
|
}
|
2011-03-07 08:02:35 +01:00
|
|
|
return mu + z * sigma;
|
2010-08-08 10:06:34 +02:00
|
|
|
}
|
|
|
|
|
2011-03-07 08:02:35 +01:00
|
|
|
int ntimespprob(int n, double p, double mod)
|
2010-08-08 10:06:34 +02:00
|
|
|
{
|
|
|
|
int count = 0;
|
|
|
|
int i;
|
|
|
|
|
2011-03-07 08:02:35 +01:00
|
|
|
for (i = 0; i < n && p > 0; i++) {
|
|
|
|
if (rng_double() < p) {
|
2010-08-08 10:06:34 +02:00
|
|
|
count++;
|
|
|
|
p += mod;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return count;
|
|
|
|
}
|
|
|
|
|
2012-06-24 07:41:07 +02:00
|
|
|
bool chance(double x)
|
2010-08-08 10:06:34 +02:00
|
|
|
{
|
2011-03-07 08:02:35 +01:00
|
|
|
if (x >= 1.0)
|
|
|
|
return true;
|
2010-08-08 10:06:34 +02:00
|
|
|
return rng_double() < x;
|
|
|
|
}
|