*=======================================================================
* Random number generator R1279
* Usage
* 1. input an initial seed ISEED
* 2. warm up the generator: CALL  RINITIALIZE(ISEED)
* 3. generate a number with the function RANF()
*=======================================================================

      implicit none
      integer i,iseed
      real*8 ranf

      iseed = 136667
      call RINITIALIZE(ISEED)
      do i=1,10
         write(*,*)ranf()
      enddo
      stop
      end


      SUBROUTINE RINITIALIZE(ISEED)
      IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
      PARAMETER(NAB3=101280)
      COMMON/CNUM/NUM
      NUM=NAB3-1280
      CALL RANINI(ISEED)
      END


      FUNCTION RANF()
      IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)
      COMMON/CNUM/NUM
      ranf =  RANDA(NUM)
      NUM=NUM+1
      END

C***********************************************************************

      SUBROUTINE RANINI(ISeed)
      IMPLICIT NONE

      INTEGER NAB3
      PARAMETER(NAB3 = 101280)

      INTEGER ISeed,IMod,I,J,K
      INTEGER RanVec,IMax
      REAL*8 RMod,PMod,DMaxI,RanVec2
      COMMON/dom/RanVec(NAB3),Ranvec2(NAB3)

      IMax = 2147483647
      DMaxI = 1.0D0/2147483647.0D0
      RMod = DBLE(ISeed)
      PMod = DBLE(IMax)

      DO I = 1,1000
        RMod = RMod*16807.0D0
        IMod = RMod*DMaxI
        RMod = RMod - PMod*IMod
      END DO

      DO I = 1,1279
        RanVec(I) = 0
        DO J = 0,30
          DO K = 1,36
            RMod = RMod*16807.0D0
            IMod = RMod*DMaxi
            RMod = RMod - PMod*IMod
          END DO
          RMod = RMod*16807.0D0
          IMod = RMod*DMaxi
          RMod = RMod - PMod*IMod
          IF (RMod .GT. 0.5D0*PMod) RanVec(I) = IBSET(RanVec(I),J)
        END DO
      END DO

C**  Generate 1000 random numbers to warm up the generator
      CALL RANDOM(1000)

      RETURN
      END

C***********************************************************************

      SUBROUTINE RANDOM(Number)

      IMPLICIT NONE
      INTEGER NAB3
      PARAMETER(NAB3 = 101280)
      real*8 RanVec2
      INTEGER RanVec,Number,I
      COMMON/dom/RanVec(NAB3),RanVec2(NAB3)

C**  This works because Number will always be a multiple of four for this
C**  program
C**  Unroll this loop for extra speed

      DO I = 1,Number,4
        RanVec(I+1279) = IEOR(RanVec(I),    RanVec(I+216))
        RanVec(I+1280) = IEOR(RanVec(I+1),  RanVec(I+217))
        RanVec(I+1281) = IEOR(RanVec(I+2),  RanVec(I+218))
        RanVec(I+1282) = IEOR(RanVec(I+3),  RanVec(I+219))
      END DO

C**  Copy the final 1279 elements to the beginning for use on the next call
C**  Unroll this loop for extra speed

      DO I = 1,1276,4
        RanVec(I)   = RanVec(I+Number)
        RanVec(I+1)   = RanVec(I+1+Number)
        RanVec(I+2)   = RanVec(I+2+Number)
        RanVec(I+3)   = RanVec(I+3+Number)
      END DO

      RanVec(1277) = RanVec(1277 + Number)
      RanVec(1278) = RanVec(1278 + Number)
      RanVec(1279) = RanVec(1279 + Number)

      do I = 1,number
        RanVec2(I) = dble(RanVec(I))*4.656612875D-10
      end do

      RETURN
      END

      FUNCTION RANDA(NUM)
      INTEGER NUM,ranvec,NAB3
      REAL*8 RANVEC2,RANDA
      PARAMETER(NAB3=101280)
      COMMON/dom/RanVec(NAB3),Ranvec2(NAB3)

        if(NUM.ge.NAB3-1280)then
          CALL RANDOM(NAB3-1280)
          RANDA = ranvec2(1)
          NUM = 2
        else
          RANDA = ranvec2(NUM)
          NUM = NUM + 1
        endif
      end