#include #include "Random.h" #ifndef NULL #define NULL 0 #endif /* static const int mdig = 32; */ #define MDIG 32 /* static const int one = 1; */ #define ONE 1 static const int m1 = (ONE << (MDIG-2)) + ((ONE << (MDIG-2) )-ONE); static const int m2 = ONE << MDIG/2; /* For mdig = 32 : m1 = 2147483647, m2 = 65536 For mdig = 64 : m1 = 9223372036854775807, m2 = 4294967296 */ /* move to initialize() because */ /* compiler could not resolve as */ /* a constant. */ static /*const*/ double dm1; /* = 1.0 / (double) m1; */ /* private methods (defined below, but not in Random.h */ static void initialize(Random R, int seed); Random new_Random_seed(int seed) { Random R = (Random) malloc(sizeof(Random_struct)); initialize(R, seed); R->left = 0.0; R->right = 1.0; R->width = 1.0; R->haveRange = 0 /*false*/; return R; } Random new_Random(int seed, double left, double right) { Random R = (Random) malloc(sizeof(Random_struct)); initialize(R, seed); R->left = left; R->right = right; R->width = right - left; R->haveRange = 1; /* true */ return R; } void Random_delete(Random R) { free(R); } /* Returns the next random number in the sequence. */ double Random_nextDouble(Random R) { int k; int I = R->i; int J = R->j; int *m = R->m; k = m[I] - m[J]; if (k < 0) k += m1; R->m[J] = k; if (I == 0) I = 16; else I--; R->i = I; if (J == 0) J = 16 ; else J--; R->j = J; if (R->haveRange) return R->left + dm1 * (double) k * R->width; else return dm1 * (double) k; } /*-------------------------------------------------------------------- PRIVATE METHODS ----------------------------------------------------------------- */ static void initialize(Random R, int seed) { int jseed, k0, k1, j0, j1, iloop; dm1 = 1.0 / (double) m1; R->seed = seed; if (seed < 0 ) seed = -seed; /* seed = abs(seed) */ jseed = (seed < m1 ? seed : m1); /* jseed = min(seed, m1) */ if (jseed % 2 == 0) --jseed; k0 = 9069 % m2; k1 = 9069 / m2; j0 = jseed % m2; j1 = jseed / m2; for (iloop = 0; iloop < 17; ++iloop) { jseed = j0 * k0; j1 = (jseed / m2 + j0 * k1 + j1 * k0) % (m2 / 2); j0 = jseed % m2; R->m[iloop] = j0 + m2 * j1; } R->i = 4; R->j = 16; } double *RandomVector(int N, Random R) { int i; double *x = (double *) malloc(sizeof(double)*N); for (i=0; i