CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Random number generator RLFSR113  - real version
C
C Following a suggestion of Pierre L'Ecuyer 1997
C "Tables of maximally equidistributed combined LFSR generators"
C see http://www.iro.umontreal.ca/~lecuyer/papers.html
C
C A call to rlfsr113() gives one random real in the open
C interval (0,1).
C
C Before using rlfsr113 call lfsrinit(seed) to initialize
C the generator by random integers produced by Park/Millers
C minimal standard LCG.
C Seed should be any positive integer.
C 
C FORTRAN version by Thomas Vojta, vojta@physik.tu-chemnitz.de
C 
C History:
C  04 Feb 1998    v0.9    first FORTRAN implementation
C  05 Feb 1998    v0.91   corrected multiplicator am to 1/(2^31)
C  15 Apr 1999    v0.92   added real*8 rlfs113 in lfsrinit
C  10 Feb 2000    v0.93   changed lfsrinit to allow for CONSTANT arguments
C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


      FUNCTION rlfsr113 ()
      integer*4 b,z1,z2,z3,z4
      real*8 am,rlfsr113
      parameter (am=4.656612873077e-10)
      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)
      integer*4 idum,ia,im,iq,ir,iinit
      integer*4 k,z1,z2,z3,z4,c1,c2,c3,c4
      real*8    rlfsr113
      parameter (ia=16807,im=2147483647,iq=127773,ir=2836)
      common /lfsrcom/z1,z2,z3,z4

C Check whether the FORTRAN integers can be used as unsigned long !

      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

C Initialize z1,z2,z3,z4

      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

C Make a single call to rlfsr113() to achieve a valid state

      idum=rlfsr113()
      
      return
      end