com.c

#include "randlib.h"
#include <stdio.h>
#include <stdlib.h>

long Xm1, Xm2, Xa1, Xa2, Xcg1[32], Xcg2[32], Xqanti[32];
long Xa1w, Xa2w, Xig1[32], Xig2[32], Xlg1[32], Xlg2[32], Xa1vw, Xa2vw;

void advnst(long k)
/*
**********************************************************************
     void advnst(long k)
               ADV-a-N-ce ST-ate
     Advances the state  of  the current  generator  by 2^K values  and
     resets the initial seed to that value.
     This is  a  transcription from   Pascal to  Fortran    of  routine
     Advance_State from the paper
     L'Ecuyer, P. and  Cote, S. "Implementing  a  Random Number Package
     with  Splitting   Facilities."  ACM  Transactions  on Mathematical
     Software, 17:98-111 (1991)
                              Arguments
     k -> The generator is advanced by2^K values
**********************************************************************
*/
{
#define numg 32L
    extern void gsrgs(long getset, long *qvalue);
    extern void gscgn(long getset, long *g);
    extern long Xm1, Xm2, Xa1, Xa2, Xcg1[], Xcg2[];
    static long g, i, ib1, ib2;
    static long qrgnin;
/*
     Abort unless random number generator initialized
*/
    gsrgs(0L, &qrgnin);
    if (qrgnin)
        goto S10;
    fputs(" ADVNST called before random generator initialized - ABORT\n", stderr);
    exit(1);
  S10:
    gscgn(0L, &g);
    ib1 = Xa1;
    ib2 = Xa2;
    for (i = 1; i <= k; i++) {
        ib1 = mltmod(ib1, ib1, Xm1);
        ib2 = mltmod(ib2, ib2, Xm2);
    }
    setsd(mltmod(ib1, *(Xcg1 + g - 1), Xm1), mltmod(ib2, *(Xcg2 + g - 1), Xm2));
/*
     NOW, IB1 = A1**K AND IB2 = A2**K
*/
#undef numg
}

void getsd(long *iseed1, long *iseed2)
/*
**********************************************************************
     void getsd(long *iseed1,long *iseed2)
               GET SeeD
     Returns the value of two integer seeds of the current generator
     This  is   a  transcription from  Pascal   to  Fortran  of routine
     Get_State from the paper
     L'Ecuyer, P. and  Cote,  S. "Implementing a Random Number  Package
     with   Splitting Facilities."  ACM  Transactions   on Mathematical
     Software, 17:98-111 (1991)
                              Arguments
     iseed1 <- First integer seed of generator G
     iseed2 <- Second integer seed of generator G
**********************************************************************
*/
{
#define numg 32L
    extern void gsrgs(long getset, long *qvalue);
    extern void gscgn(long getset, long *g);
    extern long Xcg1[], Xcg2[];
    static long g;
    static long qrgnin;
/*
     Abort unless random number generator initialized
*/
    gsrgs(0L, &qrgnin);
    if (qrgnin)
        goto S10;
    fprintf(stderr, "%s\n",
            " GETSD called before random number generator  initialized -- abort!");
    exit(1);
  S10:
    gscgn(0L, &g);
    *iseed1 = *(Xcg1 + g - 1);
    *iseed2 = *(Xcg2 + g - 1);
#undef numg
}

long ignlgi(void)
/*
**********************************************************************
     long ignlgi(void)
               GeNerate LarGe Integer
     Returns a random integer following a uniform distribution over
     (1, 2147483562) using the current generator.
     This is a transcription from Pascal to Fortran of routine
     Random from the paper
     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
     with Splitting Facilities." ACM Transactions on Mathematical
     Software, 17:98-111 (1991)
**********************************************************************
*/
{
#define numg 32L
    extern void gsrgs(long getset, long *qvalue);
    extern void gssst(long getset, long *qset);
    extern void gscgn(long getset, long *g);
    extern void inrgcm(void);
    extern long Xm1, Xm2, Xa1, Xa2, Xcg1[], Xcg2[];
    extern long Xqanti[];
    static long ignlgi, curntg, k, s1, s2, z;
    static long qqssd, qrgnin;
/*
     IF THE RANDOM NUMBER PACKAGE HAS NOT BEEN INITIALIZED YET, DO SO.
     IT CAN BE INITIALIZED IN ONE OF TWO WAYS : 1) THE FIRST CALL TO
     THIS ROUTINE  2) A CALL TO SETALL.
*/
    gsrgs(0L, &qrgnin);
    if (!qrgnin)
        inrgcm();
    gssst(0, &qqssd);
    if (!qqssd)
        setall(1234567890L, 123456789L);
/*
     Get Current Generator
*/
    gscgn(0L, &curntg);
    s1 = *(Xcg1 + curntg - 1);
    s2 = *(Xcg2 + curntg - 1);
    k = s1 / 53668L;
    s1 = Xa1 * (s1 - k * 53668L) - k * 12211;
    if (s1 < 0)
        s1 += Xm1;
    k = s2 / 52774L;
    s2 = Xa2 * (s2 - k * 52774L) - k * 3791;
    if (s2 < 0)
        s2 += Xm2;
    *(Xcg1 + curntg - 1) = s1;
    *(Xcg2 + curntg - 1) = s2;
    z = s1 - s2;
    if (z < 1)
        z += (Xm1 - 1);
    if (*(Xqanti + curntg - 1))
        z = Xm1 - z;
    ignlgi = z;
    return ignlgi;
#undef numg
}

