From:	CSBVAX::CSBVAX::MRGATE::"SMTP::CRVAX.SRI.COM::RELAY-INFO-VAX" 24-FEB-1989 14:52
To:	MRGATE::"ARISIA::EVERHART"
Subj:	Random Number Generators


Received: From KL.SRI.COM by CRVAX.SRI.COM with TCP; Fri, 24 FEB 89 10:31:28 PDT
Received: from mis.arizona.edu (misvax.mis.arizona.edu) by KL.SRI.COM with TCP; Fri, 24 Feb 89 10:16:03 PST
Received: by mis.arizona.edu id <20A004D9141@mis.arizona.edu> ;
       Fri, 24 Feb 89 11:12:20 MST
Date: Fri, 24 Feb 89 11:10:54 MST
From: Mail is *not* all it's cracked up to be. <JMS@mis.arizona.edu>
Subject: Random Number Generators
To: info-vax@kl.sri.com
X-VMS-Mail-To: EXOS%"info-vax@kl.sri.com"
Message-ID: <890224111054.20A004D9141@mis.arizona.edu>

Here is a routine which we use for Random Number Generation in a
simulation system:

jms

Joel M Snyder, Univ of Arizona Dep't of MIS, Tucson, AZ, 85721  
BITNET: jms@arizmis  Internet: jms@mis.arizona.edu Phone: 602.621.2748
"The 3 B's of Spring Break ... Beer, Beach, and Babes." - Miller Beer Ad
--------------------------------------------------------------

DOCUMENTATION OF SIMLIB RANDOM NUMBER GENERATOR ROUTINES.
=========================================================

There are 4 functions associated with generating uniform ( 0, 1 )
random numbers in SIMLIB.  They are :

function   arguments    argement purpose   function purpose
--------   ----------   ----------------   -------------------------
 RANUN     (r) DUMMY    nothing - dummy    generate a uniform from
                                           stream 1

 RANUM     (i) STREAM   stream desired,    generate a uniform from
                        between 1 and 10   desired stream

 GETSED    (i) STREAM   stream desired,    retrieve the next integer
                        between 1 and 10   which will be used to 
           (i) SEED     integer which      generate the next uniform
                        will be used as    from the desired stream
                        a seed for the
                        stream, output

 SETSED    (i) STREAM   stream desired,    set the next integer which
                        between 1 and 10   will be used to generate
           (i) SEED     integer which      the next uniform from the
                        will be used as    desired stream
                        a seed for the 
                        stream, input

Note : The mode of an argument is given as (i) : integer or (r) : real.
       GETSED and SETSED return positive values if the operation is 
       a success and negative values if the operation is a failure.
       RANUN and RANUM need no initialization.  Further, no action is
       necessary to continue the production of new random uniforms if
       a do-loop is set up around the simulation in your main program ;
       e.g., the following will work :

           DO 1000 SIMRUM = 1, NUMRUN

              CALL INITLK 
              ...
       100    CALL TIMING
              ...

              ELSE IF ( NEXT .EQ. IENDSM ) THEN
                 CALL REPORT
                 GO TO 1000
              ENDIF
              GO TO 100

      1000 CONTINUE

        
Example usage where a do-loop is NOT used around the simulation logic :

C  ---Attempt to retrieve the seed for stream 1 from a data file and then
C  ---initialize the random number generator.

      OPEN  ( UNIT = 3,       FILE = 'SEED.DAT', 
     -        STATUS = 'OLD', ERR  = 20          )

C     ---The file exists ;  the model must have been run before and
C     ---the seed written to a file.

         READ  ( 3, * ) ISEED
         CLOSE ( UNIT = 3 )

         IF ( SETSED ( 1, ISEED ) .LT. 0. ) THEN
            STOP 'ERROR IN SETSED : BAD SEED VALUE'
         ENDIF

C  ---If the file did not exist then there was nothing to do.

   20 CONTINUE

...

      ELSE IF ( NEXT .EQ. IENDSM ) THEN

C     ---End of simulation.  Write report and record seed for next
C     ---simulation run.

         CALL REPORT
         CLOSE ( UNIT = 2 )

         IF ( GETSEED ( 1, ISEED ) .LT. 0. ) THEN
            STOP 'ERROR IN GETSEED'
         ELSE
            OPEN  ( UNIT = 3, FILE = 'SEED.DAT', STATUS = 'NEW' )
            WRITE ( 3, * ) ISEED
            CLOSE ( UNIT = 3 )
         ENDIF

         STOP ' '


--------------------------------------------------

      FUNCTION RAND(ISTRM) 
  
*     Prime modulus multiplicative linear congruential generator 
*     Z(I) = (630360016 * Z(I - 1)) (MOD(2**31 - 1)), based on Marse 
*     and Roberts' portable random-number generator UNIRAN. 
*     Multiple (100) streams are supported, with seeds spaced 100,000 
*     apart.  Input argument ISTRM must be a positive integer constant
*     or variable giving the desired stream number. 
  
