modules/Math-Random/com.c

   1 #include "randlib.h"
   2 #include <stdio.h>
   3 #include <stdlib.h>
   4 static long Xm1,Xm2,Xa1,Xa2,Xcg1[32],Xcg2[32],Xa1w,Xa2w,Xig1[32],Xig2[32],
   5 Xlg1[32],Xlg2[32],Xa1vw,Xa2vw;
   6 static long Xqanti[32];
   7 void advnst(long k)
   8 /*
   9 **********************************************************************
  10      void advnst(long k)
  11                ADV-a-N-ce ST-ate
  12      Advances the state  of  the current  generator  by 2^K values  and
  13      resets the initial seed to that value.
  14      This is  a  transcription from   Pascal to  Fortran    of  routine
  15      Advance_State from the paper
  16      L'Ecuyer, P. and  Cote, S. "Implementing  a  Random Number Package
  17      with  Splitting   Facilities."  ACM  Transactions  on Mathematical
  18      Software, 17:98-111 (1991)
  19                               Arguments
  20      k -> The generator is advanced by2^K values
  21 **********************************************************************
  22 */
  23 {
  24 #define numg 32L
  25 extern void gsrgs(long getset,long *qvalue);
  26 extern void gscgn(long getset,long *g);
  27 extern long Xm1,Xm2,Xa1,Xa2,Xcg1[],Xcg2[];
  28 static long g,i,ib1,ib2;
  29 static long qrgnin;
  30 /*
  31      Abort unless random number generator initialized
  32 */
  33     gsrgs(0L,&qrgnin);
  34     if(qrgnin) goto S10;
  35     fputs(" ADVNST called before random generator initialized - ABORT\n",
  36           stderr);
  37     exit(1);
  38 S10:
  39     gscgn(0L,&g);
  40     ib1 = Xa1;
  41     ib2 = Xa2;
  42     for(i=1; i<=k; i++) {
  43         ib1 = mltmod(ib1,ib1,Xm1);
  44         ib2 = mltmod(ib2,ib2,Xm2);
  45     }
  46     setsd(mltmod(ib1,*(Xcg1+g-1),Xm1),mltmod(ib2,*(Xcg2+g-1),Xm2));
  47 /*
  48      NOW, IB1 = A1**K AND IB2 = A2**K
  49 */
  50 #undef numg
  51 }
  52 void getsd(long *iseed1,long *iseed2)
  53 /*
  54 **********************************************************************
  55      void getsd(long *iseed1,long *iseed2)
  56                GET SeeD
  57      Returns the value of two integer seeds of the current generator
  58      This  is   a  transcription from  Pascal   to  Fortran  of routine
  59      Get_State from the paper
  60      L'Ecuyer, P. and  Cote,  S. "Implementing a Random Number  Package
  61      with   Splitting Facilities."  ACM  Transactions   on Mathematical
  62      Software, 17:98-111 (1991)
  63                               Arguments
  64      iseed1 <- First integer seed of generator G
  65      iseed2 <- Second integer seed of generator G
  66 **********************************************************************
  67 */
  68 {
  69 #define numg 32L
  70 extern void gsrgs(long getset,long *qvalue);
  71 extern void gscgn(long getset,long *g);
  72 extern long Xcg1[],Xcg2[];
  73 static long g;
  74 static long qrgnin;
  75 /*
  76      Abort unless random number generator initialized
  77 */
  78     gsrgs(0L,&qrgnin);
  79     if(qrgnin) goto S10;
  80     fprintf(stderr,"%s\n",
  81       " GETSD called before random number generator  initialized -- abort!");
  82     exit(0);
  83 S10:
  84     gscgn(0L,&g);
  85     *iseed1 = *(Xcg1+g-1);
  86     *iseed2 = *(Xcg2+g-1);
  87 #undef numg
  88 }
  89 long ignlgi(void)
  90 /*
  91 **********************************************************************
  92      long ignlgi(void)
  93                GeNerate LarGe Integer
  94      Returns a random integer following a uniform distribution over
  95      (1, 2147483562) using the current generator.
  96      This is a transcription from Pascal to Fortran of routine
  97      Random from the paper
  98      L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
  99      with Splitting Facilities." ACM Transactions on Mathematical
 100      Software, 17:98-111 (1991)
 101 **********************************************************************
 102 */
 103 {
 104 #define numg 32L
 105 extern void gsrgs(long getset,long *qvalue);
 106 extern void gssst(long getset,long *qset);
 107 extern void gscgn(long getset,long *g);
 108 extern void inrgcm(void);
 109 extern long Xm1,Xm2,Xa1,Xa2,Xcg1[],Xcg2[];
 110 extern long Xqanti[];
 111 static long ignlgi,curntg,k,s1,s2,z;
 112 static long qqssd,qrgnin;
 113 /*
 114      IF THE RANDOM NUMBER PACKAGE HAS NOT BEEN INITIALIZED YET, DO SO.
 115      IT CAN BE INITIALIZED IN ONE OF TWO WAYS : 1) THE FIRST CALL TO
 116      THIS ROUTINE  2) A CALL TO SETALL.
 117 */
 118     gsrgs(0L,&qrgnin);
 119     if(!qrgnin) inrgcm();
 120     gssst(0,&qqssd);
 121     if(!qqssd) setall(1234567890L,123456789L);
 122 /*
 123      Get Current Generator
 124 */
 125     gscgn(0L,&curntg);
 126     s1 = *(Xcg1+curntg-1);
 127     s2 = *(Xcg2+curntg-1);
 128     k = s1/53668L;
 129     s1 = Xa1*(s1-k*53668L)-k*12211;
 130     if(s1 < 0) s1 += Xm1;
 131     k = s2/52774L;
 132     s2 = Xa2*(s2-k*52774L)-k*3791;
 133     if(s2 < 0) s2 += Xm2;
 134     *(Xcg1+curntg-1) = s1;
 135     *(Xcg2+curntg-1) = s2;
 136     z = s1-s2;
 137     if(z < 1) z += (Xm1-1);
 138     if(*(Xqanti+curntg-1)) z = Xm1-z;
 139     ignlgi = z;
 140     return ignlgi;
 141 #undef numg
 142 }
 143 void initgn(long isdtyp)
 144 /*
 145 **********************************************************************
 146      void initgn(long isdtyp)
 147           INIT-ialize current G-e-N-erator
 148      Reinitializes the state of the current generator
 149      This is a transcription from Pascal to C of routine
 150      Init_Generator from the paper
 151      L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
 152      with Splitting Facilities." ACM Transactions on Mathematical
 153      Software, 17:98-111 (1991)
 154                               Arguments
 155      isdtyp -> The state to which the generator is to be set
 156           isdtyp = -1  => sets the seeds to their initial value
 157           isdtyp =  0  => sets the seeds to the first value of
 158                           the current block
 159           isdtyp =  1  => sets the seeds to the first value of
 160                           the next block
 161      WGR, 12/19/00: replaced S10, S20, etc. with C blocks {} per
 162      original paper.
 163 **********************************************************************
 164 */
 165 {
 166 #define numg 32L
 167 extern void gsrgs(long getset,long *qvalue);
 168 extern void gscgn(long getset,long *g);
 169 extern long Xm1,Xm2,Xa1w,Xa2w,Xig1[],Xig2[],Xlg1[],Xlg2[],Xcg1[],Xcg2[];
 170 static long g;
 171 static long qrgnin;
 172 /*
 173      Abort unless random number generator initialized
 174 */
 175     gsrgs(0L,&qrgnin);
 176  if (! qrgnin) {
 177     fprintf(stderr,"%s\n",
 178       " INITGN called before random number generator  initialized -- abort!");
 179     exit(1);
 180  }
 181
 182  gscgn(0L,&g);
 183  if(isdtyp == -1) {             /* Initial seed */
 184     *(Xlg1+g-1) = *(Xig1+g-1);
 185     *(Xlg2+g-1) = *(Xig2+g-1);
 186  }
 187  else if (isdtyp == 0) { ; }    /* Last seed */
 188  else if (isdtyp == 1) {        /* New seed */
 189     *(Xlg1+g-1) = mltmod(Xa1w,*(Xlg1+g-1),Xm1);
 190     *(Xlg2+g-1) = mltmod(Xa2w,*(Xlg2+g-1),Xm2);
 191  } else {
 192     fprintf(stderr,"%s\n","isdtyp not in range in INITGN");
 193     exit(1);
 194  }
 195
 196  *(Xcg1+g-1) = *(Xlg1+g-1);
 197  *(Xcg2+g-1) = *(Xlg2+g-1);
 198 #undef numg
 199 }
 200
 201 void inrgcm(void)
 202 /*
 203 **********************************************************************
 204      void inrgcm(void)
 205           INitialize Random number Generator CoMmon
 206                               Function
 207      Initializes common area  for random number  generator.  This saves
 208      the  nuisance  of  a  BLOCK DATA  routine  and the  difficulty  of
 209      assuring that the routine is loaded with the other routines.
 210 **********************************************************************
 211 */
 212 {
 213 #define numg 32L
 214 extern void gsrgs(long getset,long *qvalue);
 215 extern long Xm1,Xm2,Xa1,Xa2,Xa1w,Xa2w,Xa1vw,Xa2vw;
 216 extern long Xqanti[];
 217 static long T1;
 218 static long i;
 219 /*
 220      V=20;                            W=30;
 221      A1W = MOD(A1**(2**W),M1)         A2W = MOD(A2**(2**W),M2)
 222      A1VW = MOD(A1**(2**(V+W)),M1)    A2VW = MOD(A2**(2**(V+W)),M2)
 223    If V or W is changed A1W, A2W, A1VW, and A2VW need to be recomputed.
 224     An efficient way to precompute a**(2*j) MOD m is to start with
 225     a and square it j times modulo m using the function MLTMOD.
 226 */
 227     Xm1 = 2147483563L;
 228     Xm2 = 2147483399L;
 229     Xa1 = 40014L;
 230     Xa2 = 40692L;
 231     Xa1w = 1033780774L;
 232     Xa2w = 1494757890L;
 233     Xa1vw = 2082007225L;
 234     Xa2vw = 784306273L;
 235     for(i=0; i<numg; i++) *(Xqanti+i) = 0;
 236     T1 = 1;
 237 /*
 238      Tell the world that common has been initialized
 239 */
 240     gsrgs(1L,&T1);
 241 #undef numg
 242 }
 243 void setall(long iseed1,long iseed2)
 244 /*
 245 **********************************************************************
 246      void setall(long iseed1,long iseed2)
 247                SET ALL random number generators
 248      Sets the initial seed of generator 1 to ISEED1 and ISEED2. The
 249      initial seeds of the other generators are set accordingly, and
 250      all generators states are set to these seeds.
 251      This is a transcription from Pascal to Fortran of routine
 252      Set_Initial_Seed from the paper
 253      L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
 254      with Splitting Facilities." ACM Transactions on Mathematical
 255      Software, 17:98-111 (1991)
 256                               Arguments
 257      iseed1 -> First of two integer seeds
 258      iseed2 -> Second of two integer seeds
 259 **********************************************************************
 260 */
 261 {
 262 #define numg 32L
 263 extern void gsrgs(long getset,long *qvalue);
 264 extern void gssst(long getset,long *qset);
 265 extern void gscgn(long getset,long *g);
 266 extern long Xm1,Xm2,Xa1vw,Xa2vw,Xig1[],Xig2[];
 267 static long T1;
 268 static long g,ocgn;
 269 static long qrgnin;
 270     T1 = 1;
 271 /*
 272      TELL IGNLGI, THE ACTUAL NUMBER GENERATOR, THAT THIS ROUTINE
 273       HAS BEEN CALLED.
 274 */
 275     gssst(1,&T1);
 276     gscgn(0L,&ocgn);
 277 /*
 278      Initialize Common Block if Necessary
 279 */
 280     gsrgs(0L,&qrgnin);
 281     if(!qrgnin) inrgcm();
 282     *Xig1 = iseed1;
 283     *Xig2 = iseed2;
 284     initgn(-1L);
 285     for(g=2; g<=numg; g++) {
 286         *(Xig1+g-1) = mltmod(Xa1vw,*(Xig1+g-2),Xm1);
 287         *(Xig2+g-1) = mltmod(Xa2vw,*(Xig2+g-2),Xm2);
 288         gscgn(1L,&g);
 289         initgn(-1L);
 290     }
 291     gscgn(1L,&ocgn);
 292 #undef numg
 293 }
 294 void setant(long qvalue)
 295 /*
 296 **********************************************************************
 297      void setant(long qvalue)
 298                SET ANTithetic
 299      Sets whether the current generator produces antithetic values.  If
 300      X   is  the value  normally returned  from  a uniform [0,1] random
 301      number generator then 1  - X is the antithetic  value. If X is the
 302      value  normally  returned  from a   uniform  [0,N]  random  number
 303      generator then N - 1 - X is the antithetic value.
 304      All generators are initialized to NOT generate antithetic values.
 305      This is a transcription from Pascal to Fortran of routine
 306      Set_Antithetic from the paper
 307      L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
 308      with Splitting Facilities." ACM Transactions on Mathematical
 309      Software, 17:98-111 (1991)
 310                               Arguments
 311      qvalue -> nonzero if generator G is to generating antithetic
 312                     values, otherwise zero
 313 **********************************************************************
 314 */
 315 {
 316 #define numg 32L
 317 extern void gsrgs(long getset,long *qvalue);
 318 extern void gscgn(long getset,long *g);
 319 extern long Xqanti[];
 320 static long g;
 321 static long qrgnin;
 322 /*
 323      Abort unless random number generator initialized
 324 */
 325     gsrgs(0L,&qrgnin);
 326     if(qrgnin) goto S10;
 327     fprintf(stderr,"%s\n",
 328       " SETANT called before random number generator  initialized -- abort!");
 329     exit(1);
 330 S10:
 331     gscgn(0L,&g);
 332     Xqanti[g-1] = qvalue;
 333 #undef numg
 334 }
 335 void setsd(long iseed1,long iseed2)
 336 /*
 337 **********************************************************************
 338      void setsd(long iseed1,long iseed2)
 339                SET S-ee-D of current generator
 340      Resets the initial  seed of  the current  generator to  ISEED1 and
 341      ISEED2. The seeds of the other generators remain unchanged.
 342      This is a transcription from Pascal to Fortran of routine
 343      Set_Seed from the paper
 344      L'Ecuyer, P. and Cote, S. "Implementing a Random Number Package
 345      with Splitting Facilities." ACM Transactions on Mathematical
 346      Software, 17:98-111 (1991)
 347                               Arguments
 348      iseed1 -> First integer seed
 349      iseed2 -> Second integer seed
 350 **********************************************************************
 351 */
 352 {
 353 #define numg 32L
 354 extern void gsrgs(long getset,long *qvalue);
 355 extern void gscgn(long getset,long *g);
 356 extern long Xig1[],Xig2[];
 357 static long g;
 358 static long qrgnin;
 359 /*
 360      Abort unless random number generator initialized
 361 */
 362     gsrgs(0L,&qrgnin);
 363     if(qrgnin) goto S10;
 364     fprintf(stderr,"%s\n",
 365       " SETSD called before random number generator  initialized -- abort!");
 366     exit(1);
 367 S10:
 368     gscgn(0L,&g);
 369     *(Xig1+g-1) = iseed1;
 370     *(Xig2+g-1) = iseed2;
 371     initgn(-1L);
 372 #undef numg
 373 }