!=======================================================================
! 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(16) :: ranf

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


      SUBROUTINE RINITIALIZE(ISEED)
      IMPLICIT REAL(16)(A-H,O-Z)
      IMPLICIT 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)
      IMPLICIT REAL(16)(A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      COMMON/CNUM/NUM
      ranf =  RANDA(NUM)
      NUM=NUM+1
      END

!***********************************************************************

      SUBROUTINE RANINI(ISeed)
      IMPLICIT NONE

      INTEGER :: NAB3
      PARAMETER(NAB3 = 101280)

      INTEGER :: ISeed,IMod,I,J,K
      INTEGER :: RanVec,IMax
      REAL(16) :: 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

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

      RETURN
      END

!***********************************************************************

      SUBROUTINE RANDOM(Number)

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

!**  This works because Number will always be a multiple of four for this
!**  program
!**  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

!**  Copy the final 1279 elements to the beginning for use on the next call
!**  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(16) :: 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