C ALGORITHM 488 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN COMM. ACM, VOL. 17, NO. 12, C P. 704. FUNCTION GRAND(N) GRA 10 C EXCEPT ON THE FIRST CALL GRAND RETURNS A C PSEUDO-RANDOM NUMBER HAVING A GAUSSIAN (I.E. C NORMAL) DISTRIBUTION WITH ZERO MEAN AND UNIT C STANDARD DEVIATION. THUS, THE DENSITY IS F(X) = C EXP(-0.5*X**2)/SQRT(2.0*PI). THE FIRST CALL C INITIALIZES GRAND AND RETURNS ZERO. C THE PARAMETER N IS DUMMY. C GRAND CALLS A FUNCTION RAND, AND IT IS ASSUMED THAT C SUCCESSIVE CALLS TO RAND(0) GIVE INDEPENDENT C PSEUDO- RANDOM NUMBERS DISTRIBUTED UNIFORMLY ON (0, C 1), POSSIBLY INCLUDING 0 (BUT NOT 1). C THE METHOD USED WAS SUGGESTED BY VON NEUMANN, AND C IMPROVED BY FORSYTHE, AHRENS, DIETER AND BRENT. C ON THE AVERAGE THERE ARE 1.37746 CALLS OF RAND FOR C EACH CALL OF GRAND. C WARNING - DIMENSION AND DATA STATEMENTS BELOW ARE C MACHINE-DEPENDENT. C DIMENSION OF D MUST BE AT LEAST THE NUMBER OF BITS C IN THE FRACTION OF A FLOATING-POINT NUMBER. C THUS, ON MOST MACHINES THE DATA STATEMENT BELOW C CAN BE TRUNCATED. C IF THE INTEGRAL OF SQRT(2.0/PI)*EXP(-0.5*X**2) FROM C A(I) TO INFINITY IS 2**(-I), THEN D(I) = A(I) - C A(I-1). DIMENSION D(60) DATA D(1), D(2), D(3), D(4), D(5), D(6), D(7), * D(8), D(9), D(10), D(11), D(12), D(13), * D(14), D(15), D(16), D(17), D(18), D(19), * D(20), D(21), D(22), D(23), D(24), D(25), * D(26), D(27), D(28), D(29), D(30), D(31), * D(32) /0.674489750,0.475859630,0.383771164, * 0.328611323,0.291142827,0.263684322, * 0.242508452,0.225567444,0.211634166, * 0.199924267,0.189910758,0.181225181, * 0.173601400,0.166841909,0.160796729, * 0.155349717,0.150409384,0.145902577, * 0.141770033,0.137963174,0.134441762, * 0.131172150,0.128125965,0.125279090, * 0.122610883,0.120103560,0.117741707, * 0.115511892,0.113402349,0.111402720, * 0.109503852,0.107697617/ DATA D(33), D(34), D(35), D(36), D(37), D(38), * D(39), D(40), D(41), D(42), D(43), D(44), * D(45), D(46), D(47), D(48), D(49), D(50), * D(51), D(52), D(53), D(54), D(55), D(56), * D(57), D(58), D(59), D(60) * /0.105976772,0.104334841,0.102766012, * 0.101265052,0.099827234,0.098448282, * 0.097124309,0.095851778,0.094627461, * 0.093448407,0.092311909,0.091215482, * 0.090156838,0.089133867,0.088144619, * 0.087187293,0.086260215,0.085361834, * 0.084490706,0.083645487,0.082824924, * 0.082027847,0.081253162,0.080499844, * 0.079766932,0.079053527,0.078358781, * 0.077681899/ C END OF MACHINE-DEPENDENT STATEMENTS C U MUST BE PRESERVED BETWEEN CALLS. DATA U /0.0/ C INITIALIZE DISPLACEMENT A AND COUNTER I. A = 0.0 I = 0 C INCREMENT COUNTER AND DISPLACEMENT IF LEADING BIT C OF U IS ONE. 10 U = U + U IF (U.LT.1.0) GO TO 20 U = U - 1.0 I = I + 1 A = A - D(I) GO TO 10 C FORM W UNIFORM ON 0 .LE. W .LT. D(I+1) FROM U. 20 W = D(I+1)*U C FORM V = 0.5*((W-A)**2 - A**2). NOTE THAT 0 .LE. V C .LT. LOG(2). V = W*(0.5*W-A) C GENERATE NEW UNIFORM U. 30 U = RAND(0) C ACCEPT W AS A RANDOM SAMPLE IF V .LE. U. IF (V.LE.U) GO TO 40 C GENERATE RANDOM V. V = RAND(0) C LOOP IF U .GT. V. IF (U.GT.V) GO TO 30 C REJECT W AND FORM A NEW UNIFORM U FROM V AND U. U = (V-U)/(1.0-U) GO TO 20 C FORM NEW U (TO BE USED ON NEXT CALL) FROM U AND V. 40 U = (U-V)/(1.0-V) C USE FIRST BIT OF U FOR SIGN, RETURN NORMAL VARIATE. U = U + U IF (U.LT.1.0) GO TO 50 U = U - 1.0 GRAND = W - A RETURN 50 GRAND = A - W RETURN END