summaryrefslogtreecommitdiff
path: root/numlib/rand.c
diff options
context:
space:
mode:
Diffstat (limited to 'numlib/rand.c')
-rwxr-xr-x[-rw-r--r--]numlib/rand.c128
1 files changed, 86 insertions, 42 deletions
diff --git a/numlib/rand.c b/numlib/rand.c
index 436de52..9e58de3 100644..100755
--- 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;
}
}