void initgn(long isdtyp)
/*
**********************************************************************
     void initgn(long isdtyp)
          INIT-ialize current G-e-N-erator
     Reinitializes the state of the current generator
     This is a transcription from Pascal to Fortran of routine
     Init_Generator from the paper
     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
     with Splitting Facilities." ACM Transactions on Mathematical
     Software, 17:98-111 (1991)
                              Arguments
     isdtyp -> The state to which the generator is to be set
          isdtyp = -1  => sets the seeds to their initial value
          isdtyp =  0  => sets the seeds to the first value of
                          the current block
          isdtyp =  1  => sets the seeds to the first value of
                          the next block
**********************************************************************
*/
{
#define numg 32L
    extern void gsrgs(long getset, long *qvalue);
    extern void gscgn(long getset, long *g);
    extern long Xm1, Xm2, Xa1w, Xa2w, Xig1[], Xig2[], Xlg1[], Xlg2[], Xcg1[], Xcg2[];
    static long g;
    static long qrgnin;
/*
     Abort unless random number generator initialized
*/
    gsrgs(0L, &qrgnin);
    if (qrgnin)
        goto S10;
    fprintf(stderr, "%s\n",
            " INITGN called before random number generator  initialized -- abort!");
    exit(1);
  S10:
    gscgn(0L, &g);
    if (-1 != isdtyp)
        goto S20;
    *(Xlg1 + g - 1) = *(Xig1 + g - 1);
    *(Xlg2 + g - 1) = *(Xig2 + g - 1);
    goto S50;
  S20:
    if (0 != isdtyp)
        goto S30;
    goto S50;
  S30:
/*
     do nothing
*/
    if (1 != isdtyp)
        goto S40;
    *(Xlg1 + g - 1) = mltmod(Xa1w, *(Xlg1 + g - 1), Xm1);
    *(Xlg2 + g - 1) = mltmod(Xa2w, *(Xlg2 + g - 1), Xm2);
    goto S50;
  S40:
    fprintf(stderr, "%s\n", "isdtyp not in range in INITGN");
    exit(1);
  S50:
    *(Xcg1 + g - 1) = *(Xlg1 + g - 1);
    *(Xcg2 + g - 1) = *(Xlg2 + g - 1);
#undef numg
}

void inrgcm(void)
/*
**********************************************************************
     void inrgcm(void)
          INitialize Random number Generator CoMmon
                              Function
     Initializes common area  for random number  generator.  This saves
     the  nuisance  of  a  BLOCK DATA  routine  and the  difficulty  of
     assuring that the routine is loaded with the other routines.
**********************************************************************
*/
{
#define numg 32L
    extern void gsrgs(long getset, long *qvalue);
    extern long Xm1, Xm2, Xa1, Xa2, Xa1w, Xa2w, Xa1vw, Xa2vw;
    extern long Xqanti[];
    static long T1;
    static long i;
/*
     V=20;                            W=30;
     A1W = MOD(A1**(2**W),M1)         A2W = MOD(A2**(2**W),M2)
     A1VW = MOD(A1**(2**(V+W)),M1)    A2VW = MOD(A2**(2**(V+W)),M2)
   If V or W is changed A1W, A2W, A1VW, and A2VW need to be recomputed.
    An efficient way to precompute a**(2*j) MOD m is to start with
    a and square it j times modulo m using the function MLTMOD.
*/
    Xm1 = 2147483563L;
    Xm2 = 2147483399L;
    Xa1 = 40014L;
    Xa2 = 40692L;
    Xa1w = 1033780774L;
    Xa2w = 1494757890L;
    Xa1vw = 2082007225L;
    Xa2vw = 784306273L;
    for (i = 0; i < numg; i++)
        *(Xqanti + i) = 0;
    T1 = 1;
/*
     Tell the world that common has been initialized
*/
    gsrgs(1L, &T1);
#undef numg
}

