Performance/bin/Random.c

0001
0002
0003 #include <stdlib.h>
0004
0005 #include "Random.h"
0006
0007 #ifndef NULL
0008 #define NULL 0
0009 #endif
0010
0011 /* static const int mdig = 32; */
0012 #define MDIG 32
0013
0014 /* static const int one = 1; */
0015 #define ONE 1
0016
0017 static const int m1 = (ONE << (MDIG - 2)) + ((ONE << (MDIG - 2)) - ONE);
0018 static const int m2 = ONE << MDIG / 2;
0019
0020 /* For mdig = 32 : m1 =          2147483647, m2 =      65536
0021      For mdig = 64 : m1 = 9223372036854775807, m2 = 4294967296
0022   */
0023
0024 /* move to initialize() because  */
0025 /* compiler could not resolve as */
0026 /*   a constant.                 */
0027
0028 static /*const*/ double dm1; /*  = 1.0 / (double) m1; */
0029
0030 /* private methods (defined below, but not in Random.h */
0031
0032 static void initialize(Random R, int seed);
0033
0034 Random new_Random_seed(int seed) {
0035   Random R = (Random)malloc(sizeof(Random_struct));
0036
0037   initialize(R, seed);
0038   R->left = 0.0;
0039   R->right = 1.0;
0040   R->width = 1.0;
0041   R->haveRange = 0 /*false*/;
0042
0043   return R;
0044 }
0045
0046 Random new_Random(int seed, double left, double right) {
0047   Random R = (Random)malloc(sizeof(Random_struct));
0048
0049   initialize(R, seed);
0050   R->left = left;
0051   R->right = right;
0052   R->width = right - left;
0053   R->haveRange = 1; /* true */
0054
0055   return R;
0056 }
0057
0058 void Random_delete(Random R) { free(R); }
0059
0060 /* Returns the next random number in the sequence.  */
0061
0062 double Random_nextDouble(Random R) {
0063   int k;
0064
0065   int I = R->i;
0066   int J = R->j;
0067   int *m = R->m;
0068
0069   k = m[I] - m[J];
0070   if (k < 0)
0071     k += m1;
0072   R->m[J] = k;
0073
0074   if (I == 0)
0075     I = 16;
0076   else
0077     I--;
0078   R->i = I;
0079
0080   if (J == 0)
0081     J = 16;
0082   else
0083     J--;
0084   R->j = J;
0085
0086   if (R->haveRange)
0087     return R->left + dm1 * (double)k * R->width;
0088   else
0089     return dm1 * (double)k;
0090 }
0091
0092 /*--------------------------------------------------------------------
0093                            PRIVATE METHODS
0094   ----------------------------------------------------------------- */
0095
0096 static void initialize(Random R, int seed) {
0097   int jseed, k0, k1, j0, j1, iloop;
0098
0099   dm1 = 1.0 / (double)m1;
0100
0101   R->seed = seed;
0102
0103   if (seed < 0)
0104     seed = -seed;                  /* seed = abs(seed) */
0105   jseed = (seed < m1 ? seed : m1); /* jseed = min(seed, m1) */
0106   if (jseed % 2 == 0)
0107     --jseed;
0108   k0 = 9069 % m2;
0109   k1 = 9069 / m2;
0110   j0 = jseed % m2;
0111   j1 = jseed / m2;
0112   for (iloop = 0; iloop < 17; ++iloop) {
0113     jseed = j0 * k0;
0114     j1 = (jseed / m2 + j0 * k1 + j1 * k0) % (m2 / 2);
0115     j0 = jseed % m2;
0116     R->m[iloop] = j0 + m2 * j1;
0117   }
0118   R->i = 4;
0119   R->j = 16;
0120 }
0121
0122 double *RandomVector(int N, Random R) {
0123   int i;
0124   double *x = (double *)malloc(sizeof(double) * N);
0125
0126   for (i = 0; i < N; i++)
0127     x[i] = Random_nextDouble(R);
0128
0129   return x;
0130 }
0131
0132 double **RandomMatrix(int M, int N, Random R) {
0133   int i;
0134   int j;
0135
0136   /* allocate matrix */
0137
0138   double **A = (double **)malloc(sizeof(double *) * M);
0139
0140   if (A == NULL)
0141     return NULL;
0142
0143   for (i = 0; i < M; i++) {
0144     A[i] = (double *)malloc(sizeof(double) * N);
0145     if (A[i] == NULL) {
0146       free(A);
0147       return NULL;
0148     }
0149     for (j = 0; j < N; j++)
0150       A[i][j] = Random_nextDouble(R);
0151   }
0152   return A;
0153 }