guassian distribution has no endless loops.

This commit is contained in:
Enno Rehling 2017-05-03 21:02:30 +02:00
parent 5778bc2c93
commit 26795ae717
2 changed files with 19 additions and 14 deletions

View File

@ -14,6 +14,7 @@
#define NO_STRDUP #define NO_STRDUP
#define NO_MKDIR #define NO_MKDIR
#define _CRT_SECURE_NO_WARNINGS #define _CRT_SECURE_NO_WARNINGS
#define _USE_MATH_DEFINES
#pragma warning(disable: 4710 4820) #pragma warning(disable: 4710 4820)
#pragma warning(disable: 4100) // unreferenced formal parameter #pragma warning(disable: 4100) // unreferenced formal parameter
#pragma warning(disable: 4456) // declaration hides previous #pragma warning(disable: 4456) // declaration hides previous

View File

@ -36,23 +36,27 @@ int lovar(double xpct_x2)
return (rng_int() % n + rng_int() % n) / 1000; return (rng_int() % n + rng_int() % n) / 1000;
} }
/* NormalRand aus python, random.py geklaut, dort ist Referenz auf /* gaussian distribution
* den Algorithmus. mu = Mittelwert, sigma = Standardabweichung. * taken from http://c-faq.com/lib/gaussian.html
* http://de.wikipedia.org/wiki/Standardabweichung#Diskrete_Gleichverteilung.2C_W.C3.BCrfel */
*/
double normalvariate(double mu, double sigma) double normalvariate(double mu, double sigma)
{ {
static const double NV_MAGICCONST = 1.7155277699214135; static double U, V;
double z; static int phase = 0;
for (;;) { double Z;
double u1 = rng_double();
double u2 = rng_double(); if (phase == 0) {
z = NV_MAGICCONST * (u1 - 0.5) / u2; U = (rng_int() + 1.) / (RNG_RAND_MAX + 2.);
if (z * z / 4.0 <= -log(u2)) { V = rng_int() / (RNG_RAND_MAX + 1.);
break; Z = sqrt(-2 * log(U)) * sin(2 * M_PI * V);
} }
else {
Z = sqrt(-2 * log(U)) * cos(2 * M_PI * V);
} }
return mu + z * sigma; phase = 1 - phase;
return mu + Z *sigma;
} }
int ntimespprob(int n, double p, double mod) int ntimespprob(int n, double p, double mod)