void setall(long iseed1, long iseed2)
/*
**********************************************************************
     void setall(long iseed1,long iseed2)
               SET ALL random number generators
     Sets the initial seed of generator 1 to ISEED1 and ISEED2. The
     initial seeds of the other generators are set accordingly, and
     all generators states are set to these seeds.
     This is a transcription from Pascal to Fortran of routine
     Set_Initial_Seed from the paper
     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
     with Splitting Facilities." ACM Transactions on Mathematical
     Software, 17:98-111 (1991)
                              Arguments
     iseed1 -> First of two integer seeds
     iseed2 -> Second of two integer seeds
**********************************************************************
*/
{
#define numg 32L
    extern void gsrgs(long getset, long *qvalue);
    extern void gssst(long getset, long *qset);
    extern void gscgn(long getset, long *g);
    extern long Xm1, Xm2, Xa1vw, Xa2vw, Xig1[], Xig2[];
    static long T1;
    static long g, ocgn;
    static long qrgnin;
    T1 = 1;
/*
     TELL IGNLGI, THE ACTUAL NUMBER GENERATOR, THAT THIS ROUTINE
      HAS BEEN CALLED.
*/
    gssst(1, &T1);
    gscgn(0L, &ocgn);
/*
     Initialize Common Block if Necessary
*/
    gsrgs(0L, &qrgnin);
    if (!qrgnin)
        inrgcm();
    *Xig1 = iseed1;
    *Xig2 = iseed2;
    initgn(-1L);
    for (g = 2; g <= numg; g++) {
        *(Xig1 + g - 1) = mltmod(Xa1vw, *(Xig1 + g - 2), Xm1);
        *(Xig2 + g - 1) = mltmod(Xa2vw, *(Xig2 + g - 2), Xm2);
        gscgn(1L, &g);
        initgn(-1L);
    }
    gscgn(1L, &ocgn);
#undef numg
}

void setant(long qvalue)
/*
**********************************************************************
     void setant(long qvalue)
               SET ANTithetic
     Sets whether the current generator produces antithetic values.  If
     X   is  the value  normally returned  from  a uniform [0,1] random
     number generator then 1  - X is the antithetic  value. If X is the
     value  normally  returned  from a   uniform  [0,N]  random  number
     generator then N - 1 - X is the antithetic value.
     All generators are initialized to NOT generate antithetic values.
     This is a transcription from Pascal to Fortran of routine
     Set_Antithetic from the paper
     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
     with Splitting Facilities." ACM Transactions on Mathematical
     Software, 17:98-111 (1991)
                              Arguments
     qvalue -> nonzero if generator G is to generating antithetic
                    values, otherwise zero
**********************************************************************
*/
{
#define numg 32L
    extern void gsrgs(long getset, long *qvalue);
    extern void gscgn(long getset, long *g);
    extern long Xqanti[];
    static long g;
    static long qrgnin;
/*
     Abort unless random number generator initialized
*/
    gsrgs(0L, &qrgnin);
    if (qrgnin)
        goto S10;
    fprintf(stderr, "%s\n",
            " SETANT called before random number generator  initialized -- abort!");
    exit(1);
  S10:
    gscgn(0L, &g);
    Xqanti[g - 1] = qvalue;
#undef numg
}

void setsd(long iseed1, long iseed2)
/*
**********************************************************************
     void setsd(long iseed1,long iseed2)
               SET S-ee-D of current generator
     Resets the initial  seed of  the current  generator to  ISEED1 and
     ISEED2. The seeds of the other generators remain unchanged.
     This is a transcription from Pascal to Fortran of routine
     Set_Seed from the paper
     L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
     with Splitting Facilities." ACM Transactions on Mathematical
     Software, 17:98-111 (1991)
                              Arguments
     iseed1 -> First integer seed
     iseed2 -> Second integer seed
**********************************************************************
*/
{
#define numg 32L
    extern void gsrgs(long getset, long *qvalue);
    extern void gscgn(long getset, long *g);
    extern long Xig1[], Xig2[];
    static long g;
    static long qrgnin;
/*
     Abort unless random number generator initialized
*/
    gsrgs(0L, &qrgnin);
    if (qrgnin)
        goto S10;
    fprintf(stderr, "%s\n",
            " SETSD called before random number generator  initialized -- abort!");
    exit(1);
  S10:
    gscgn(0L, &g);
    *(Xig1 + g - 1) = iseed1;
    *(Xig2 + g - 1) = iseed2;
    initgn(-1L);
#undef numg
}