

      FUNCTION GRN(KSEED)

C  GRN - Gaussian (Normal) Random Number Generator
C
C  Returns a normal random number with mean zero
C  and standard deviation one.
C
C  Algorithm 488:  Collected Algorithms ACM
C
C  Reference:
C
C    Brent: A Gaussian pseudo-random number
C    generator, Comm ACM, vol 17 #12, Dec 1974,
C    704-706.
C
C  Implemented:
C
C    R Kamen, Code 3873, Sept 81
C
C  Modified:
C
C    K Lindblom, Code 3873, per
C    L Lucas, Code 3873, 21 Jan 82
C
C  GRN is Brent's GRAND modified as follows:
C
C    1. The dummy argument N is replaced by
C       the argument KSEED, which is used to
C       hold the 'memory' of the uniform
C       generator URN.
C
C    2. References to RAND(0) are replaced by
C       references to URN(KSEED).
C
C    3. Dimension of D is reduced from 60 to 30.
C
C  Except on the first call to GRN...
C  GRN returns a pseudo-random number having a
C  Gaussian (i.e. normal) distribution with zero
C  mean and unit variance.  Thus, the probability
C  density is
C
C    f(x) = exp(- x**2 / 2.) / sqrt(2.*pi)
C
C  The first call initializes GRN and returns 0.
C
C  When GRN is called, KSEED must be the name of
C  an integer variable in the calling program.
C  Initially, KSEED must be an integer between
C  1 and 2147483646.  KSEED should not be changed
C  between calls to GRN.
C
C  GRN assumes that successive calls to URN yield
C  independent, pseudo-random numbers distributed
C  uniformly on (0,1) - possibly including 0 but
C  NOT 1.
C
C  The method used was suggested by von Newmann
C  and improved by Forsythe, Ahrens, Dieter and
C  Brent.
C
C  On the average, there are 1.37746 calls to URN
C  for each call to GRN.
C
C
C
C  Warning - dimensions and data statements below
C            are machine-dependent.
C
C  Dimension of D must be at least the number of
C  bits in the fraction of a floating-point number.
C  Thus, on most machines the DATA statement below
C  may be truncated.
C
C  If the integral of sqrt(2./pi)*exp(- x**2 / 2)
C  from A(I) to infinity is 2**(-I), then
C
C       D(I) = A(I) - A(I-1)
C
      EXTERNAL URN
c
c     DIMENSION D(60)
      DIMENSION  D(30)
      DATA D(1)   /0.674489750/
      DATA D(2)   /0.475859630/
      DATA D(3)   /0.383771164/
      DATA D(4)   /0.328611323/
      DATA D(5)   /0.291142827/
      DATA D(6)   /0.263684322/
      DATA D(7)   /0.242508452/
      DATA D(8)   /0.225567444/
      DATA D(9)   /0.211634166/
      DATA D(10)  /0.199924267/
      DATA D(11)  /0.189910758/
      DATA D(12)  /0.181225181/
      DATA D(13)  /0.173601400/
      DATA D(14)  /0.166841909/
      DATA D(15)  /0.160796729/
      DATA D(16)  /0.155349717/
      DATA D(17)  /0.150409384/
      DATA D(18)  /0.145902577/
      DATA D(19)  /0.141770033/
      DATA D(20)  /0.137963174/
      DATA D(21)  /0.134441762/
      DATA D(22)  /0.131172150/
      DATA D(23)  /0.128125965/
      DATA D(24)  /0.125279090/
      DATA D(25)  /0.122610883/
      DATA D(26)  /0.120103560/
      DATA D(27)  /0.117741707/
      DATA D(28)  /0.115511892/
      DATA D(29)  /0.113402349/
      DATA D(30)  /0.111402720/
c     DATA D(31)  /0.109503852/
c     DATA D(32)  /0.107697617/
c     DATA D(33)  /0.105976772/
c     DATA D(34)  /0.104334841/
c     DATA D(35)  /0.102766012/
c     DATA D(36)  /0.101265052/
c     DATA D(37)  /0.099827234/
c     DATA D(38)  /0.098448282/
c     DATA D(39)  /0.097124309/
c     DATA D(40)  /0.095851778/
c     DATA D(41)  /0.094627461/
c     DATA D(42)  /0.093448407/
c     DATA D(43)  /0.092311909/
c     DATA D(44)  /0.091215482/
c     DATA D(45)  /0.090156838/
c     DATA D(46)  /0.089133867/
c     DATA D(47)  /0.088144619/
c     DATA D(48)  /0.087187293/
c     DATA D(49)  /0.086260215/
c     DATA D(50)  /0.085361834/
c     DATA D(51)  /0.084490706/
c     DATA D(52)  /0.083645487/
c     DATA D(53)  /0.082824924/
c     DATA D(54)  /0.082027847/
c     DATA D(55)  /0.081253162/
c     DATA D(56)  /0.080499844/
c     DATA D(57)  /0.079766932/
c     DATA D(58)  /0.079053527/
c     DATA D(59)  /0.078358781/
c     DATA D(60)  /0.077681899/
C
C  End of machine-dependent statements.
C
C  U must be preserved between calls.
C
       DATA        U / 0.0 /
C
C  Initialize displacement A and counter I
C
       A = 0.0
       I = 0
C
C  Increment counter and displacement if
C  leading bit of U is one.
C
   10  U = U + U
       IF ( U .LT. 1.0 )  GOTO 20
       U = U - 1.0
       I = I + 1
       A = A - D(I)
       GOTO 10
C
C  Form W uniform on 0 .LE. W .LT. D(I+1) from U
C
   20  W = D(I+1) * U
C
C  Form V = 0.5 * ((W - A)**2 - A**2)
C
C  Note:             0 .LE. V .LT. LOG(2)
C
       V = W * (0.5*W - A)
C
C  Generate new uniform U
C
   30  U = URN(KSEED)
C
C  Accept W as a random sample if V .LE. U
C
       IF ( V .LE. U )  GOTO 40
C
C  Generate random V
C
       V = URN(KSEED)
C
C  Loop if U .GT. V
C
       IF ( U .GT. V )  GOTO 30
C
C  Reject W and form a new uniform U from V and U
C
       U = (V - U) / (1.0 - U)
       GOTO 20
C
C  Form new U (to be used on next call) from U and V
C
   40  U = (U - V) / (1.0 - V)
C
C  Use first bit of U for sign, return normal variate
C
       U = U + U
       IF ( U .LT. 1.0 )  GOTO 50
       U = U - 1.0
       GRN = W - A
       RETURN
C
   50  GRN = A - W
       RETURN
       END





