CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Random number generator ILFSR113  - integer 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 ilfsr113() gives one random integer in the 
C interval (1,2147483647).
C
C Before using ilfsr113 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  10 Feb 2000    v0.91   changed lfsrinit to allow for CONSTANT arguments
C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


      FUNCTION ilfsr113 ()
      integer*4 b,z1,z2,z3,z4,ilfsr113
      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)

      ilfsr113=ishft( ieor(ieor(ieor(z1,z2),z3),z4) , -1)

      return
      end


      SUBROUTINE lfsrinit(iinit)
      integer*4 idum,ia,im,iq,ir,ilfsr113,iinit
      integer*4 k,z1,z2,z3,z4,c1,c2,c3,c4
      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

      idum=ilfsr113()
      
      return
      end