| |
| |
| #include <stdlib.h> |
| #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<N; i++) |
| x[i] = Random_nextDouble(R); |
| |
| return x; |
| } |
| |
| |
| double **RandomMatrix(int M, int N, Random R) |
| { |
| int i; |
| int j; |
| |
| /* allocate matrix */ |
| |
| double **A = (double **) malloc(sizeof(double*)*M); |
| |
| if (A == NULL) return NULL; |
| |
| for (i=0; i<M; i++) |
| { |
| A[i] = (double *) malloc(sizeof(double)*N); |
| if (A[i] == NULL) |
| { |
| free(A); |
| return NULL; |
| } |
| for (j=0; j<N; j++) |
| A[i][j] = Random_nextDouble(R); |
| } |
| return A; |
| } |
| |
| |
| |