DOUBLE PRECISION FUNCTION UNI()
C***BEGIN PROLOGUE  UNI
C***DATE WRITTEN   810915 (YYMMDD)
C***REVISION DATE  871210 (YYMMDD)
C***CATEGORY NO.  L6A21
C***KEYWORDS  RANDOM NUMBERS, UNIFORM RANDOM NUMBERS
C***AUTHOR    KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS
C             MARSAGLIA, GEORGE, SUPERCOMPUTER RES. INST., FLORIDA ST. U.
C
C***PURPOSE  THIS ROUTINE GENERATES REAL (SINGLE PRECISION) UNIFORM 
C             RANDOM NUMBERS ON [0,1)
C***DESCRIPTION
C        Computes real (single precision) uniform numbers on [0,1).
C           From the book, "Numerical Methods and Software" by
C                D. Kahaner, C. Moler, S. Nash
C                Prentice Hall, 1988
C 
C       USAGE: 
C              To initialize the generator 
C                   USEED = USTART(ISEED) 
C               where: ISEED is any NONZERO integer 
C                  will return floating point value of ISEED. 
C 
C               Subsequently 
C                       U = UNI() 
C                  will return a real uniform on [0,1) 
C 
C                One initialization is necessary, but any number of evaluations 
C                  of  UNI in any order, are allowed.
C 
C           Note: Depending upon the value of K (see below), the output
C                       of UNI may differ from one machine to another.
C 
C           Typical usage:
C
C               REAL U,UNI,USTART,USEED
C               INTEGER ISEED
CC                 Set seed
C               ISEED = 305 
C               USEED = USTART(ISEED)
C               DO 1 I = 1,1000 
C                   U = UNI()
C             1 CONTINUE 
CC                 NOTE: If K=24 (the default, see below) the output value of 
CC                           U will be 0.1570390462475...
C               WRITE(*,*) U 
C               END 
C
C          NOTE ON PORTABILITY: Users can choose to run UNI in its default
C               mode (requiring NO user action) which will generate the same 
C               sequence of numbers on any computer supporting floating point
C               numbers with at least 24 bit mantissas, or in a mode that 
C               will generate numbers with a longer period on computers with 
C               larger mantissas.
C          TO EXERCISE THIS OPTION:  B E F O R E  invoking USTART insert
C               the instruction        UBITS = UNIB(K)      K >= 24
C               where K is the number of bits in the mantissa of your floating
C               point word (K=48 for Cray, Cyber 205). UNIB returns the
C               floating point value of K that it actually used.
C                    K input as .LE. 24, then UBITS=24.
C                    K input as .GT. 24, then UBITS=FLOAT(K)
C               If K>24 the sequence of numbers generated by UNI may differ
C               from one computer to another.
C               
C             
C
C***REFERENCES  MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM 
C                 NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U.
C***ROUTINES CALLED  (NONE)
C***END PROLOGUE UNI
      PARAMETER(
     *    CSAVE=362436./16777216.  ,
     *    CD=7654321./16777216.,
     *    CM=16777213./16777216.  )
C                            2**24=16777216
      DOUBLE PRECISION U(17),S,T,USTART,C,UNIB
      INTEGER I,J,II,JJ,K,KK,I1,J1,K1,L1,M1,ISEED
C      
      SAVE U,I,J,K,C
C      Load data array in case user forgets to initialize.
C      This array is the result of calling UNI 100000 times
C         with ISEED=305 and K=64.
      DATA U/
     *0.8668672834288,  0.3697986366357,  0.8008968294805,
     *0.4173889774680,  0.8254561579836,  0.9640965269077,
     *0.4508667414265,  0.6451309529668,  0.1645456024730,
     *0.2787901807898,  0.06761531340295, 0.9663226330820,
     *0.01963343943798, 0.02947398211399, 0.1636231515294,
     *0.3976343250467,  0.2631008574685/
      DATA I,J,K,C/17,5,24,CSAVE/
C 
C   Basic generator is Fibonacci
C 
      UNI = U(I)-U(J)
      IF(UNI.LT.0.0)UNI = UNI+1.0
      U(I) = UNI
      I = I-1
      IF(I.EQ.0)I = 17
      J = J-1
      IF(J.EQ.0)J = 17
C
C   Second generator is congruential
C
      C = C-CD
      IF(C.LT.0.0) C=C+CM
C
C   Combination generator
C
      UNI = UNI-C
      IF(UNI.LT.0.0)UNI = UNI+1.0
      RETURN
C
      ENTRY USTART(ISEED)
C 
C          Set up ...
C          Convert ISEED to four smallish positive integers.
C
        I1 = MOD(ABS(ISEED),177)+1
        J1 = MOD(ABS(ISEED),167)+1
        K1 = MOD(ABS(ISEED),157)+1
        L1 = MOD(ABS(ISEED),147)+1
C
C              Generate random bit pattern in array based on given seed.
C 
        DO 2 II = 1,17
          S = 0.0
          T = 0.5
C             Do for each of the bits of mantissa of word 
C             Loop  over K bits, where K is defaulted to 24 but can
C               be changed by user call to UNIB(K)
          DO 3 JJ = 1,K
                  M1 = MOD(MOD(I1*J1,179)*K1,179)
                  I1 = J1
                  J1 = K1
                  K1 = M1
                  L1 = MOD(53*L1+1,169)
                  IF(MOD(L1*M1,64).GE.32)S=S+T
    3             T = .5*T
    2   U(II) = S
        USTART = FLOAT(ISEED)
        RETURN
C
      ENTRY UNIB(KK)
        IF(KK.LE.24)THEN
             K=24
        ELSE
             K=KK
        ENDIF
        UNIB=FLOAT(K)
      END