]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1996/04/01 15:02:55 mclareni | |
6 | * Mathlib gen | |
7 | * | |
8 | * | |
9 | #include "gen/pilot.h" | |
10 | SUBROUTINE RNORML(DEVIAS,NDEV) | |
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(*) | |
26 | DATA S, T, A, B / 0.449871, -0.386595, 0.19600, 0.25472/ | |
27 | DATA R1, R2/ 0.27597, 0.27846/ | |
28 | C generate pair of uniform deviates | |
29 | DO 200 IDEV = 1, NDEV | |
30 | 50 CALL RANMAR(U,2) | |
31 | V = 1.7156 * (U(2) - 0.5) | |
32 | X = U(1) - S | |
33 | Y = ABS(V) - T | |
34 | Q = X**2 + Y*(A*Y - B*X) | |
35 | C accept P if inside inner ellipse | |
36 | IF (Q .LT. R1) GO TO 100 | |
37 | C reject P if outside outer ellipse | |
38 | IF (Q .GT. R2) GO TO 50 | |
39 | C reject P if outside acceptance region | |
40 | IF (V**2 .GT. -4.0 *ALOG(U(1)) *U(1)**2) GO TO 50 | |
41 | C ratio of P's coordinates is normal deviate | |
42 | 100 DEVIAT = V/U(1) | |
43 | 200 DEVIAS(IDEV) = DEVIAT | |
44 | RETURN | |
45 | END |