DOUBLE PRECISION FUNCTION DUNI() C***BEGIN PROLOGUE DUNI C***DATE WRITTEN 880714 (YYMMDD) C***REVISION DATE 880714 (YYMMDD) C***CATEGORY NO. L6A21 C***KEYWORDS RANDOM NUMBERS, UNIFORM RANDOM NUMBERS C***AUTHOR KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS C MARSAGLIA, GEORGE, SUPERCOMPUTER RES. INST., FLORIDA ST. U. C C***PURPOSE THIS ROUTINE GENERATES DOUBLE PRECISION UNIFORM C RANDOM NUMBERS ON [0,1) C***DESCRIPTION C COMPUTES DOUBLE PRECISION UNIFORM NUMBERS ON [0,1). C FROM THE BOOK, "NUMERICAL METHODS AND SOFTWARE" BY C D. KAHANER, C. MOLER, S. NASH C PRENTICE HALL, 1988 C C USAGE: C TO INITIALIZE THE GENERATOR C USEED = DUSTAR(ISEED) C WHERE: ISEED IS ANY NONZERO INTEGER C WILL RETURN FLOATING POINT VALUE OF ISEED. C C SUBSEQUENTLY C U = DUNI() C WILL RETURN A REAL UNIFORM ON [0,1) C C ONE INITIALIZATION IS NECESSARY, BUT ANY NUMBER OF EVALUATIONS C OF DUNI IN ANY ORDER, ARE ALLOWED. C C NOTE: DEPENDING UPON THE VALUE OF K (SEE BELOW), THE OUTPUT C OF DUNI MAY DIFFER FROM ONE MACHINE TO ANOTHER. C C TYPICAL USAGE: C C DOUBLE PRECISION U,DUNI,DUSTAR,USEED C INTEGER ISEED CC SET SEED C ISEED = 305 C USEED = DUSTAR(ISEED) C DO 1 I = 1,1000 C U = DUNI() C 1 CONTINUE CC NOTE: IF K=47 (THE DEFAULT, SEE BELOW) THE OUTPUT VALUE OF CC U WILL BE 0.812053811384E-01... C WRITE(*,*) U C END C C NOTE ON PORTABILITY: USERS CAN CHOOSE TO RUN DUNI IN ITS DEFAULT C MODE (REQUIRING NO USER ACTION) WHICH WILL GENERATE THE SAME C SEQUENCE OF NUMBERS ON ANY COMPUTER SUPPORTING FLOATING POINT C NUMBERS WITH AT LEAST 47 BIT MANTISSAS, OR IN A MODE THAT C WILL GENERATE NUMBERS WITH A LONGER PERIOD ON COMPUTERS WITH C LARGER MANTISSAS. C TO EXERCISE THIS OPTION: B E F O R E INVOKING DUSTAR INSERT C THE INSTRUCTION UBITS = DUNIB(K) K >= 47 C WHERE K IS THE NUMBER OF BITS IN THE MANTISSA OF YOUR FLOATING C POINT WORD (K=96 FOR CRAY, CYBER 205). DUNIB RETURNS THE C FLOATING POINT VALUE OF K THAT IT ACTUALLY USED. C K INPUT AS .LE. 47, THEN UBITS=47. C K INPUT AS .GT. 47, THEN UBITS=FLOAT(K) C IF K>47 THE SEQUENCE OF NUMBERS GENERATED BY DUNI MAY DIFFER C FROM ONE COMPUTER TO ANOTHER. C C C C***REFERENCES MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM C NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U. C***ROUTINES CALLED (NONE) C***END PROLOGUE DUNI DOUBLE PRECISION CSAVE,CD,CM PARAMETER( * CSAVE=0.9162596898123D+13/0.140737488355328D+15, * CD=0.76543212345678D+14/0.140737488355328D+15, * CM=0.140737488355213D+15/0.140737488355328D+15) C 2**47=0.140737488355328D+15 DOUBLE PRECISION U(17),S,T,DUSTAR,C,DUNIB INTEGER I,J,II,JJ,K,KK,I1,J1,K1,L1,M1,ISEED C SAVE U,I,J,K,C C LOAD DATA ARRAY IN CASE USER FORGETS TO INITIALIZE. C THIS ARRAY IS THE RESULT OF CALLING DUNI 100000 TIMES C WITH ISEED=305 AND K=96. DATA U/ *0.471960981577884755837789724978D+00, *0.930323453205669578433639632431D+00, *0.110161790933730836587127944899D+00, *0.571501996273139518362638757010D-01, *0.402467554779738266237538503137D+00, *0.451181953427459489458279456915D+00, *0.296076152342721102174129954053D+00, *0.128202189325888116466879622359D-01, *0.314274693850973603980853259266D+00, *0.335521366752294932468163594171D-02, *0.488685045200439371607850367840D+00, *0.195470426865656758693860613516D+00, *0.864162706791773556901599326053D+00, *0.335505955815259203596381170316D+00, *0.377190200199058085469526470541D+00, *0.400780392114818314671676525916D+00, *0.374224214182207466262750307281D+00/ DATA I,J,K,C/17,5,47,CSAVE/ C C BASIC GENERATOR IS FIBONACCI C DUNI = U(I)-U(J) IF(DUNI.LT.0.0D0)DUNI = DUNI+1.0D0 U(I) = DUNI I = I-1 IF(I.EQ.0)I = 17 J = J-1 IF(J.EQ.0)J = 17 C C SECOND GENERATOR IS CONGRUENTIAL C C = C-CD IF(C.LT.0.0D0) C=C+CM C C COMBINATION GENERATOR C DUNI = DUNI-C IF(DUNI.LT.0.0D0)DUNI = DUNI+1.0D0 RETURN C ENTRY DUSTAR(ISEED) C C SET UP ... C CONVERT ISEED TO FOUR SMALLISH POSITIVE INTEGERS. C I1 = MOD(ABS(ISEED),177)+1 J1 = MOD(ABS(ISEED),167)+1 K1 = MOD(ABS(ISEED),157)+1 L1 = MOD(ABS(ISEED),147)+1 C C GENERATE RANDOM BIT PATTERN IN ARRAY BASED ON GIVEN SEED. C DO 2 II = 1,17 S = 0.0D0 T = 0.5D0 C DO FOR EACH OF THE BITS OF MANTISSA OF WORD C LOOP OVER K BITS, WHERE K IS DEFAULTED TO 47 BUT CAN C BE CHANGED BY USER CALL TO DUNIB(K) DO 3 JJ = 1,K M1 = MOD(MOD(I1*J1,179)*K1,179) I1 = J1 J1 = K1 K1 = M1 L1 = MOD(53*L1+1,169) IF(MOD(L1*M1,64).GE.32)S=S+T 3 T = 0.5D0*T 2 U(II) = S DUSTAR = FLOAT(ISEED) RETURN C ENTRY DUNIB(KK) IF(KK.LE.47)THEN K=47 ELSE K=KK ENDIF DUNIB=FLOAT(K) END DOUBLE PRECISION FUNCTION DRNOR() C***BEGIN PROLOGUE DRNOR C***DATE WRITTEN 810915 (YYMMDD) C***REVISION DATE 870419 (YYMMDD) C***CATEGORY NO. L6A14 C***KEYWORDS RANDOM NUMBERS, NORMAL DEVIATES C***AUTHOR KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS C MARSAGLIA, GEORGE, SUPERCOMPUTER RES. INST., FLORIDA ST. U. C C***PURPOSE GENERATES NORMAL RANDOM NUMBERS, WITH MEAN ZERO AND C UNIT STANDARD DEVIATION, OFTEN DENOTED N(0,1). C C***DESCRIPTION C C DRNOR GENERATES NORMAL RANDOM NUMBERS WITH ZERO MEAN AND C UNIT STANDARD DEVIATION, OFTEN DENOTED N(0,1). C FROM THE BOOK, "NUMERICAL METHODS AND SOFTWARE" BY C D. KAHANER, C. MOLER, S. NASH C PRENTICE HALL, 1988 C USE C FIRST TIME.... C Z = DSTART(ISEED) C HERE ISEED IS ANY N O N - Z E R O INTEGER. C THIS CAUSES INITIALIZATION OF THE PROGRAM. C DSTART RETURNS A DOUBLE PRECISION ECHO OF ISEED. C C SUBSEQUENT TIMES... C Z = DRNOR() C CAUSES THE NEXT DOUBLE PRECISION RANDOM NUMBER C TO BE RETURNED AS Z. C C..................................................................... C TYPICAL USAGE C C DOUBLE PRECISION DSTART,DRNOR,Z C INTEGER ISEED,I C ISEED = 305 C Z = DSTART(ISEED) C DO 1 I = 1,10 C Z = DRNOR() C WRITE(*,'(1X,D20.15)') Z C 1 CONTINUE C END C C C***REFERENCES MARSAGLIA & TSANG, "A FAST, EASILY IMPLEMENTED C METHOD FOR SAMPLING FROM DECREASING OR C SYMMETRIC UNIMODAL DENSITY FUNCTIONS", C PUBLISHED IN SIAM J SISC, JUNE 1984. C***ROUTINES CALLED (NONE) C***END PROLOGUE DRNOR DOUBLE PRECISION AA,B,C,C1,C2,PC,X,Y,XN,V(65),DSTART,U(17), *S,T,UN,VNI INTEGER J,IA,IB,IC,II,JJ,ID,III,JJJ,L SAVE U,II,JJ C DATA AA,B,C/0.123758602991705622657860318807D+02 +,0.487899177760378968003825536710D+00 +,0.126770580788654778410032042685D+02/ DATA C1,C2,PC,XN/0.9689279D0,1.301198D0 +,0.195830333955546914251231662871D-01 +,0.277699426966287548981739308903D+01/ DATA (V(L), L=1,18) +/0.340945028703977890881422029377D+00 +,0.457314591866933536170242641762D+00 +,0.539779281611666939596931167481D+00 +,0.606242679653048904661174385559D+00 +,0.663169062764524854743428299352D+00 +,0.713697459056025891183276150202D+00 +,0.759612474933920604605610034675D+00 +,0.802035600355531312751497342081D+00 +,0.841722667978955421276418428136D+00 +,0.879210223208313639290346470191D+00 +,0.914894804386750587541168254518D+00 +,0.949079113753090250631877133376D+00 +,0.982000481239888200218207508382D+00 +,0.101384923802994173461911276018D+01 +,0.104478103674017351825847605485D+01 +,0.107492538202855351339149779813D+01 +,0.110439170226812581109973656162D+01 +,0.113327377624394079212251428682D+01/ DATA (V(L), L=19,37) +/0.116165303013393172842858957666D+01 +,0.118960104083873798956793871425D+01 +,0.121718147070087121538258873613D+01 +,0.124445158789824683238161436879D+01 +,0.127146348057211969375402099579D+01 +,0.129826504188319751190626947962D+01 +,0.132490078218086109537654808436D+01 +,0.135141250993337129690631764473D+01 +,0.137783991287001181384096757263D+01 +,0.140422106355997540689484486002D+01 +,0.143059286850269131403410180874D+01 +,0.145699147613767157824869156623D+01 +,0.148345265660321931119703498108D+01 +,0.151001216431851991531882612256D+01 +,0.153670609335952099134607533122D+01 +,0.156357123503769104042967185962D+01 +,0.159064544701425352365935513885D+01 +,0.161796804367444698360816323707D+01 +,0.164558021836908161542865488149D+01/ DATA (V(L), L=38,55) +/0.167352550956703867146009214486D+01 +,0.170185032506274055254533570699D+01 +,0.173060454131778319060975251429D+01 +,0.175984219903830120010946138955D+01 +,0.178962232156657450014351836107D+01 +,0.182000989013069176863519209140D+01 +,0.185107702023027589942628767312D+01 +,0.188290439759287281399927405628D+01 +,0.191558305194303202395065401364D+01 +,0.194921657491636060191700129156D+01 +,0.198392392890568577258506733664D+01 +,0.201984305290623555306662745946D+01 +,0.205713555999009616804474181513D+01 +,0.209599295624939161758467989658D+01 +,0.213664502254438986524966622832D+01 +,0.217937134039813565892460111431D+01 +,0.222451750721601784110056845259D+01 +,0.227251855485014779874266158018D+01/ DATA (V(L), L=56,65) +/0.232393382009430256940425938218D+01 +,0.237950077408282829688673722776D+01 +,0.244022179797994340264326423618D+01 +,0.250751170186531701106382130475D+01 +,0.258346583522542956831304962942D+01 +,0.267139159032083601869533973173D+01 +,0.277699426966286466722522163764D+01 +,0.277699426966286466722522163764D+01 +,0.277699426966286466722522163764D+01 +,0.277699426966286466722522163764D+01/ C LOAD DATA ARRAY IN CASE USER FORGETS TO INITIALIZE. C THIS ARRAY IS THE RESULT OF CALLING DUNI 100000 TIMES C WITH SEED 305. DATA U/ *0.47196 09815 77884 75583 77897 24978D+00, *0.93032 34532 05669 57843 36396 32431D+00, *0.11016 17909 33730 83658 71279 44899D+00, *0.57150 19962 73139 51836 26387 57010D-01, *0.40246 75547 79738 26623 75385 03137D+00, *0.45118 19534 27459 48945 82794 56915D+00, *0.29607 61523 42721 10217 41299 54053D+00, *0.12820 21893 25888 11646 68796 22359D-01, *0.31427 46938 50973 60398 08532 59266D+00, *0.33552 13667 52294 93246 81635 94171D-02, *0.48868 50452 00439 37160 78503 67840D+00, *0.19547 04268 65656 75869 38606 13516D+00, *0.86416 27067 91773 55690 15993 26053D+00, *0.33550 59558 15259 20359 63811 70316D+00, *0.37719 02001 99058 08546 95264 70541D+00, *0.40078 03921 14818 31467 16765 25916D+00, *0.37422 42141 82207 46626 27503 07281D+00/ C DATA II,JJ / 17, 5 / C C***FIRST EXECUTABLE STATEMENT DRNOR C C FAST PART... C C C BASIC GENERATOR IS FIBONACCI C UN = U(II)-U(JJ) IF(UN.LT.0.0D0) UN = UN+1. U(II) = UN C U(II) AND UN ARE UNIFORM ON [0,1) C VNI IS UNIFORM ON [-1,1) VNI = UN + UN -1. II = II-1 IF(II.EQ.0)II = 17 JJ = JJ-1 IF(JJ.EQ.0)JJ = 17 C INT(UN(II)*128) IN RANGE [0,127], J IS IN RANGE [1,64] J = MOD(INT(U(II)*128),64)+1 C PICK SIGN AS VNI IS POSITIVE OR NEGATIVE DRNOR = VNI*V(J+1) IF(ABS(DRNOR).LE.V(J))RETURN C C SLOW PART; AA IS A*F(0) X = (ABS(DRNOR)-V(J))/(V(J+1)-V(J)) C Y IS UNIFORM ON [0,1) Y = U(II)-U(JJ) IF(Y.LT.0.0D0) Y = Y+1. U(II) = Y II = II-1 IF(II.EQ.0)II = 17 JJ = JJ-1 IF(JJ.EQ.0)JJ = 17 C S = X+Y IF(S.GT.C2)GO TO 11 IF(S.LE.C1)RETURN IF(Y.GT.C-AA*EXP(-.5D0*(B-B*X)**2))GO TO 11 IF(EXP(-.5D0*V(J+1)**2)+Y*PC/V(J+1).LE.EXP(-.5D0*DRNOR**2)) *RETURN C C TAIL PART; .36010157... IS 1.0D0/XN C Y IS UNIFORM ON [0,1) 22 Y = U(II)-U(JJ) IF(Y.LE.0.0D0) Y = Y+1. U(II) = Y II = II-1 IF(II.EQ.0)II = 17 JJ = JJ-1 IF(JJ.EQ.0)JJ = 17 C X = 0.360101571301190680192994239651D+00*LOG(Y) C Y IS UNIFORM ON [0,1) Y = U(II)-U(JJ) IF(Y.LE.0.0D0) Y = Y+1. U(II) = Y II = II-1 IF(II.EQ.0)II = 17 JJ = JJ-1 IF(JJ.EQ.0)JJ = 17 IF( -2.0D0*LOG(Y).LE.X**2 )GO TO 22 DRNOR = SIGN(XN-X,DRNOR) RETURN 11 DRNOR = SIGN(B-B*X,DRNOR) RETURN C C C FILL ENTRY DSTART(ISEED) IF(ISEED.NE.0) THEN C C SET UP ... C GENERATE RANDOM BIT PATTERN IN ARRAY BASED ON GIVEN SEED C II = 17 JJ = 5 IA = MOD(ABS(ISEED),32707) IB = 1111 IC = 1947 DO 2 III = 1,17 S = 0.0D0 T = 0.5D0 C DO FOR EACH OF THE BITS OF MANTISSA OF WORD C LOOP OVER 95 BITS, ENOUGH FOR MOST MACHINES C IN DOUBLE PRECISION. DO 3 JJJ = 1,95 ID = IC-IA IF(ID.GE.0)GOTO 4 ID = ID+32707 S = S+T 4 IA = IB IB = IC IC = ID 3 T = 0.5D0*T 2 U(III) = S ENDIF C RETURN FLOATING ECHO OF ISEED DSTART=ISEED RETURN END