src/jnt/scimark2/Random.java

package jnt.scimark2;

/* Random.java based on Java Numerical Toolkit (JNT) Random.UniformSequence
      class.  We do not use Java's own java.util.Random so that we can compare
      results with equivalent C and Fortran coces.
*/

public class Random {


/* ------------------------------------------------------------------------------
                                 CLASS VARIABLES
     ------------------------------------------------------------------------------ */

    int seed = 0;

    private int m[];
    private int i = 4;
    private int j = 16;

    private final int mdig = 32;
    private final int one = 1;
    private final int m1 = (one << mdig-2) + ((one << mdig-2)-one);
    private final int m2 = one << mdig/2;

    /* For mdig = 32 : m1 =          2147483647, m2 =      65536
       For mdig = 64 : m1 = 9223372036854775807, m2 = 4294967296 
    */

    private double dm1 = 1.0 / (double) m1;

    private boolean haveRange = false;
    private double left  = 0.0;
    private double right = 1.0;
    private double width = 1.0;


/* ------------------------------------------------------------------------------
                                  CONSTRUCTORS
     ------------------------------------------------------------------------------ */

/**
     Initializes a sequence of uniformly distributed quasi random numbers with a
     seed based on the system clock.
*/
    public Random () {
        initialize( (int) System.currentTimeMillis());
    }

/**
     Initializes a sequence of uniformly distributed quasi random numbers on a
     given half-open interval [left,right) with a seed based on the system
     clock.

@param <B>left</B> (double)<BR>

         The left endpoint of the half-open interval [left,right).

@param <B>right</B> (double)<BR>

         The right endpoint of the half-open interval [left,right).
*/
    public Random ( double left, double right) {
        initialize( (int) System.currentTimeMillis() );
        this.left = left;
        this.right = right;
        width = right - left;
        haveRange = true;
    }

/**
     Initializes a sequence of uniformly distributed quasi random numbers with a
     given seed.

@param <B>seed</B> (int)<BR>

         The seed of the random number generator.  Two sequences with the same
         seed will be identical.
*/
    public Random (int seed) {
        initialize( seed);
    }

/**
     Initializes a sequence of uniformly distributed quasi random numbers
     with a given seed on a given half-open interval [left,right).

@param <B>seed</B> (int)<BR>

         The seed of the random number generator.  Two sequences with the same
         seed will be identical.

@param <B>left</B> (double)<BR>

         The left endpoint of the half-open interval [left,right).

@param <B>right</B> (double)<BR>

         The right endpoint of the half-open interval [left,right).
*/
    public Random (int seed, double left, double right) {
        initialize( seed);
        this.left = left;
        this.right = right;
        width = right - left;
        haveRange = true;
    }

/* ------------------------------------------------------------------------------
                               PUBLIC METHODS
     ------------------------------------------------------------------------------ */

/**
     Returns the next random number in the sequence.
*/
    public final synchronized double nextDouble () {

        int k;
        double nextValue;

        k = m[i] - m[j];
        if (k < 0) k += m1;
            m[j] = k;

        if (i == 0) 
            i = 16;
        else i--;

        if (j == 0) 
            j = 16 ;
        else j--;

        if (haveRange) 
            return  left +  dm1 * (double) k * width;
        else
            return dm1 * (double) k;

    } 

/**
     Returns the next N random numbers in the sequence, as
     a vector.
*/
    public final synchronized void nextDoubles (double x[]) {

        int N = x.length;
        int remainder = N & 3;        // N mod 4

        if (haveRange) {
            for (int count=0; count<N; count++) {
                int k = m[i] - m[j];

                if (i == 0) i = 16;
                else i--;

                if (k < 0) k += m1;
                m[j] = k;

                if (j == 0) j = 16;
                else j--;

                x[count] = left + dm1 * (double) k * width;
            }
        }
        else
        {

            for (int count=0; count<remainder; count++) {
                int k = m[i] - m[j];

                if (i == 0) i = 16;
                else i--;

                if (k < 0) k += m1;
                m[j] = k;

                if (j == 0) j = 16;
                else j--;


                x[count] = dm1 * (double) k;
            }

            for (int count=remainder; count<N; count+=4) {
                int k = m[i] - m[j];
                if (i == 0) i = 16;
                else i--;
                if (k < 0) k += m1;
                m[j] = k;
                if (j == 0) j = 16;
                else j--;
                x[count] = dm1 * (double) k;


                k = m[i] - m[j];
                if (i == 0) i = 16;
                else i--;
                if (k < 0) k += m1;
                m[j] = k;
                if (j == 0) j = 16;
                else j--;
                x[count+1] = dm1 * (double) k;


                k = m[i] - m[j];
                if (i == 0) i = 16;
                else i--;
                if (k < 0) k += m1;
                m[j] = k;
                if (j == 0) j = 16;
                else j--;
                x[count+2] = dm1 * (double) k;


                k = m[i] - m[j];
                if (i == 0) i = 16;
                else i--;
                if (k < 0) k += m1;
                m[j] = k;
                if (j == 0) j = 16;
                else j--;
                x[count+3] = dm1 * (double) k;
            }
        }
    }


/*----------------------------------------------------------------------------
                             PRIVATE METHODS
    ------------------------------------------------------------------------ */

    private void initialize (int seed) {

        int jseed, k0, k1, j0, j1, iloop;

        this.seed = seed;

        m = new int[17];

        jseed = Math.min(Math.abs(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;
            m[iloop] = j0 + m2 * j1;
        }
        i = 4;
        j = 16;

    }

}