!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ! ! Random number generator RLFSR113 - real version ! ! Following a suggestion of Pierre L'Ecuyer 1997 ! "Tables of maximally equidistributed combined LFSR generators" ! see http://www.iro.umontreal.ca/~lecuyer/papers.html ! ! A call to rlfsr113() gives one random real in the open ! interval (0,1). ! ! Before using rlfsr113 call lfsrinit(seed) to initialize ! the generator by random integers produced by Park/Millers ! minimal standard LCG. ! Seed should be any positive integer. ! ! ! FORTRAN version by Thomas Vojta, vojta@physik.tu-chemnitz.de ! ! History: ! 04 Feb 1998 v0.9 first FORTRAN implementation ! 05 Feb 1998 v0.91 corrected multiplicator am to 1/(2^31) ! 15 Apr 1999 v0.92 added real*8 rlfs113 in lfsrinit ! 10 Feb 2000 v0.93 changed lfsrinit to allow for CONSTANT arguments ! !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC FUNCTION rlfsr113 () implicit none integer,parameter :: i4b= kind(2147483647) integer,parameter :: r8b= kind(1.D200) real(r8b),parameter :: am=4.656612873077d-10 real(r8b) rlfsr113 integer(i4b) b,z1,z2,z3,z4 common /lfsrcom/z1,z2,z3,z4 b = ishft(ieor(ishft(z1,6),z1),-13) z1 = ieor(ishft(iand(z1,-2),18),b) b = ishft(ieor(ishft(z2,2),z2),-27) z2 = ieor(ishft(iand(z2,-8),2),b) b = ishft(ieor(ishft(z3,13),z3),-21) z3 = ieor(ishft(iand(z3,-16),7),b) b = ishft(ieor(ishft(z4,3),z4),-12) z4 = ieor(ishft(iand(z4,-128),13),b) rlfsr113=ishft( ieor(ieor(ieor(z1,z2),z3),z4) , -1)*am return end SUBROUTINE lfsrinit(iinit) implicit none integer,parameter :: r8b= kind(1.D200) integer,parameter :: i4b= kind(2147483647) integer(i4b) idum,ia,im,iq,ir,iinit integer(i4b) k,z1,z2,z3,z4,c1,c2,c3,c4 real(r8b) rlfsr113,rdum parameter (ia=16807,im=2147483647,iq=127773,ir=2836) common /lfsrcom/z1,z2,z3,z4 data c1 /B'11111111111111111111111111111110'/ data c2 /B'11111111111111111111111111111000'/ data c3 /B'11111111111111111111111111110000'/ data c4 /B'11111111111111111111111110000000'/ if ((c1.ne.-2).or.(c2.ne.-8).or.(c3.ne.-16).or.(c4.ne.-128)) then print *,'Nonstandard integer representation. Stoped.' stop endif idum=iinit if (idum.le.0) idum=1 k=(idum)/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum = idum + IM if (idum.lt.2) then z1=idum+2 else z1=idum endif k=(idum)/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum = idum + IM if (idum.lt.8) then z2=idum+8 else z2=idum endif k=(idum)/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum = idum + IM if (idum.lt.16) then z3=idum+16 else z3=idum endif k=(idum)/IQ idum=IA*(idum-k*IQ)-IR*k if (idum.lt.0) idum = idum + IM if (idum.lt.128) then z4=idum+128 else z4=idum endif rdum=rlfsr113() return end