!****************************frng************************** B Bunk 8/2010 
!                                                           rev    5/2014
!  Multiplicative lagged Fibonacci random number generator
!
!  Implementation with 64 bits and recursion
!     x(n) = x(n-13) * x(n-31) mod 2^64
!
!  Entries: 
!  - frng()    real
!  - dfrng()   double precision
!  - ifrng()   integer(4)  (31 most significant bits)
!  - lfrng()   integer(8)  (63 most significant bits)
!  - calls which fill a random vector:
!     frngv(n, rvec), dfrngv(n, dvec), ifrngv(n, ivec), lfrngv(n, lvec)
!  - frnget(iseed), frnset(iseed) : get/set register
!           integer(8) iseed(31)    (63 most significant bits)
!  - frnini(num) : initialise register (num = 0 .. 714024)
!
!  There are no unsigned integers in Fortran90.
!  64bit integers with the most significant bit=1 appear to be negative.
!  Anyway - multiplication mod 2^64 works as if they were unsigned ints.

module frngmod

  public frng, dfrng, ifrng, lfrng
  public frngv, dfrngv, ifrngv, lfrngv
  public frnget, frnset, frnini

  private

  integer, parameter :: k4=selected_int_kind(9)
  integer, parameter :: k8=selected_int_kind(18)

  integer(k8) :: ireg(0:30) = &
    (/ 7569555838333416573_k8, 3279508493814418161_k8,-3845546956793070485_k8 &
     ,-7089769803070399727_k8,-6913536167697213337_k8,   48396622165758097_k8 &
     ,-8647237902909285879_k8,-7950859126646081777_k8,-2684521842256023033_k8 &
     ,  592708693962015337_k8,  249471268379597569_k8, -153673632038031339_k8 &
     , 1095135063400818497_k8,-6172267299408869197_k8,-5490224961767833701_k8 &
     , 6258116919374872021_k8, 7125927167643344901_k8, 3558223830886234463_k8 &
     , 6629248264238368005_k8, 9006395648747532825_k8,-8631203895171811737_k8 &
     , 7397323557512132751_k8,-1521152214184682699_k8,-7422728869875354831_k8 &
     ,-7439117882532651161_k8, 7857178218845272995_k8,  406276102641255367_k8 &
     , 6196339901492533131_k8, 7804281187626479095_k8,-2907244538723958819_k8 &
     , 4020080149601329621_k8/)
  integer :: i0 = 0, i1 = 18
  integer :: iup(0:30) = (/(i, i=1,30), 0/)

  double precision, parameter :: &
       dnorm63 = 1.08420217248550443400745280086994171142578125d-19 ! 2^(-63)
  real, parameter :: rnorm63 = dnorm63
  real, parameter :: rnorm31 = 4.656612873077392578125e-10          ! 2^(-31)

contains

real function frng()
  ireg(i0) = ireg(i0) * ireg(i1)
  !!!frng = rnorm31 * ishft(ireg(i0), -33)
  frng = rnorm63 * ishft(ireg(i0), -1)
  i0 = iup(i0)
  i1 = iup(i1)
end function frng

double precision function dfrng()
  ireg(i0) = ireg(i0) * ireg(i1)
  dfrng = dnorm63 * ishft(ireg(i0), -1)
  i0 = iup(i0)
  i1 = iup(i1)
end function dfrng

function ifrng()
  integer(k4) ifrng           ! k4 is defined in the module...
  ireg(i0) = ireg(i0) * ireg(i1)
  ifrng = ishft(ireg(i0), -33)
  i0 = iup(i0)
  i1 = iup(i1)
end function ifrng

function lfrng()
  integer(k8) lfrng           ! k8 is defined in the module...
  ireg(i0) = ireg(i0) * ireg(i1)
  lfrng = ishft(ireg(i0), -1)
  i0 = iup(i0)
  i1 = iup(i1)
end function lfrng

subroutine frngv(n, rvec)
  real rvec(n)
  do i=1,n
    ireg(i0) = ireg(i0) * ireg(i1)
    !!!rvec(i) = rnorm31 * ishft(ireg(i0), -33)
    rvec(i) = rnorm63 * ishft(ireg(i0), -1)
    i0 = iup(i0)
    i1 = iup(i1)
  enddo
end subroutine frngv

subroutine dfrngv(n, dvec)
  double precision dvec(n)
  do i=1,n
    ireg(i0) = ireg(i0) * ireg(i1)
    dvec(i) = dnorm63 * ishft(ireg(i0), -1)
    i0 = iup(i0)
    i1 = iup(i1)
  enddo
end subroutine dfrngv

subroutine ifrngv(n, ivec)
  integer(k4) ivec(n)
  do i=1,n
    ireg(i0) = ireg(i0) * ireg(i1)
    ivec(i) = ishft(ireg(i0), -33)
    i0 = iup(i0)
    i1 = iup(i1)
  enddo
end subroutine ifrngv

subroutine lfrngv(n, lvec)
  integer(k8) lvec(n)
  do i=1,n
    ireg(i0) = ireg(i0) * ireg(i1)
    lvec(i) = ishft(ireg(i0), -1)
    i0 = iup(i0)
    i1 = iup(i1)
  enddo
end subroutine lfrngv

subroutine frnget(iseed)
  integer(k8) iseed(0:30)
  ipt = i0
  do i=0,30
    iseed(i)   = ishft(ireg(ipt), -1)
    ipt = iup(ipt)
  enddo
end subroutine frnget

subroutine frnset(iseed)
  integer(k8) iseed(0:30)
  do i=0,30
    ireg(i) = 2*iseed(i) + 1
  enddo
  i0 = 0
  i1 = 18
end subroutine frnset

subroutine frnini(num)
  integer, parameter ::  mr=714025, ia=1366, ic=150889, mrhalf=(mr+1)/2

  irand = mod(abs(num), mr)
  do i=0,30
    ireg(i) = 0
    do ibit=0,62
      irand = mod(ia*irand+ic, mr)
      ireg(i) = 2*ireg(i) + irand/mrhalf
    enddo
    ireg(i) = 2*ireg(i) + 1               ! all odd!
  enddo

  i0 = 0
  i1 = 18
end subroutine frnini

end module frngmod
