Random.c

#include <malloc.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;
}