      FUNCTION GRAND(N)
C
C     NORMAL RANDOM NUMBER GENERATOR (SINGLE-PRECISION VERSION)
C     ******************************
C
C     WHEN CALLED WITH N = 0 GRAND RETURNS A PSEUDO-RANDOM NUMBER
C     HAVING A NORMAL (I.E., GAUSSIAN) DISTRIBUTION WITH MEAN ZERO
C     AND UNIT STANDARD DEVIATION.  THUS, THE PROBABILITY
C     DENSITY IS F(X) = (2.PI)**(-1/2)*EXP(-X*X/2).
C
C     THE FIRST CALL SHOULD HAVE N NONZERO TO INITIALIZE THE
C     GENERATOR.  IF THE FIRST CALL HAS N = 0 A DEFAULT
C     INITIALIZATION IS USED.
C
C     GRAND CALLS A UNIFORM RANDOM NUMBER GENERATOR RAND(N)
C     WHICH SHOULD RETURN A NUMBER UNIFORMLY DISTRIBUTED IN
C     (0, 1) (POSSIBLY INCLUDING 0), WITH THE SAME INITIALIZATION
C     CONVENTIONS.    GRAND USES THE METHOD OF FORSYTHE AND
C     VON NEUMANN (IMPROVED BY AHRENS, BRENT & DIETER), AND
C     CALLS RAND 1.37746 TIMES ON THE AVERAGE.
C
C     R. BRENT, 7/11/1972, LAST MODIFIED 2/7/1973.
C
      IMPLICIT REAL*4 (A-H, O-Z), INTEGER (I-N)
C     DIMENSION OF D IS AT LEAST THE NUMBER OF BITS (IW) IN THE
C     FLOATING-POINT FRACTION.
      DIMENSION D(27)
C     IF THE INTEGRAL OF (2/PI)**0.5*EXP(-X**2/2) FROM A(I) TO
C     INFINITY IS 2**(-I), THEN D(I) = A(I) - A(I-1), 0 < I <= IW.
      DATA D /0.674489750, 0.475859630, 0.383771164, 0.328611323,
     2        0.291142827, 0.263684322, 0.242508452, 0.225567444,
     3        0.211634166, 0.199924267, 0.189910758, 0.181225181,
     4        0.173601400, 0.166841909, 0.160796729, 0.155349717,
     5        0.150409384, 0.145902577, 0.141770033, 0.137963174,
     6        0.134441762, 0.131172150, 0.128125965, 0.125279090,
     7        0.122610883, 0.120103560, 0.117741707/
      DATA U /-1E0/
C     INITIALIZE U IF NECESSARY.
      IF ((U.LT.0E0).OR.(N.NE.0)) U = RAND(N)
C     INITIALIZE DISPLACEMENT A AND COUNTER I.
      A = 0E0
      I = 0
C     INCREMENT COUNTER AND DISPLACEMENT IF LEADING BIT OF U IS 1.
   10 U = U + U
      IF (U.LT.1E0) GO TO 30
   20 U = U - 1E0
      I = I + 1
      A = A - D(I)
      GO TO 10
C     FORM W UNIFORM ON 0 <= W < D(I+1) FROM U.
   30 W = D(I+1)*U
C     FORM V = ((W-A)**2 - A**2)/2.  (NOTE THAT 0 <= V < LN2.)
      V = W*(0.5E0*W - A)
C     GENERATE NEW UNIFORM U.
   40 U = RAND(0)
C     ACCEPT W AS A RANDOM SAMPLE IF V <= U.
      IF (V.LE.U) GO TO 50
C     GENERATE RANDOM V.
      V = RAND(0)
C     LOOP IF U > V.
      IF (U.GT.V) GO TO 40
C     REJECT W AND FORM A NEW UNIFORM U FROM V AND U.
      U = (V - U)/(1E0 - U)
      GO TO 30
C     FORM NEW U (TO BE USED ON NEXT CALL) FROM U AND V.
   50 U = (U - V)/(1E0 - V)
C     USE FIRST BIT OF U FOR SIGN, RETURN NORMAL VARIATE.
      U = U + U
      IF (U.LT.1E0) GO TO 60
      U = U - 1E0
      GRAND = W - A
      RETURN
   60 GRAND = A - W
      RETURN
      END
