#include "randlib.h"
#include <stdio.h>
#include <stdlib.h>
static long Xm1,Xm2,Xa1,Xa2,Xcg1[32],Xcg2[32],Xa1w,Xa2w,Xig1[32],Xig2[32],
Xlg1[32],Xlg2[32],Xa1vw,Xa2vw;
static long Xqanti[32];
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(0);
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
}