[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / v / rnormx.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:55  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
      SUBROUTINE RNORMX(DEVIAS,NDEV,ROUTIN)
C        Generator of a vector of independent Gaussian-distributed 
C        (pseudo-)random numbers, of mean zero and variance one,
C        making use of a uniform pseudo-random generator (RANMAR).
C        The algorithm for converting uniform numbers to Gaussian
C        is that of "Ratio of Uniforms with Quadratic Bounds."  The
C        method is in principle exact (apart from rounding errors),
C        and is based on the variant published by Joseph Leva in
C        ACM TOMS vol. 18(1992), page 449 for the method and 454 for
C        the Fortran algorithm (ACM No. 712).
C        It requires at least 2 and on average 2.74 uniform deviates
C        per Gaussian (normal) deviate.
C   WARNING -- The uniform generator should not produce exact zeroes,
C   since the pair (0.0, 0.5) provokes a floating point exception.
      SAVE  S, T, A, B, R1, R2
      DIMENSION U(2), DEVIAS(*)
      EXTERNAL ROUTIN
      DATA  S, T, A, B / 0.449871, -0.386595, 0.19600, 0.25472/
      DATA  R1, R2/ 0.27597, 0.27846/
C         generate pair of uniform deviates
      DO 200 IDEV = 1, NDEV
   50 CALL ROUTIN(U,2)
      V = 1.7156 * (U(2) - 0.5)
      X = U(1) - 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(1)) *U(1)**2)  GO TO 50
C           ratio of P's coordinates is normal deviate
  100 DEVIAT = V/U(1)
  200 DEVIAS(IDEV) = DEVIAT
      RETURN
      END
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 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(AY - BX)
	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