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