C----------------------------------------------------------------------- c Last date of change 25.08.03 V.Uzhinsky SUBROUTINE RNDMGET1(DSEED) common/rndu/DS DOUBLE PRECISION DS(2), DSEED, DX24 DATA DX24 / 1677 7216.D0 / SAVE DX24 DSEED = DS(2)*DX24 + DS(1) RETURN END SUBROUTINE RNDMSET1(DSEED) common/rndu/DS DOUBLE PRECISION DS(2), DSEED, DX24 DATA DX24 / 1677 7216.D0 / SAVE DX24 DS(2) = DINT(DSEED/DX24) DS(1) = DSEED - DS(2)*DX24 CALL sgrnd(IDINT(DSEED)) CALL init_genrand(IDINT(DSEED)) RETURN END FUNCTION RNDM1(IX) c1 RNDM1 = RANF1(1) 1 RNDM1 = genrand_res53() c print *, RNDM1 if(RNDM1.eq.0.) go to 1 RETURN END C----------------------------------------------------------------------- FUNCTION RANF1(IX) C UNIFORM RANDOM NUMBER GENERATOR FROM CERN LIBRARY C----------------------------------------------------------------------- DOUBLE PRECISION DS(2), DM(2) DOUBLE PRECISION DX24, DX48 DOUBLE PRECISION DL, DC, DU, DR common/rndu/DS DATA DM / 1518 4245.D0, 265 1554.D0 / DATA DX24 / 1677 7216.D0 / DATA DX48 / 281 4749 7671 0656.D0 / SAVE DM, DX24, DX48 10 DL = DS(1) * DM(1) DC = DINT(DL/DX24) DL = DL - DC*DX24 DU = DS(1)*DM(2) + DS(2)*DM(1) + DC DS(2) = DU - DINT(DU/DX24)*DX24 DS(1) = DL DR = (DS(2)*DX24 + DS(1)) / DX48 RANF1 = SNGL(DR) * RANF1 = DR RETURN END ************************************************************************ subroutine sgrnd(seed) * implicit integer(a-z) * * Period parameters parameter(N = 624) * dimension mt(0:N-1) * the array for the state vector common /block/mti,mt save /block/ * * setting initial seeds to mt[N] using * the generator Line 25 of Table 1 in * [KNUTH 1981, The Art of Computer Programming * Vol. 2 (2nd Ed.), pp102] * mt(0)= iand(seed,-1) do 1000 mti=1,N-1 mt(mti) = iand(69069 * mt(mti-1),-1) 1000 continue * return end ************************************************************************ double precision function grnd() * implicit integer(a-z) * * Period parameters parameter(N = 624) parameter(N1 = N+1) parameter(M = 397) parameter(MATA = -1727483681) * constant vector a parameter(UMASK = -2147483648) * most significant w-r bits parameter(LMASK = 2147483647) * least significant r bits * Tempering parameters parameter(TMASKB= -1658038656) parameter(TMASKC= -272236544) * dimension mt(0:N-1) * the array for the state vector common /block/mti,mt save /block/ data mti/N1/ * mti==N+1 means mt[N] is not initialized * dimension mag01(0:1) data mag01/0, MATA/ save mag01 * mag01(x) = x * MATA for x=0,1 * TSHFTU(y)=ishft(y,-11) TSHFTS(y)=ishft(y,7) TSHFTT(y)=ishft(y,15) TSHFTL(y)=ishft(y,-18) * if(mti.ge.N) then * generate N words at one time if(mti.eq.N+1) then * if sgrnd() has not been called, call sgrnd(4357) * a default initial seed is used endif * do 1000 kk=0,N-M-1 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1))) 1000 continue do 1100 kk=N-M,N-2 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1))) 1100 continue y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK)) mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1))) mti = 0 endif * y=mt(mti) mti=mti+1 y=ieor(y,TSHFTU(y)) y=ieor(y,iand(TSHFTS(y),TMASKB)) y=ieor(y,iand(TSHFTT(y),TMASKC)) y=ieor(y,TSHFTL(y)) * if(y.lt.0) then grnd=(dble(y)+2.0d0**32)/(2.0d0**32-1.0d0) else grnd=dble(y)/(2.0d0**32-1.0d0) endif * return end c----------------------------------------------------------------------- c initialize mt(0:N-1) with a seed c----------------------------------------------------------------------- subroutine init_genrand(s) integer s integer N integer DONE integer ALLBIT_MASK parameter (N=624) parameter (DONE=123456789) integer mti,initialized integer mt(0:N-1) common /mt_state1/ mti,initialized common /mt_state2/ mt common /mt_mask1/ ALLBIT_MASK c call mt_initln mt(0)=iand(s,ALLBIT_MASK) do 100 mti=1,N-1 mt(mti)=1812433253* & ieor(mt(mti-1),ishft(mt(mti-1),-30))+mti mt(mti)=iand(mt(mti),ALLBIT_MASK) 100 continue initialized=DONE c return end c----------------------------------------------------------------------- c initialize by an array with array-length c init_key is the array for initializing keys c key_length is its length c----------------------------------------------------------------------- subroutine init_by_array(init_key,key_length) integer init_key(0:*) integer key_length integer N integer ALLBIT_MASK integer TOPBIT_MASK parameter (N=624) integer i,j,k integer mt(0:N-1) common /mt_state2/ mt common /mt_mask1/ ALLBIT_MASK common /mt_mask2/ TOPBIT_MASK c call init_genrand(19650218) i=1 j=0 do 100 k=max(N,key_length),1,-1 mt(i)=ieor(mt(i),ieor(mt(i-1),ishft(mt(i-1),-30))*1664525) & +init_key(j)+j mt(i)=iand(mt(i),ALLBIT_MASK) i=i+1 j=j+1 if(i.ge.N)then mt(0)=mt(N-1) i=1 endif if(j.ge.key_length)then j=0 endif 100 continue do 200 k=N-1,1,-1 mt(i)=ieor(mt(i),ieor(mt(i-1),ishft(mt(i-1),-30))*1566083941)-i mt(i)=iand(mt(i),ALLBIT_MASK) i=i+1 if(i.ge.N)then mt(0)=mt(N-1) i=1 endif 200 continue mt(0)=TOPBIT_MASK c return end c----------------------------------------------------------------------- c generates a random number on [0,0xffffffff]-interval c----------------------------------------------------------------------- function genrand_int32() integer genrand_int32 integer N,M integer DONE integer UPPER_MASK,LOWER_MASK,MATRIX_A integer T1_MASK,T2_MASK parameter (N=624) parameter (M=397) parameter (DONE=123456789) integer mti,initialized integer mt(0:N-1) integer y,kk integer mag01(0:1) common /mt_state1/ mti,initialized common /mt_state2/ mt common /mt_mask3/ UPPER_MASK,LOWER_MASK,MATRIX_A,T1_MASK,T2_MASK common /mt_mag01/ mag01 c if(initialized.ne.DONE)then call init_genrand(21641) endif c if(mti.ge.N)then do 100 kk=0,N-M-1 y=ior(iand(mt(kk),UPPER_MASK),iand(mt(kk+1),LOWER_MASK)) mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1))) 100 continue do 200 kk=N-M,N-1-1 y=ior(iand(mt(kk),UPPER_MASK),iand(mt(kk+1),LOWER_MASK)) mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1))) 200 continue y=ior(iand(mt(N-1),UPPER_MASK),iand(mt(0),LOWER_MASK)) mt(kk)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1))) mti=0 endif c y=mt(mti) mti=mti+1 c y=ieor(y,ishft(y,-11)) y=ieor(y,iand(ishft(y,7),T1_MASK)) y=ieor(y,iand(ishft(y,15),T2_MASK)) y=ieor(y,ishft(y,-18)) c genrand_int32=y return end c----------------------------------------------------------------------- c generates a random number on [0,0x7fffffff]-interval c----------------------------------------------------------------------- function genrand_int31() integer genrand_int31 integer genrand_int32 genrand_int31=int(ishft(genrand_int32(),-1)) return end c----------------------------------------------------------------------- c generates a random number on [0,1]-real-interval c----------------------------------------------------------------------- function genrand_real1() double precision genrand_real1,r integer genrand_int32 r=dble(genrand_int32()) if(r.lt.0.d0)r=r+2.d0**32 genrand_real1=r/4294967295.d0 return end c----------------------------------------------------------------------- c generates a random number on [0,1)-real-interval c----------------------------------------------------------------------- function genrand_real2() double precision genrand_real2,r integer genrand_int32 r=dble(genrand_int32()) if(r.lt.0.d0)r=r+2.d0**32 genrand_real2=r/4294967296.d0 return end c----------------------------------------------------------------------- c generates a random number on (0,1)-real-interval c----------------------------------------------------------------------- function genrand_real3() double precision genrand_real3,r integer genrand_int32 r=dble(genrand_int32()) if(r.lt.0.d0)r=r+2.d0**32 genrand_real3=(r+0.5d0)/4294967296.d0 return end c----------------------------------------------------------------------- c generates a random number on [0,1) with 53-bit resolution c----------------------------------------------------------------------- function genrand_res53() double precision genrand_res53 integer genrand_int32 double precision a,b a=dble(ishft(genrand_int32(),-5)) b=dble(ishft(genrand_int32(),-6)) if(a.lt.0.d0)a=a+2.d0**32 if(b.lt.0.d0)b=b+2.d0**32 genrand_res53=(a*67108864.d0+b)/9007199254740992.d0 return end c----------------------------------------------------------------------- c initialize large number (over 32-bit constant number) c----------------------------------------------------------------------- subroutine mt_initln integer ALLBIT_MASK integer TOPBIT_MASK integer UPPER_MASK,LOWER_MASK,MATRIX_A,T1_MASK,T2_MASK integer mag01(0:1) common /mt_mask1/ ALLBIT_MASK common /mt_mask2/ TOPBIT_MASK common /mt_mask3/ UPPER_MASK,LOWER_MASK,MATRIX_A,T1_MASK,T2_MASK common /mt_mag01/ mag01 CC TOPBIT_MASK = Z'80000000' CC ALLBIT_MASK = Z'ffffffff' CC UPPER_MASK = Z'80000000' CC LOWER_MASK = Z'7fffffff' CC MATRIX_A = Z'9908b0df' CC T1_MASK = Z'9d2c5680' CC T2_MASK = Z'efc60000' TOPBIT_MASK=1073741824 TOPBIT_MASK=ishft(TOPBIT_MASK,1) ALLBIT_MASK=2147483647 ALLBIT_MASK=ior(ALLBIT_MASK,TOPBIT_MASK) UPPER_MASK=TOPBIT_MASK LOWER_MASK=2147483647 MATRIX_A=419999967 MATRIX_A=ior(MATRIX_A,TOPBIT_MASK) T1_MASK=489444992 T1_MASK=ior(T1_MASK,TOPBIT_MASK) T2_MASK=1875247104 T2_MASK=ior(T2_MASK,TOPBIT_MASK) mag01(0)=0 mag01(1)=MATRIX_A return end