MINICERN/mathlib/gen/v/rnormx.F

   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 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(*)
  26       EXTERNAL ROUTIN
  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
  30       DO 200 IDEV = 1, NDEV
  31    50 CALL ROUTIN(U,2)
  32       V = 1.7156 * (U(2) - 0.5)
  33       X = U(1) - S
  34       Y = ABS(V) - T
  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
  43   100 DEVIAT = V/U(1)
  44   200 DEVIAS(IDEV) = DEVIAT
  45       RETURN
  46       END