!=============================================================================!
!                          A L G O R I T H M S                                !
!=============================================================================!
!-----------------------------------------------------------------------------!
! Simple, general purpose algorithms that may be required at various places   !
! within the code, e.g. portable random number generators, etc.               !
!-----------------------------------------------------------------------------!
! This cut-down version MIJP 05/10/2013                                       !
!=============================================================================!

module algor
  use constants, only : dp                      !Minimal useage where practical

  implicit none                                 !Impose strong typing

  private                                       !Everything is private ...

  !---------------------------------------------------------------------------!
  !                       P u b l i c   R o u t i n e s                       !
  !---------------------------------------------------------------------------!
  public :: algor_set_random_seed               !... unless exposed here.
  public :: algor_uniform_random
  public :: algor_erf
  public :: algor_erfc
  public :: algor_sizeof

  interface algor_sizeof
     module procedure algor_sizeof_logical
     module procedure algor_sizeof_int
     module procedure algor_sizeof_real
     module procedure algor_sizeof_cmplx    
  end interface algor_sizeof

  !---------------------------------------------------------------------------!
  !                        P u b l i c   V a r i a b l e s                    !
  !---------------------------------------------------------------------------!

  !---------------------------------------------------------------------------!
  !                      P r i v a t e   R o u t i n e s                      !
  !---------------------------------------------------------------------------!

  !---------------------------------------------------------------------------!
  !                      P r i v a t e   V a r i a b l e s                    !
  !---------------------------------------------------------------------------!
  !For the uniform random number generator, we presume a 32-bit integer in order
  !to do bit-shuffling, so define and store here.
  integer, parameter         :: bit_32=selected_int_kind(9)
  integer(kind=bit_32), save :: ix,iy                     !intermediate randoms
  integer, save              :: algor_seed=223257167      !start random sequence
  logical, save              :: random_initialised=.false.

