! ranlxd.f90
!
! Copyright 2008 Claude Dion
! Copyright 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
!
! This program is free software; you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation; either version 3 of the License, or (at
! your option) any later version.
! 
! This program is distributed in the hope that it will be useful, but
! WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
! General Public License for more details.
! 
! You should have received a copy of the GNU General Public License
! along with this program; if not, write to the Free Software
! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

!************************************************************************
!
! Fortran 90 adaptation of the ranlxd2 random number generator of the
! GNU Scientific Library, which is an implementation of Martin Luescher's 
! second generation double-precision (48-bit) version of the RANLUX 
! generator. 
!
! Only the highest level of luxury, p = 397, is available here.
! See the GSL documentation for more information and references.
!
! C. M. Dion, Umea University, July 2008
!
!************************************************************************

module ranlxd2
  integer, parameter :: dp = kind(1.d0), long = kind(z'FFFFFFFF'), luxury = 397
  integer, dimension(12), parameter :: next = (/ 2,3,4,5,6,7,8,9,10,11,12,1 /)
  real(dp), parameter :: one_bit = 1.d0 / 281474976710656.d0 ! 1/2^48 

  integer :: ir, jr, ir_old
  real(dp) :: carry
  real(dp), dimension(12) :: xdbl 
  logical :: initialize = .true.
end module ranlxd2


#define RANLUX_STEP(x1,x2,i1,i2,i3) \
x1 = xdbl(i1) - xdbl(i2); \
if (x2 < 0.d0) then; \
   x1 = x1 - one_bit; \
   x2 = x2 + 1.d0; \
end if; \
xdbl(i3) = x2


subroutine increment_state()
  !
  ! Used internally to skip over a series of values ("luxury"),
  ! reducing correlations
  !
  use ranlxd2
  implicit none

  integer :: k, kk, kmax
  real(dp) :: y1, y2, y3

  k = 0
  do
     if (ir <= 1) exit

     y1 = xdbl(jr) - xdbl(ir)
     y2 = y1 - carry
     if (y2 < 0.d0) then
        carry = one_bit
        y2 = y2 + 1
     else
        carry = 0.d0
     end if
     xdbl(ir) = y2
     ir = next(ir)
     jr = next(jr)
     
     k = k+1
  end do
  
  kmax = luxury - 12;
  
  do kk = k,kmax,12
     y1 = xdbl(8) - xdbl(1) - carry
     
     RANLUX_STEP (y2, y1, 9, 2, 1)
     RANLUX_STEP (y3, y2, 10, 3, 2)
     RANLUX_STEP (y1, y3, 11, 4, 3)
     RANLUX_STEP (y2, y1, 12, 5, 4)
     RANLUX_STEP (y3, y2, 1, 6, 5)
     RANLUX_STEP (y1, y3, 2, 7, 6)
     RANLUX_STEP (y2, y1, 3, 8, 7)
     RANLUX_STEP (y3, y2, 4, 9, 8)
     RANLUX_STEP (y1, y3, 5, 10, 9)
     RANLUX_STEP (y2, y1, 6, 11, 10)
     RANLUX_STEP (y3, y2, 7, 12, 11)
     
     if (y3 < 0.d0) then
        carry = one_bit
        y3 = y3 + 1
     else
        carry = 0.d0
     end if
     xdbl(12) = y3
  end do

  kmax = luxury

  do k = kk,kmax-1
     y1 = xdbl(jr) - xdbl(ir)
     y2 = y1 - carry
     if (y2 < 0.d0) then
        carry = one_bit
        y2 = y2 + 1
     else
        carry = 0.d0
     end if
      xdbl(ir) = y2
      ir = next(ir)
      jr = next(jr)
   end do

  ir_old = ir

  return
end subroutine increment_state


real(kind(1.d0)) function ranlxd()
  !
  ! Generate a double precision floating point pseudo-random number
  ! 
  ! Subroutine ranlxd_set must be called once before any call to ranlxd
  ! to set up the random number generator
  !
  use ranlxd2
  implicit none

  if(initialize) call ranlxd_set(0)

  ir = next(ir)

  if (ir == ir_old) call increment_state

  ranlxd = xdbl(ir)

  return 
end function ranlxd


subroutine ranlxd_set(s)
  !
  ! Initialize the random number generator ranlxd
  !
  ! On input, s contains the seed of the generator
  ! While the seed is internally an unsigned 32 bit number, requiring thus 
  ! more than 32 bits to represent it as a Fortran signed integer, for 
  ! convenience a "regular" 32 bit integer is used to initialize the seed
  !
  use ranlxd2
  implicit none

  integer(long), parameter :: two32 = z'FFFFFFFF' ! 2^32

  integer, intent(in) :: s

  integer :: ibit, jbit, k, l
  integer(long) :: i, seed
  integer, dimension(31) :: xbit
  real(dp) :: x, y

  if (s <= 0) then
     seed = 1 ! default seed is 1 
  else
     seed = s
  end if

  i = iand(seed, two32)

  do k = 1,31
     xbit(k) = mod(i, 2_long)
     i = i/2_long
  end do

  ibit = 1
  jbit = 19

  do k = 1,12
     x = 0.d0
     
     do l = 1,48
        y = real(mod(xbit(ibit) + 1, 2),dp)
        x = 2.d0*x + y
        xbit(ibit) = mod(xbit(ibit) + xbit(jbit), 2)

        ibit = ibit + 1
        if(ibit == 32) ibit = 1
        jbit = jbit + 1
        if(jbit == 32) jbit = 1
     end do

     xdbl(k) = one_bit * x

  end do

  carry = 0
  ir = 12
  jr = 8
  ir_old = 1
  initialize = .false.

  return
end subroutine ranlxd_set