*     Usage: (Three options)  
*      
*            1. To obtain the next U(0,1) random number from stream  
*               ISTRM, execute:  
*                  U = RAND(ISTRM)  
*               The REAL variable U will contain the next random number.  
*      
*            2. To set the seed for stream ISTRM to a desired value 
*               IZIN, execute:  
*                  DUMMY = RANDST(IZIN,ISTRM)  
*               where IZIN must be an INTEGER constant or variable set 
*               to a number between 1 and 2147483646 (inclusive). 
*               Default seeds for all 100 streams are given in the code.  
*      
*            3. To obtain the current (most recently used) integer in  
*               the sequence being generated for stream ISTRM, execute:  
*                  DUMMY = RANDGT(IZOUT,ISTRM) 
*               Then the INTEGER variable IZOUT is returned with the  
*               current integer in the sequence being used for stream  
*               ISTRM. 
 
      INTEGER B2E15,B2E16,HI15,HI31,LOW15,LOWPRD,MODLUS, 
     1        MULT1,MULT2,OVFLOW,IZ(100) 
 
*     Force saving of IZ between calls. 
 
      SAVE IZ 
 
*     Define the constants. 
 
      DATA MULT1,MULT2/24112,26143/ 
      DATA B2E15,B2E16,MODLUS/32768,65536,2147483647/ 
  
*     Set the default seeds for all 100 streams. 
  
      DATA IZ/ 1973272912, 281629770,  20006270,1280689831,2096730329, 
     1         1933576050, 913566091, 246780520,1363774876, 604901985, 
     2         1511192140,1259851944, 824064364, 150493284, 242708531, 
     3           75253171,1964472944,1202299975, 233217322,1911216000, 
     4          726370533, 403498145, 993232223,1103205531, 762430696, 
     5         1922803170,1385516923,  76271663, 413682397, 726466604, 
     6          336157058,1432650381,1120463904, 595778810, 877722890, 
     7         1046574445,  68911991,2088367019, 748545416, 622401386, 
     8         2122378830, 640690903,1774806513,2132545692,2079249579, 
     9           78130110, 852776735,1187867272,1351423507,1645973084, 
     A         1997049139, 922510944,2045512870, 898585771, 243649545, 
     B         1004818771, 773686062, 403188473, 372279877,1901633463, 
     C          498067494,2087759558, 493157915, 597104727,1530940798, 
     D         1814496276, 536444882,1663153658, 855503735,  67784357, 
     E         1432404475, 619691088, 119025595, 880802310, 176192644, 
     F         1116780070, 277854671,1366580350,1142483975,2026948561, 
     G         1053920743, 786262391,1792203830,1494667770,1923011392, 
     H         1433700034,1244184613,1147297105, 539712780,1545929719, 
     I          190641742,1645390429, 264907697, 620389253,1502074852, 
     J          927711160, 364849192,2049576050, 638580085, 547070247 / 
  
*     Generate the next random number. 
 
      IZI    = IZ(ISTRM) 
      HI15   = IZI / B2E16 
      LOWPRD = (IZI - HI15 * B2E16) * MULT1 
      LOW15  = LOWPRD / B2E16 
      HI31   = HI15 * MULT1 + LOW15 
      OVFLOW = HI31 / B2E15 
      IZI    = (((LOWPRD - LOW15 * B2E16) - MODLUS) + 
     1          (HI31 - OVFLOW * B2E15) * B2E16) + OVFLOW 
      IF (IZI .LT. 0) IZI = IZI + MODLUS 
      HI15   = IZI / B2E16 
      LOWPRD = (IZI - HI15 * B2E16) * MULT2 
      LOW15  = LOWPRD / B2E16 
      HI31   = HI15 * MULT2 + LOW15 
      OVFLOW = HI31 / B2E15 
      IZI    = (((LOWPRD - LOW15 * B2E16) - MODLUS) + 
     1          (HI31 - OVFLOW * B2E15) * B2E16) + OVFLOW 
      IF (IZI .LT. 0) IZI = IZI + MODLUS 
      RAND   = (2 * (IZI / 256) + 1) / 16777216.0 
      IZ(ISTRM) = IZI 
      RETURN 
 
*     Set the current IZ for stream ISTRM to IZIN. 
 
      ENTRY RANDST(IZIN,ISTRM) 
      IZ(ISTRM) = IZIN 
      RETURN 
 
*     Return the current IZ for stream ISTRM in IZOUT.  
 
      ENTRY RANDGT(IZOUT,ISTRM) 
      IZOUT = IZ(ISTRM) 
      RETURN  
 
      END 