contains

  subroutine algor_set_random_seed(seed, over_ride)
    !=========================================================================!
    ! Set module internal random number generator seed.                       !
    !-------------------------------------------------------------------------!
    ! Arguments:                                                              !
    !   seed intent(in)=> deterministic seed.                                 !
    !   over_ride intent(in), optional: allows resetting of random_seed       !
    !-------------------------------------------------------------------------!
    ! Parent module variables used:                                           !
    !   iseed, ix, iy and random_initialised are set here.                    !
    !-------------------------------------------------------------------------!
    ! Key Internal Variables:                                                 !
    !   iseed as starting 32-bit integer seed for all generators              !
    !   ix as seed 32-bit integer in Marsaglia generator                      !
    !   iy as seed 32-bit integer in Park-Miller generator                    !
    !-------------------------------------------------------------------------!
    ! Written by Matt Probert, v1.1, 18/04/2007                               !
    !=========================================================================!
    implicit none
    integer, intent(inout) :: seed             !may or may not be 32 bit...
    logical, intent(in), optional :: over_ride

    !32 bit integer version of seed
    integer(kind=bit_32) :: iseed              !random number seed

    !F90 random number stuff
    integer :: rand_size
    integer, allocatable, dimension(:) :: rand_seed
    real (kind=dp) :: rtime

    integer        :: ierr

    if (bit_size(iseed)/=32) stop'Error in algor_set_random_seed - 32-bit integer kind type problem'

    !pick up optional over_ride
    if (present(over_ride)) then
     !Over-ride an already set seed?
      if (over_ride) random_initialised=.false.
    end if

    !We only want to set the seed if it has not already been set
    !and the flag is set to .false. by the compiler
    if (.not.random_initialised) then

       rtime = real(abs(seed),kind=dp)                   !convert seed to real

       !and then convert to bit_32 size integer
       iseed = int(rtime,kind=bit_32)                    !must fit within huge(1)

       ix=ieor(777755555_bit_32,iseed)                   !Marsaglia generator
       iy=ior(ieor(888889999_bit_32,iseed),1_bit_32)     !Parks-Miller generator

       !Nah - lets not give away too many secrets - use F90 random number generator instead
       call random_seed(size=rand_size)
       allocate(rand_seed(1:rand_size),stat=ierr)
       if (ierr/=0) stop 'Error in allocating rand_seed'
       rand_seed=algor_seed
       call random_seed(put=rand_seed)
       deallocate(rand_seed,stat=ierr)
       if (ierr/=0) stop 'Error in deallocating rand_seed'
       !now random number generator is ready to use ...

       random_initialised=.true.                         !set flag
       seed = int(rtime)                                 !return the seed that was used in default integer
    end if

    return
  end subroutine algor_set_random_seed

  function algor_uniform_random()
    !=========================================================================!
    ! Return a single random deviate ~ uniform [0,1].                         !
    ! Based on Park-Miller "minimal standard" generator with Schrage's method !
    !  to do 32-bit multiplication without requiring higher precision, plus   !
    !  additional Marsaglia shift to suppress any weaknesses & correlations.  !
    ! Using two independent methods greatly increases the period of the       !
    !  generator, s.t. resulting period ~2*10^18                              !
    ! NB Routine is only set to work with 32 bit integers!                    !
    !-------------------------------------------------------------------------!
    ! References:                                                             !
    !   S.K. Park and K.W. Miller, Commun. ACM, 31, p1192-1201 (1988)         !
    !   L. Schrage, ACM Trans. Math. Soft., 5, p132-138 (1979)                !
    !   G. Marsaglia, Linear Algebra and its Applications, 67, p147-156 (1985)!
    !-------------------------------------------------------------------------!
    ! Arguments:                                                              !
    !   none                                                                  !
    !-------------------------------------------------------------------------!
    ! Return value:                                                           !
    !   algor_uniform_random => required random deviate                       !
    !-------------------------------------------------------------------------!
    ! Parent module variables used:                                           !
    !   ix as next 32-bit integer in Marsaglia generator (updated)            !
    !   iy as next 32-bit integer in Park-Miller generator (updated)          !
    !-------------------------------------------------------------------------!
    ! Modules used:                                                           !
    !   none                                                                  !
    !-------------------------------------------------------------------------!
    ! Key Internal Variables:                                                 !
    !   none                                                                  !
    !-------------------------------------------------------------------------!
    ! Necessary conditions:                                                   !
    !   none                                                                  !
    !-------------------------------------------------------------------------!
    ! Written by Matt Probert, v0.1, 01/07/2000                               !
    !=========================================================================!
    implicit none
    real(kind=dp) :: algor_uniform_random

    if (.not.random_initialised) call algor_set_random_seed(algor_seed)

    !do Marsaglia shift sequence, period 2^32-1, to get new ix
    !Nah - lets not give away too many secrets - use standard fortran generator instead

    !Park-Miller sequence, period iy_max-1, to get new iy
    !Nah - lets not give away too many secrets - use standard fortran generator instead

    !Combine ix and iy to get new random number, rescale onto range [0,1]
    !Nah - lets not give away too many secrets - use standard fortran generator instead

    call random_number(algor_uniform_random)           !0<random_number<1

    return
  end function algor_uniform_random

  function algor_erf(x)
    !=========================================================================!
    ! Calculate an accurage approximation to the error function, erf(x),      !
    !    erf(x)=(2/sqrt(pi))*integral(0->x)[exp(-t^2)]dt.                     !
    !    Based upon parameterization given in NSWC Mathematics Library.       !
    !    NB This is machine portable and not dependent on a system function.  !
    !-------------------------------------------------------------------------!
    ! Arguments:                                                              !
    !   x, intent=in, argument of erf to evaluate.                            !
    !-------------------------------------------------------------------------!
    ! Parent module variables used:                                           !
    !   none.                                                                 !
    !-------------------------------------------------------------------------!
    ! Modules used:                                                           !
    !   none.                                                                 !
    !-------------------------------------------------------------------------!
    ! Key Internal Variables:                                                 !
    !   none.                                                                 !
    !-------------------------------------------------------------------------!
    ! Necessary conditions:                                                   !
    !   none.                                                                 !
    !-------------------------------------------------------------------------!
    ! Written by Matt Probert, v0.1, 25/09/2001                               !
    !=========================================================================!

    implicit none

    real(kind=dp), intent(in) :: x
    real(kind=dp)             :: algor_erf

    !expansion parameters
    real(kind=dp), dimension(1:5) :: a=(/ 0.771058495001320E-04_dp, -0.133733772997339E-02_dp,   &
         & 0.323076579225834E-01_dp, 0.479137145607681E-01_dp, 0.128379167095513E+00_dp /)
    real(kind=dp), dimension(1:3) :: b=(/ 0.301048631703895E-02_dp,  0.538971687740286E-01_dp,   &
         & 0.375795757275549E+00_dp /)
    real(kind=dp), dimension(1:8) :: p=(/-1.36864857382717E-07_dp,  5.64195517478974E-01_dp, &
         & 7.21175825088309E+00_dp,  4.31622272220567E+01_dp,  1.52989285046940E+02_dp, &
         & 3.39320816734344E+02_dp,  4.51918953711873E+02_dp,  3.00459261020162E+02_dp /)
    real(kind=dp), dimension(1:8) :: q=(/ 1.00000000000000E+00_dp,  1.27827273196294E+01_dp, &
         & 7.70001529352295E+01_dp,  2.77585444743988E+02_dp,  6.38980264465631E+02_dp,  &
         & 9.31354094850610E+02_dp,  7.90950925327898E+02_dp,  3.00459260956983E+02_dp /)
    real(kind=dp), dimension(1:5) :: r=(/ 2.10144126479064E+00_dp,  2.62370141675169E+01_dp, &
         & 2.13688200555087E+01_dp,  4.65807828718470E+00_dp,  2.82094791773523E-01_dp /)
    real(kind=dp), dimension(1:4) :: s=(/ 9.41537750555460E+01_dp,  1.87114811799590E+02_dp, &
         & 9.90191814623914E+01_dp,  1.80124575948747E+01_dp /)
    real(kind=dp)                 :: c=0.564189583547756_dp

    !local vars
    real(kind=dp) :: ax, bot, t, top, x2

    ax=abs(x)
    x2=x*x
    
    if (ax<0.5_dp) then
       t=x2
       top=((((a(1)*t + a(2))*t + a(3))*t + a(4))*t + a(5)) + 1.0_dp
       bot=((  b(1)*t + b(2))*t + b(3))*t + 1.0_dp
       algor_erf=ax*(top/bot)

    else if (ax<4.0_dp) then
       top=((((((p(1) *ax + p(2))*ax + p(3))*ax + p(4))*ax + p(5))*ax  &
            &  + p(6))*ax + p(7))*ax + p(8)
       bot=((((((q(1) *ax + q(2))*ax + q(3))*ax + q(4))*ax + q(5))*ax  &
            &  + q(6))*ax + q(7))*ax + q(8)
       algor_erf=0.5_dp + (0.5_dp - exp(-x2)*top/bot)

    else if (ax<5.8_dp) then
       t=1.0_dp/x2
       top=(((r(1)*t + r(2))*t + r(3))*t + r(4))*t + r(5)
       bot=(((s(1)*t + s(2))*t + s(3))*t + s(4))*t + 1.0_dp
       algor_erf=0.5_dp + (0.5_dp - exp(-x2)*((c - top/(x2*bot))/ax))

    else
       algor_erf=1.0_dp            !large |x| limit

    end if

    !now put the sign back in
    algor_erf=sign(algor_erf,x)

    return
  end function algor_erf

  function algor_erfc(x)
    !=========================================================================!
    ! Calculate an accurage approximation to the complemenary error function, !
    !    erfc(x)=(2/sqrt(pi))*integral(x->infinity)[exp(-t^2)]dt.             !
    !    Based upon parameterization given in NSWC Mathematics Library.       !
    !    NB This is machine portable and not dependent on a system function.  !
    !-------------------------------------------------------------------------!
    ! Arguments:                                                              !
    !   x, intent=in, argument of erfc to evaluate.                           !
    !-------------------------------------------------------------------------!
    ! Parent module variables used:                                           !
    !   none.                                                                 !
    !-------------------------------------------------------------------------!
    ! Modules used:                                                           !
    !   none.                                                                 !
    !-------------------------------------------------------------------------!
    ! Key Internal Variables:                                                 !
    !   none.                                                                 !
    !-------------------------------------------------------------------------!
    ! Necessary conditions:                                                   !
    !   none.                                                                 !
    !-------------------------------------------------------------------------!
    ! Written by Matt Probert, v0.1, 25/09/2001                               !
    !=========================================================================!
    implicit none

    real(kind=dp), intent(in) :: x
    real(kind=dp)             :: algor_erfc

    !expansion parameters
    real(kind=dp), dimension(1:5) :: a=(/  0.771058495001320E-04_dp, -0.133733772997339E-02_dp,  &
         &  0.323076579225834E-01_dp,  0.479137145607681E-01_dp,  0.128379167095513E+00_dp /)
    real(kind=dp), dimension(1:3) :: b=(/  0.301048631703895E-02_dp,  0.538971687740286E-01_dp,  &
         &  0.375795757275549E+00_dp /)
    real(kind=dp), dimension(1:8) :: p=(/ -1.36864857382717E-07_dp,  5.64195517478974E-01_dp, &
         &  7.21175825088309E+00_dp,  4.31622272220567E+01_dp,  1.52989285046940E+02_dp,  &
         &  3.39320816734344E+02_dp,  4.51918953711873E+02_dp,  3.00459261020162E+02_dp /)
    real(kind=dp), dimension(1:8) :: q=(/  1.00000000000000E+00_dp,  1.27827273196294E+01_dp, &
         &  7.70001529352295E+01_dp, 2.77585444743988E+02_dp,  6.38980264465631E+02_dp, &
         &  9.31354094850610E+02_dp, 7.90950925327898E+02_dp,  3.00459260956983E+02_dp /)
    real(kind=dp), dimension(1:5) :: r=(/  2.10144126479064E+00_dp,  2.62370141675169E+01_dp, &
         &  2.13688200555087E+01_dp, 4.65807828718470E+00_dp,  2.82094791773523E-01_dp /)
    real(kind=dp), dimension(1:4) :: s=(/ 9.41537750555460E+01_dp, 1.87114811799590E+02_dp, &
         &  9.90191814623914E+01_dp, 1.80124575948747E+01_dp /)
    real(kind=dp)                 :: c = .564189583547756_dp

    ! Local variables
    real(kind=dp) :: ax, x2, t, top, bot

    ax=abs(x)
    x2=x*x
    if (ax<0.5_dp) then
       t=x2
       top=((((a(1)*t + a(2))*t + a(3))*t + a(4))*t + a(5)) + 1.0_dp
       bot=((  b(1)*t + b(2))*t + b(3)) * t + 1.0_dp
       algor_erfc=0.5_dp + (0.5_dp-x*(top/bot))

    else if (ax<4.0_dp) then
       top=((((((p(1) *ax + p(2))*ax + p(3))*ax + p(4))*ax + p(5))*ax + p(6))*ax  &
            &  + p(7))*ax + p(8)
       bot=((((((q(1) *ax + q(2))*ax + q(3))*ax + q(4))*ax + q(5))*ax + q(6))*ax  &
            &  + q(7))*ax + q(8)
       algor_erfc=exp(-x2)*top/bot
       if (x<0.0_dp) algor_erfc=2.0_dp - algor_erfc

    else
       if (x<=-5.6_dp) then
          algor_erfc=2.0_dp    !large negative x limit

       else if (x>100.0_dp) then
          algor_erfc=0.0_dp    !large positive x limit

       else
          t=1.0_dp/x2
          top=(((r(1)*t + r(2))*t + r(3))*t + r(4)) * t + r(5)
          bot=(((s(1)*t + s(2))*t + s(3))*t + s(4)) * t + 1.0_dp
          algor_erfc=exp(-x2)*(c-t*top/bot)/ax
          if (x<0.0_dp) algor_erfc=2.0_dp - algor_erfc

       end if
    end if

    return
  end function algor_erfc

  integer function algor_sizeof_logical(bool)
    !===========================================================================!
    ! This function take an logical bool, and returns the number of bytes       !
    ! required to store it in memory.                                           !
    !---------------------------------------------------------------------------!
    ! Written by Phil Hasnip, version 0.2, 04/11/05                             !
    !===========================================================================!
    implicit none

    logical, intent(in) :: bool

    !we presume that a FORTRAN logical is always the same size as a default integer

    algor_sizeof_logical = algor_sizeof(1)

    return

  end function algor_sizeof_logical

  integer function algor_sizeof_int(scalar)
    !===========================================================================!
    ! This function take an integer scalar, and returns the number of bytes     !
    ! required to store it in memory.                                           !
    !---------------------------------------------------------------------------!
    ! Written by Phil Hasnip, version 0.2, 27/02/01                             !
    !===========================================================================!
    implicit none

    integer, intent(in) :: scalar

    algor_sizeof_int = (7+bit_size(scalar))/8

    return

  end function algor_sizeof_int

  integer function algor_sizeof_real(scalar)
    !===========================================================================!
    ! This function take a real scalar, and returns the number of bytes         !
    ! required to store it in memory.                                           !
    !---------------------------------------------------------------------------!
    ! Written by Phil Hasnip, version 0.2, 27/02/01                             !
    !===========================================================================!
    implicit none

    real(kind=dp), intent(in) :: scalar

    real(kind=dp) :: bits,b,p,emax,emin
    integer       :: bytes

    ! assume a real no. has a sign bit, p digits in base b in the mantissa, and log_2(e_max-e_min) 
    ! digits in the exponent. In binary there may be a hidden bit.

    bits=0.0_dp

    b=real(radix(scalar),kind=dp)
    p=real(digits(scalar),kind=dp)

    bits=bits+p*log(b)/log(2.0_dp)

    emax=real(maxexponent(scalar),kind=dp)
    emin=real(minexponent(scalar),kind=dp)

    bits=bits+log(emax-emin)/log(2.0_dp)

    if(b.eq.2) then
       bits=bits-1.0_dp           ! assume hidden bit if binary
    end if

    bytes=ceiling(bits/8.0_dp)   ! we round up to whole bytes

    algor_sizeof_real = bytes

    return

  end function algor_sizeof_real

  integer function algor_sizeof_cmplx(scalar)
    !===========================================================================!
    ! This function take a complex scalar, and returns the number of bytes      !
    ! required to store it in memory.                                           !
    !---------------------------------------------------------------------------!
    ! Written by Phil Hasnip, version 0.2, 27/02/01                             !
    !===========================================================================!
    implicit none

    complex(kind=dp), intent(in) :: scalar

    ! Assumes complex numbers use the same storage as 2 real numbers

    algor_sizeof_cmplx = 2*algor_sizeof(1.0_dp)

    return 

  end function algor_sizeof_cmplx

end module algor
