C ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
FUNCTION RANDN()
C The function RANDN() returns a normally distributed pseudo-random
C number with zero mean and unit variance. Calls are made to a
C function subprogram RANDU() which must return independent random
C numbers uniform in the interval (0,1).
C
C The algorithm uses the ratio of uniforms method of A.J. Kinderman
C and J.F. Monahan augmented with quadratic bounding curves.
C
DATA S,T,A,B / 0.449871, -0.386595, 0.19600, 0.25472/
DATA R1,R2/ 0.27597, 0.27846/
C
C Generate P = (u,v) uniform in rectangle enclosing acceptance region
50 U = RANDU()
V = RANDU()
V = 1.7156 * (V - 0.5)
C Evaluate the quadratic form
X = U - S
Y = ABS(V) - T
Q = X**2 + Y*(A*Y - B*X)
C Accept P if inside inner ellipse
IF (Q .LT. R1) GO TO 100
C Reject P if outside outer ellipse
IF (Q .GT. R2) GO TO 50
C Reject P if outside acceptance region
IF (V**2 .GT. -4.0*ALOG(U)*U**2) GO TO 50
C Return ratio of P's coordinates as the normal deviate
100 RANDN = V/U
RETURN
END