C     ALGORITHM 488 COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN COMM. ACM, VOL. 17, NO. 12,
C     P. 704.
      FUNCTION GRAND(N)                                                 GRA   10
C EXCEPT ON THE FIRST CALL GRAND RETURNS A
C PSEUDO-RANDOM NUMBER HAVING A GAUSSIAN (I.E.
C NORMAL) DISTRIBUTION WITH ZERO MEAN AND UNIT
C STANDARD DEVIATION.  THUS, THE DENSITY IS  F(X) =
C EXP(-0.5*X**2)/SQRT(2.0*PI). THE FIRST CALL
C INITIALIZES GRAND AND RETURNS ZERO.
C THE PARAMETER N IS DUMMY.
C GRAND CALLS A FUNCTION RAND, AND IT IS ASSUMED THAT
C SUCCESSIVE CALLS TO RAND(0) GIVE INDEPENDENT
C PSEUDO- RANDOM NUMBERS DISTRIBUTED UNIFORMLY ON (0,
C 1), POSSIBLY INCLUDING 0 (BUT NOT 1).
C THE METHOD USED WAS SUGGESTED BY VON NEUMANN, AND
C IMPROVED BY FORSYTHE, AHRENS, DIETER AND BRENT.
C ON THE AVERAGE THERE ARE 1.37746 CALLS OF RAND FOR
C EACH CALL OF GRAND.
C WARNING - DIMENSION AND DATA STATEMENTS BELOW ARE
C           MACHINE-DEPENDENT.
C DIMENSION OF D MUST BE AT LEAST THE NUMBER OF BITS
C IN THE FRACTION OF A FLOATING-POINT NUMBER.
C THUS, ON MOST MACHINES THE DATA STATEMENT BELOW
C CAN BE TRUNCATED.
C IF THE INTEGRAL OF SQRT(2.0/PI)*EXP(-0.5*X**2) FROM
C A(I) TO INFINITY IS 2**(-I), THEN D(I) = A(I) -
C A(I-1).
      DIMENSION D(60)
      DATA D(1), D(2), D(3), D(4), D(5), D(6), D(7),
     * D(8), D(9), D(10), D(11), D(12), D(13),
     * D(14), D(15), D(16), D(17), D(18), D(19),
     * D(20), D(21), D(22), D(23), D(24), D(25),
     * D(26), D(27), D(28), D(29), D(30), D(31),
     * D(32) /0.674489750,0.475859630,0.383771164,
     * 0.328611323,0.291142827,0.263684322,
     * 0.242508452,0.225567444,0.211634166,
     * 0.199924267,0.189910758,0.181225181,
     * 0.173601400,0.166841909,0.160796729,
     * 0.155349717,0.150409384,0.145902577,
     * 0.141770033,0.137963174,0.134441762,
     * 0.131172150,0.128125965,0.125279090,
     * 0.122610883,0.120103560,0.117741707,
     * 0.115511892,0.113402349,0.111402720,
     * 0.109503852,0.107697617/
      DATA D(33), D(34), D(35), D(36), D(37), D(38),
     * D(39), D(40), D(41), D(42), D(43), D(44),
     * D(45), D(46), D(47), D(48), D(49), D(50),
     * D(51), D(52), D(53), D(54), D(55), D(56),
     * D(57), D(58), D(59), D(60)
     * /0.105976772,0.104334841,0.102766012,
     * 0.101265052,0.099827234,0.098448282,
     * 0.097124309,0.095851778,0.094627461,
     * 0.093448407,0.092311909,0.091215482,
     * 0.090156838,0.089133867,0.088144619,
     * 0.087187293,0.086260215,0.085361834,
     * 0.084490706,0.083645487,0.082824924,
     * 0.082027847,0.081253162,0.080499844,
     * 0.079766932,0.079053527,0.078358781,
     * 0.077681899/
C END OF MACHINE-DEPENDENT STATEMENTS
C U MUST BE PRESERVED BETWEEN CALLS.
      DATA U /0.0/
C INITIALIZE DISPLACEMENT A AND COUNTER I.
      A = 0.0
      I = 0
C INCREMENT COUNTER AND DISPLACEMENT IF LEADING BIT
C OF U IS ONE.
   10 U = U + U
      IF (U.LT.1.0) GO TO 20
      U = U - 1.0
      I = I + 1
      A = A - D(I)
      GO TO 10
C FORM W UNIFORM ON 0 .LE. W .LT. D(I+1) FROM U.
   20 W = D(I+1)*U
C FORM V = 0.5*((W-A)**2 - A**2). NOTE THAT 0 .LE. V
C .LT. LOG(2).
      V = W*(0.5*W-A)
C GENERATE NEW UNIFORM U.
   30 U = RAND(0)
C ACCEPT W AS A RANDOM SAMPLE IF V .LE. U.
      IF (V.LE.U) GO TO 40
C GENERATE RANDOM V.
      V = RAND(0)
C LOOP IF U .GT. V.
      IF (U.GT.V) GO TO 30
C REJECT W AND FORM A NEW UNIFORM U FROM V AND U.
      U = (V-U)/(1.0-U)
      GO TO 20
C FORM NEW U (TO BE USED ON NEXT CALL) FROM U AND V.
   40 U = (U-V)/(1.0-V)
C USE FIRST BIT OF U FOR SIGN, RETURN NORMAL VARIATE.
      U = U + U
      IF (U.LT.1.0) GO TO 50
      U = U - 1.0
      GRAND = W - A
      RETURN
   50 GRAND = A - W
      RETURN
      END