! this module contains variables and methods which deal with
! pseudo-random numbers

module randommod

 use decimal, only: dop

 implicit none

 integer,allocatable :: seed(:)
 integer :: seedsize

contains

 function randuniform()

!Wrapper function around a pseudo-random number generator which
!generates a uniform distribution on the interval [0,1]
!This wrapper allows for the seamless change of random generator across
!the code
!Current implementation uses the random generator provided by the
!fortran library. It is noted that different compilers implement
!different generators so different random numbers should be expected.
!A different sequence of random numbers will be generated each time the
!program is run unless the generator has been seeded specifying a seed on
!the input run-section. This can be important for debugging purposes.

  implicit none

  real(kind=8) :: randuniform

  call random_number(randuniform)

 end function randuniform


 function randgaussian()

!This function implements the Box-Muller algorithm that from 2 uniform
!random numbers generates 2 random variables with a normal/gaussian
!distribution centered in 0 and variance equal to 1.
!The pair of gaussian distributed random numbers is output as the real
!and imaginary parts of a complex number
!
  implicit none

  complex(dop)        :: randgaussian
  real(dop)           :: x1,x2 !2 intermediate uniform random numbers
  real(dop)           :: twopi
         
  twopi=8.d0*atan(1.d0) !2*pi constant
  x1=randuniform()
  x2=randuniform()
  randgaussian=dcmplx(dsqrt(-2.d0*log(x1))*cos(twopi*x2),&
                     dsqrt(-2.d0*log(x1))*sin(twopi*x2))
  
 end function randgaussian


 subroutine sobseq(n,x,ierr)

!This subroutine generates a Sobol sequence of low-discrepancy random numbers
!this is useful for monte carlo integration schemes, since this ensures
!smoother coverage of space
!This subroutine was previously in lib/utilities/numerical.f
!and perhaps comes from NR. Maybe we want a more modern implementation?

  implicit none

  integer,intent(in) :: n
  integer,intent(out) :: ierr !ierr=0 no error; ierr=1 maxbit too small
  integer,parameter :: MAXBIT=30,MAXsDIM=6
  real(kind=4) :: x(*)
  integer :: i,im,in,ipp,j,k,l,ip(MAXsDIM),iu(MAXsDIM,MAXBIT), &
             iv(MAXBIT*MAXsDIM),ix(MAXsDIM),mdeg(MAXsDIM)
  real(kind=8) :: fac
  save :: ip,mdeg,ix,iv,in,fac
  equivalence (iv,iu)
!     DATA ip /0,1,1,2,1,4/, mdeg /1,2,3,3,4,4/, ix /6*0/
!     DATA iv /6*1,3,1,3,3,1,1,5,7,7,3,3,5,15,11,5,15,13,9,156*0/

      ierr = 0

      if (n.lt.0) then
        ip(1) = 0
        ip(2) = 1
        ip(3) = 1
        ip(4) = 2
        ip(5) = 1
        ip(6) = 4
        mdeg(1) = 1
        mdeg(2) = 2
        mdeg(3) = 3
        mdeg(4) = 3
        mdeg(5) = 4
        mdeg(6) = 4
        do i=1,6
           ix(i)=0
           iv(i)=1
        enddo
        iv(7)=3
        iv(8)=1
        iv(9)=3
        iv(10)=3
        iv(11)=1
        iv(12)=1
        iv(13)=5
        iv(14)=7
        iv(15)=7
        iv(16)=3
        iv(17)=3
        iv(18)=5
        iv(19)=15
        iv(20)=11
        iv(21)=5
        iv(22)=15
        iv(23)=13
        iv(24)=9
        do i=25,maxbit*maxsdim
           iv(i)=0
        enddo

        do 14 k=1,MAXsDIM
          do 11 j=1,mdeg(k)
            iu(k,j)=iu(k,j)*2**(MAXBIT-j)
11        continue
          do 13 j=mdeg(k)+1,MAXBIT
            ipp=ip(k)
            i=iu(k,j-mdeg(k))
            i=ieor(i,i/2**mdeg(k))
            do 12 l=mdeg(k)-1,1,-1
              if(iand(ipp,1).ne.0)i=ieor(i,iu(k,j-l))
              ipp=ipp/2
12          continue
            iu(k,j)=i
13        continue
14      continue
        fac=1./2.**MAXBIT
        in=0
      else
        im=in
        do 15 j=1,MAXBIT
          if(iand(im,1).eq.0)goto 1
          im=im/2
15      continue
!       pause 'MAXBIT too small in sobseq'
        ierr = 1
1       im=(j-1)*MAXsDIM
        do 16 k=1,min(n,MAXsDIM)
          ix(k)=ieor(ix(k),iv(im+k))
          x(k)=ix(k)*real(fac)
16      continue
        in=in+1
      endif
      return
 end subroutine sobseq


! ----------------------------------------------------------------------
! This is an old very simple random generator used in some parts of the
! code. It used to be in lib/utilities/func.f
! It is only used now in the genetic algorithm in vchfit.f90
! It could in principle be replaced by randuniform() above
!
! DRANDOM (MB, 9/97)
!
! PURPOSE:
!     Generates random numbers (more or less) evenly distributed between
!     zero and one.
! INPUT PARAMETERS:
!     seed: initializes the random number generator and hence defines
!           the sequence of random numbers (seed=4.0 is a good choice)
! FUNCTION VALUE:
!     a random number between zero and one
! ----------------------------------------------------------------------
!
 function drandom (rseed)

  implicit none

  logical                :: initial
  real(kind=8)           :: drandom,rseed,init,prev
  real(kind=8)           :: pi 

  save    init,prev,initial
  data    initial /.true./

  pi = 4.0d0*atan(1.0d0)

! --- INITIALIZE ON FIRST CALL ---

  if (initial) then
    init = rseed
    prev = rseed
    initial = .false.
  endif

! --- START NEW SEQUENCE IF SEED HAS CHANGED SINCE LAST CALL ---

  if (init .ne. rseed) then
    init = rseed
    prev = rseed
  endif
         
! --- COMPUTE RANDOM NUMBER ---

  prev = exp(5.0d0*log(prev+pi))
  prev = prev-int(prev)
  drandom = prev
      
  return
 end function drandom

end module randommod
