From f6b8e0eae4374f339487a33e3e4fe5462d5816e1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B6rg=20Frings-F=C3=BCrst?= Date: Sat, 25 Nov 2017 10:16:00 +0100 Subject: New upstream version 2.0.0 --- numlib/rand.c | 128 +++++++++++++++++++++++++++++++++++++++------------------- 1 file changed, 86 insertions(+), 42 deletions(-) mode change 100644 => 100755 numlib/rand.c (limited to 'numlib/rand.c') diff --git a/numlib/rand.c b/numlib/rand.c old mode 100644 new mode 100755 index 436de52..9e58de3 --- a/numlib/rand.c +++ b/numlib/rand.c @@ -19,84 +19,128 @@ /* (From Knuth & H.W.Lewis) */ #define PSRAND32L(S) ((S) * 1664525L + 1013904223L) +/* - - - - - - - - - - - - - - - */ +/* Global state random generator */ + +static rand_state g_rand = { 0 }; /* Use default seed for global state generator */ + /* Return a 32 bit number between 0 and 4294967295 */ /* Use Knuth shuffle to improve PSRAND32 sequence */ unsigned int rand32( /* Return 32 bit random number */ unsigned int seed /* Optional seed. Non-zero re-initialized with that seed */ ) { -#define TSIZE 2843 /* Prime */ - static unsigned int last, ran = 0x12345678; /* Default seed */ - static unsigned int pvs[TSIZE]; - static int pvs_inited = 0; - int i; - - if (seed != 0) - { - pvs_inited = 0; - ran = seed;; - } - - /* Init random storage locations */ - if (pvs_inited == 0) - { - for (i = 0; i < TSIZE; i++) - pvs[i] = ran = PSRAND32(ran); - last = ran; - pvs_inited = 1; - } - i = last % TSIZE; /* New location */ - last = pvs[i]; /* Value generated */ - pvs[i] = ran = PSRAND32(ran); /* New value */ - - return last-1; + return rand32_th(NULL, seed); } /* return a random number between 0.0 and 1.0 */ /* based on rand32 */ double ranno(void) { - return rand32(0) / 4294967295.0; + return rand32_th(NULL, 0) / 4294967295.0; } /* Return a random double in the range min to max */ double d_rand(double min, double max) { - double tt; - tt = ranno(); - return min + (max - min) * tt; + return d_rand_th(NULL, min, max); } /* Return a random integer in the range min to max inclusive */ int i_rand(int min, int max) { - double tt; - tt = ranno(); - return min + (int)floor(0.5 + ((double)(max - min)) * tt); + return i_rand_th(NULL, min, max); } - /* Return a random floating point number with a gausian/normal */ /* distribution, centered about 0.0, with standard deviation 1.0 */ /* This uses the Box-Muller transformation */ double norm_rand(void) { - static int r2 = 0; /* Can use 2nd number generated */ - static double nr2; + return norm_rand_th(NULL); +} + +/* - - - - - - - - - - - - - - - */ +/* Explicit state random generator */ + +/* Init rand_state to default */ +void rand_init(rand_state *p) { + if (p == NULL) + p = &g_rand; + p->pvs_inited = 0; + p->ran = 0; +} + +/* Return a 32 bit number between 0 and 4294967295 */ +/* Use Knuth shuffle to improve PSRAND32 sequence */ +unsigned int +rand32_th(rand_state *p, +unsigned int seed /* Optional seed. Non-zero re-initialized with that seed */ +) { + int i; + + if (p == NULL) + p = &g_rand; + + if (seed != 0) { + p->pvs_inited = 0; + p->ran = seed; + } + + /* Init random storage locations */ + if (p->pvs_inited == 0) { + if (p->ran == 0) + p->ran = RAND_SEED; + for (i = 0; i < RAND_TSIZE; i++) + p->pvs[i] = p->ran = PSRAND32(p->ran); + p->last = p->ran; + p->pvs_inited = 1; + } + i = p->last % RAND_TSIZE; /* New location */ + p->last = p->pvs[i]; /* Value generated */ + p->pvs[i] = p->ran = PSRAND32(p->ran); /* New value */ + + return p->last-1; +} + +/* return a random number between 0.0 and 1.0 */ +/* based on rand32 */ +double ranno_th(rand_state *p) { + return rand32_th(NULL, 0) / 4294967295.0; +} + +/* Return a random double in the range min to max */ +double +d_rand_th(rand_state *p, double min, double max) { + return min + (max - min) * ranno_th(p); +} + +/* Return a random integer in the range min to max inclusive */ +int +i_rand_th(rand_state *p, int min, int max) { + return min + (int)floor(0.5 + ((double)(max - min)) * ranno_th(p)); +} + +/* Return a random floating point number with a gausian/normal */ +/* distribution, centered about 0.0, with standard deviation 1.0 */ +/* This uses the Box-Muller transformation */ +double norm_rand_th(rand_state *p) { + if (p == NULL) + p = &g_rand; - if (r2 == 0) { + if (p->r2 == 0) { double v1, v2, t1, t2, r1; do { - v1 = d_rand(-1.0, 1.0); - v2 = d_rand(-1.0, 1.0); + v1 = d_rand_th(p, -1.0, 1.0); + v2 = d_rand_th(p, -1.0, 1.0); t1 = v1 * v1 + v2 * v2; } while (t1 == 0.0 || t1 >= 1.0); t2 = sqrt(-2.0 * log(t1)/t1); - nr2 = v2 * t2; /* One for next time */ - r2 = 1; + p->nr2 = v2 * t2; /* One for next time */ + p->r2 = 1; r1 = v1 * t2; return r1; } else { - r2 = 0; - return nr2; + p->r2 = 0; + return p->nr2; } } -- cgit v1.2.3