5 * Revision 1.1.1.1 1996/04/01 15:02:55 mclareni
10 SUBROUTINE RNORMX(DEVIAS,NDEV,ROUTIN)
11 C Generator of a vector of independent Gaussian-distributed
12 C (pseudo-)random numbers, of mean zero and variance one,
13 C making use of a uniform pseudo-random generator (RANMAR).
14 C The algorithm for converting uniform numbers to Gaussian
15 C is that of "Ratio of Uniforms with Quadratic Bounds." The
16 C method is in principle exact (apart from rounding errors),
17 C and is based on the variant published by Joseph Leva in
18 C ACM TOMS vol. 18(1992), page 449 for the method and 454 for
19 C the Fortran algorithm (ACM No. 712).
20 C It requires at least 2 and on average 2.74 uniform deviates
21 C per Gaussian (normal) deviate.
22 C WARNING -- The uniform generator should not produce exact zeroes,
23 C since the pair (0.0, 0.5) provokes a floating point exception.
24 SAVE S, T, A, B, R1, R2
25 DIMENSION U(2), DEVIAS(*)
27 DATA S, T, A, B / 0.449871, -0.386595, 0.19600, 0.25472/
28 DATA R1, R2/ 0.27597, 0.27846/
29 C generate pair of uniform deviates
32 V = 1.7156 * (U(2) - 0.5)
35 Q = X**2 + Y*(A*Y - B*X)
36 C accept P if inside inner ellipse
37 IF (Q .LT. R1) GO TO 100
38 C reject P if outside outer ellipse
39 IF (Q .GT. R2) GO TO 50
40 C reject P if outside acceptance region
41 IF (V**2 .GT. -4.0 *ALOG(U(1)) *U(1)**2) GO TO 50
42 C ratio of P's coordinates is normal deviate
44 200 DEVIAS(IDEV) = DEVIAT