Random.c

#include <stdlib.h>

#include "Random.h"

/* 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;

    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;

    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)
{
    double *x = (double *) malloc(sizeof(double) * N);

    for (int i = 0; i < N; i++)
        x[i] = Random_nextDouble(R);

    return x;
}

double **RandomMatrix(int M, int N, Random R)
{
    /* allocate matrix */
    double **A = (double **) malloc(sizeof(double *) * M);

    if (A == NULL)
        return NULL;

    for (int i = 0; i < M; i++) {
        A[i] = (double *) malloc(sizeof(double) * N);
        if (A[i] == NULL) {
            free(A);
            return NULL;
        }
        for (int j = 0; j < N; j++)
            A[i][j] = Random_nextDouble(R);
    }
    return A;
}