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.
PROGRAM RTEST
C This program is the top level driver to test the normal random
C number generator. The subprogram RANDN() implements the suggested
C algorithm. It uses the ratio of uniforms method with quadratic
C bounding curves. A second subprogram RANDN0() implements the ratio
C of uniforms method with no auxiliary curves. Each of these
C algorithms needs a supply of uniform random numbers. These are
C provided using table lookup with the subprogram RANDU(). Given
C the same sequence of uniform numbers, the two algorithms should
C return the same pseudo-random sequence of normal deviates.
C
C The test program generates 100 normal deviates using each of the
C algorithms with the same sequence of uniform numbers. The normal
C deviates are then compared and a message is printed regarding
C their equality. The normal deviates are finally printed.
C
C The named common RANUC is shared between RTEST and RANDU().
C The uniform generator RANDU() is initialized by setting INDX=0
DIMENSION R(100), R0(100)
COMMON /RANUC/INDX
C initialize uniform number generator and generate 100 normal
C deviates with suggested algorithm
INDX = 0
DO 20 I=1,100
20 R(I) = RANDN()
C re-initialize and repeat with alternate algorithm
INDX = 0
DO 40 I=1,100
40 R0(I) = RANDN0()
C compare sequences for equality
ICOUNT = 0
DO 60 I=1,100
IF (ABS(R(I) - R0(I)) .GT. 1.0E-6) ICOUNT=ICOUNT+1
60 CONTINUE
C print out results
WRITE(6,100)
IF (ICOUNT .EQ. 0) WRITE(6,110)
IF (ICOUNT .NE. 0) WRITE(6,120) ICOUNT
WRITE(6,130) (R(I),I=1,100)
WRITE(6,140) INDX
C
100 FORMAT(/' The routine RANDN (suggested algorithm using',
1 ' quadratic boundary'/' curves) and RANDN0 (ratio of uniforms',
2 ' method without auxiliary'/' curves) each generated 100',
3 ' normal deviates using the same'/' sequence of uniform',
4 ' numbers.')
110 FORMAT(/' As expected, the routines generated identical',
1 ' results.')
120 FORMAT(/' There is some discrepancy. The resulting',
1 ' sequences differ'/' in ', I4, ' places.')
130 FORMAT(/' The pseudo-random deviates generated by',
1 ' RANDN follow:'/ (8F8.4) )
140 FORMAT(/' The 100 normal deviates required',I4,
1 ' calls to the uniform'/' number generator.')
END
C********************
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
C********************
FUNCTION RANDN0()
C
C This function generates a normally distributed number using the
C ratio of uniforms method. No auxiliary boundary curves are used.
C Calls to the subprogram RANDU() must return independent random
C numbers uniform in the interval (0,1).
C
C Generate P = (u,v) uniform in rectangle enclosing acceptance region
50 U = RANDU()
V = RANDU()
V = 1.7156 * (V - 0.5)
C Compute coordinate ratio for candidate point
R = V/U
C Reject point if outside acceptance region
IF (R**2 .GT. -4.0*ALOG(U)) GO TO 50
C Return ratio as the normal deviate
RANDN0 = R
RETURN
END
C********************
FUNCTION RANDU()
C
C This is a dummy routine for generated uniform numbers in the
C interval (0,1). A list of stored numbers is sequentially
C accessed and returned. The routine is supplied to permit testing
C of the subprograms RANDN() and RAND0().
C
C The routine cycles through 279 values stored in the data array U.
C The variable INDX in the named common RANUC retains the index of the
C number returned. This variable can be initialized to 0 in a calling
C routine to restart the sequence.
C
DIMENSION U(279)
COMMON /RANUC/INDX
C
DATA U(1),U(2),U(3),U(4),U(5),U(6),U(7),U(8),U(9),U(10),U(11),
* U(12),U(13),U(14),U(15),U(16),U(17),U(18),U(19),U(20),U(21),
* U(22),U(23),U(24),U(25),U(26),U(27),U(28),U(29),U(30),U(31),
* U(32),U(33),U(34),U(35),U(36),U(37),U(38),U(39),U(40),U(41),
* U(42),U(43),U(44),U(45),U(46),U(47),U(48),U(49),U(50),U(51),
* U(52),U(53),U(54),U(55),U(56),U(57),U(58),U(59),U(60),U(61),
* U(62),U(63),U(64),U(65),U(66),U(67),U(68),U(69),U(70)/
* 0.5139,0.1757,0.3087,0.5345,0.9476,0.1717,0.7022,0.2264,0.4948,
* 0.1247,0.0839,0.3896,0.2772,0.3681,0.9834,0.5354,0.7657,0.6465,
* 0.7671,0.7802,0.8230,0.1519,0.6255,0.3147,0.3469,0.9172,0.5198,
* 0.4012,0.6068,0.7854,0.9315,0.8699,0.8665,0.6745,0.7584,0.5819,
* 0.3892,0.3556,0.2002,0.8269,0.4159,0.4635,0.9792,0.1264,0.2126,
* 0.9585,0.7375,0.4091,0.7801,0.7579,0.9568,0.0281,0.3187,0.7569,
* 0.2430,0.5895,0.0434,0.9560,0.3191,0.0594,0.4419,0.9150,0.5722,
* 0.1188,0.5698,0.2520,0.4959,0.2367,0.4770,0.4061/
C
DATA U(71),U(72),U(73),U(74),U(75),U(76),U(77),U(78),U(79),U(80),
* U(81),U(82),U(83),U(84),U(85),U(86),U(87),U(88),U(89),U(90),
* U(91),U(92),U(93),U(94),U(95),U(96),U(97),U(98),U(99),U(100),
* U(101),U(102),U(103),U(104),U(105),U(106),U(107),U(108),U(109),
* U(110),U(111),U(112),U(113),U(114),U(115),U(116),U(117),U(118),
* U(119),U(120),U(121),U(122),U(123),U(124),U(125),U(126),U(127),
* U(128),U(129),U(130),U(131),U(132),U(133),U(134),U(135),U(136),
* U(137),U(138),U(139),U(140)/
* 0.8730,0.4270,0.3582,0.3820,0.0432,0.1606,0.5224,0.6966,0.0971,
* 0.4008,0.7734,0.2448,0.3428,0.2300,0.2979,0.3045,0.8872,0.0367,
* 0.6511,0.3986,0.6763,0.7326,0.9378,0.2333,0.8385,0.9672,0.7786,
* 0.4315,0.6741,0.8094,0.1588,0.2799,0.1353,0.8642,0.7502,0.2080,
* 0.1400,0.2946,0.8028,0.2189,0.5631,0.7156,0.1975,0.9898,0.2500,
* 0.4306,0.7553,0.8609,0.8948,0.9781,0.3954,0.4322,0.1271,0.4577,
* 0.2378,0.9860,0.6528,0.6042,0.2419,0.4549,0.7900,0.0788,0.4764,
* 0.1526,0.2458,0.9450,0.6140,0.9882,0.4773,0.7997/
C
DATA U(141),U(142),U(143),U(144),U(145),U(146),U(147),U(148),
* U(149),U(150),U(151),U(152),U(153),U(154),U(155),U(156),U(157),
* U(158),U(159),U(160),U(161),U(162),U(163),U(164),U(165),U(166),
* U(167),U(168),U(169),U(170),U(171),U(172),U(173),U(174),U(175),
* U(176),U(177),U(178),U(179),U(180),U(181),U(182),U(183),U(184),
* U(185),U(186),U(187),U(188),U(189),U(190),U(191),U(192),U(193),
* U(194),U(195),U(196),U(197),U(198),U(199),U(200),U(201),U(202),
* U(203),U(204),U(205),U(206),U(207),U(208),U(209),U(210)/
* 0.7442,0.3807,0.4799,0.5269,0.0981,0.5942,0.3472,0.1434,0.7795,
* 0.7110,0.4461,0.7046,0.0953,0.9628,0.5513,0.7403,0.5790,0.6379,
* 0.7817,0.1879,0.3021,0.2828,0.6840,0.2929,0.5654,0.4184,0.3066,
* 0.4445,0.5657,0.4879,0.6066,0.4159,0.1304,0.2560,0.0358,0.9771,
* 0.1145,0.3781,0.6467,0.3504,0.5530,0.3584,0.5655,0.4756,0.1637,
* 0.6152,0.1722,0.5547,0.2922,0.8722,0.8351,0.8449,0.8955,0.5948,
* 0.5406,0.1682,0.6550,0.6905,0.2639,0.1067,0.8149,0.1914,0.4233,
* 0.3519,0.8392,0.1373,0.2627,0.1773,0.4799,0.3802/
C
DATA U(211),U(212),U(213),U(214),U(215),U(216),U(217),U(218),
* U(219),U(220),U(221),U(222),U(223),U(224),U(225),U(226),U(227),
* U(228),U(229),U(230),U(231),U(232),U(233),U(234),U(235),U(236),
* U(237),U(238),U(239),U(240),U(241),U(242),U(243),U(244),U(245),
* U(246),U(247),U(248),U(249),U(250),U(251),U(252),U(253),U(254),
* U(255),U(256),U(257),U(258),U(259),U(260),U(261),U(262),U(263),
* U(264),U(265),U(266),U(267),U(268),U(269),U(270),U(271),U(272),
* U(273),U(274),U(275),U(276),U(277),U(278),U(279)/
* 0.5048,0.5028,0.3519,0.5256,0.1206,0.5196,0.6071,0.7329,0.5569,
* 0.3441,0.8020,0.5910,0.2669,0.6707,0.5522,0.7889,0.8877,0.8900,
* 0.0681,0.8006,0.9074,0.6441,0.1652,0.3014,0.1663,0.2852,0.8420,
* 0.5363,0.0363,0.2072,0.0212,0.3581,0.6215,0.5200,0.5460,0.1537,
* 0.8234,0.0334,0.0260,0.3781,0.6163,0.0204,0.6266,0.9152,0.3748,
* 0.7295,0.3958,0.9823,0.5973,0.1123,0.2216,0.7992,0.8707,0.7382,
* 0.0136,0.7396,0.4184,0.3620,0.2039,0.1832,0.0763,0.1156,0.1591,
* 0.7883,0.0404,0.7906,0.5990,0.4026,0.2291/
C
INDX = INDX + 1
IF (INDX .GT. 279) INDX = 1
RANDU = U(INDX)
RETURN
END
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