C     PORTABLE RANDOM NUMBER GENERATOR (SINGLE PRECISION VERSION)
C     ********************************
C
      FUNCTION RAND (N)
C
C     A PORTABLE, GOOD AND MODERATELY FAST RANDOM NUMBER GENERATOR,
C     PROGRAMMED BY R. P. BRENT, 26/10/1972, REVISED 14/12/76.
C
C     WHEN CALLED WITH N = 0 RAND RETURNS A RANDOM NUMBER UNIFORMLY
C     DISTRIBUTED IN (0, 1)  (INCLUDING 0 BUT NOT 1).   THE GENERATOR
C     MAY BE INITIALIZED BY CALLING WITH NONZERO N (OTHERWISE A DEFAULT
C     INITIALIZATION WITH N = 2038 IS USED).
C
C     IF THE FLOATING-POINT WORD HAS A FRACTION OF IW BITS, THE
C     ALGORITHM RETURNS X(N)/(2**IW), WHERE
C     X(N) = X(N-1) + X(N-127) MODULO 2**IW.   SINCE 1 + X + X**127 IS
C     PRIMITIVE MODULO 2, THE PERIOD IS AT LEAST 2**127 - 1 > 10**38
C     (SEE KNUTH, THE ART OF COMPUTER PROGRAMMING, VOL. 2, 1969,
C     PP. 26, 32-34 AND 464).   X(N) IS STORED AS R3 = X(N)/2**IW.
C
C     THE NUMBER 127 MAY BE REPLACED BY 22, 60, 63 OR 153
C     THROUGHOUT (SEE ZIEGLER AND BRILLHART, ON PRIMITIVE TRINOMIALS
C     (MOD 2), INFORM. CONTR. 13(1968), 541-554 AND 14(1969), 566-569).
C     SUGGESTED VALUES ARE 60 IF THE TABLE SIZE OR INITIALIZATION TIME
C     IS CRITICAL, AND 153 IF THE DEGREE OF RANDOMNESS IS CRITICAL.
C
C     RESTRICTIONS ARE THAT 7*2038 DOES NOT OVERFLOW IN
C     INTEGER ARITHMETIC, AND A BINARY COMPUTER WITH REASONABLE
C     FLOATING-POINT ARITHMETIC IS USED.
C
      IMPLICIT REAL*4 (A-H, O-Z), INTEGER (I-N)
      DIMENSION R3(127)
      DATA R1/2E0/, S/0E0/, T/1E0/, RMC/1E0/, RM/0.015625E0/, IW/-1/
C
C     TEST IF INITIALIZATION IS REQUIRED.
      IF ((R1.LT.1E0).AND.(N.EQ.0)) GO TO 60
C
C     INITIALIZATION IS REQUIRED.  SEE IF IW ETC. HAVE BEEN DETERMINED.
      IF (IW.GT.0) GO TO 30
C
C     DETERMINE FRACTION LENGTH IW ETC.
   10 IW = IW + 1
      T = 0.5E0*T
      R1 = S
      S = S + T
      IF ((S.GT.R1).AND.(S.LT.1E0)) GO TO 10
      IKT = (IW - 1)/6
      IC = IW - 6*IKT
      ID = 2**(11 - IC)
C     COMPUTE RMC = 2**(-IC) EXACTLY.
      DO 20 I = 1, IC
   20 RMC = 0.5E0*RMC
C     INITIALIZE R3 INDEX.
   30 I2 = 127
C     DETERMINE STARTING VALUE FOR SMALL GENERATOR.
      IR = MOD(IABS(N), 2038) + 1
C
C     INITIALIZE R3(127), ... , R3(1).
   40 R1 = 0E0
      DO 50 I = 1, IKT
C     GENERATE RANDOM INTEGER IR IN (0, 2039).
      IR = MOD(7*IR, 2039)
C     ADD SIX RANDOM BITS TO R1.
   50 R1 = (R1 + FLOAT(IR/32))*RM
      IR = MOD (7*IR, 2039)
C     ADD LAST IC RANDOM BITS TO R1.
      R1 = (R1 + FLOAT(IR/ID))*RMC
      R3(I2) = R1
      I2 = I2 - 1
      IF (I2.GT.0) GO TO 40
C     SKIP FIRST 1234 NUMBERS TO AVOID ANY LOCAL NONRANDOMNESS
C     CAUSED BY INITIALIZATION METHOD.
      DO 55 I = 1, 1234
      IF (I2.EQ.0) I2 = 127
      T = R1 + R3(I2)
      IF (T.GE.1E0) T = (R1 - 0.5E0) + (R3(I2) - 0.5E0)
      R1 = T
      R3(I2) = R1
   55 I2 = I2 - 1
C
C     GENERATE NEXT RANDOM NUMBER
   60 IF (I2.EQ.0) I2 = 127
C     ADD LAST AND OLDEST SAVED RANDOM NUMBERS.
      T = R1 + R3(I2)
C     REDUCE MODULO ONE MAKING SURE ARITHMETIC IS EXACT.
      IF (T.GE.1E0) T = (R1 - 0.5E0) + (R3(I2) - 0.5E0)
      R1 = T
C     SAVE FOR FUTURE USE.
      R3(I2) = R1
C     UPDATE COUNTER FOR NEXT CALL.
      I2 = I2 - 1
C     RETURN THE RANDOM NUMBER GENERATED.
      RAND = R1
      RETURN
      END
