module prob

  use :: iso_c_binding
  use :: fgsl

  implicit none

  private

  integer, parameter :: dp = c_double
  integer, parameter :: i4 = c_int
  integer, parameter :: i8 = c_int64_t

  type, public :: prob_dists
    private
    type(fgsl_rng) :: r
    type(fgsl_rng_type) :: t
    logical :: initialized = .false.
  contains
    private
    procedure, pass(self) :: get_unif_rvs_scalar_double
    procedure, pass(self) :: get_unif_rvs_vector_double
    procedure, pass(self) :: get_unif_rvs_scalar_i4
    procedure, pass(self) :: get_unif_rvs_scalar_i8
    procedure, pass(self) :: get_unif_rvs_vector_i4
    procedure, pass(self) :: get_unif_rvs_vector_i8
    generic, public :: unif_rvs => get_unif_rvs_scalar_double, get_unif_rvs_vector_double, &
                                   get_unif_rvs_scalar_i4, get_unif_rvs_scalar_i8, &
                                   get_unif_rvs_vector_i4, get_unif_rvs_vector_i8
    
    procedure, pass(self) :: get_poisson_rvs_scalar_int4
    procedure, pass(self) :: get_poisson_rvs_vector_int4
    procedure, pass(self) :: get_poisson_rvs_scalar_int8
    procedure, pass(self) :: get_poisson_rvs_vector_int8
    procedure, pass(self) :: get_poisson_rvs_scalar_double
    procedure, pass(self) :: get_poisson_rvs_vector_double
    generic, public :: poisson_rvs => get_poisson_rvs_scalar_int4, get_poisson_rvs_vector_int4, &
                                      get_poisson_rvs_scalar_int8, get_poisson_rvs_vector_int8, &
                                      get_poisson_rvs_scalar_double, get_poisson_rvs_vector_double
    
    procedure, pass(self) :: get_exponential_rvs_scalar
    procedure, pass(self) :: get_exponential_rvs_vector
    generic, public :: exponential_rvs => get_exponential_rvs_scalar, get_exponential_rvs_vector
    
    procedure, pass(self) :: get_bernoulli_rvs_scalar
    procedure, pass(self) :: get_bernoulli_rvs_vector
    generic, public :: bernoulli_int_rvs => get_bernoulli_rvs_scalar, get_bernoulli_rvs_vector
    
    procedure, pass(self) :: get_bernoulli_rvs_scalar_logical
    procedure, pass(self) :: get_bernoulli_rvs_vector_logical
    generic, public :: bernoulli_rvs => get_bernoulli_rvs_scalar_logical, get_bernoulli_rvs_vector_logical

    procedure, public, pass(self) :: sample => ran_sample
    
    procedure, pass(self) :: ran_choose_scalar
    procedure, pass(self) :: ran_choose_i4
    procedure, pass(self) :: ran_choose_i8
    generic, public :: choose => ran_choose_scalar, ran_choose_i4, ran_choose_i8
    
    procedure, public, pass(self) :: get_ran_int_rvs_i4
    procedure, public, pass(self) :: get_ran_int_rvs_i8
    generic, public :: randint => get_ran_int_rvs_i4, get_ran_int_rvs_i8

    procedure, public, pass(self) :: ran_shuffle_i4
    procedure, public, pass(self) :: ran_shuffle_i8
    procedure, public, pass(self) :: ran_shuffle_double
    generic, public :: shuffle => ran_shuffle_i4, ran_shuffle_i8, ran_shuffle_double
    
    procedure, pass(self) :: get_loglogistic_rvs_scalar
    procedure, pass(self) :: get_loglogistic_rvs_vector
    generic, public :: loglogistic_rvs => get_loglogistic_rvs_scalar, get_loglogistic_rvs_vector
    
    procedure, pass(self) :: get_lognormal_rvs_scalar
    procedure, pass(self) :: get_lognormal_rvs_vector
    generic, public :: lognormal_rvs => get_lognormal_rvs_scalar, get_lognormal_rvs_vector

    procedure, pass(self) :: get_normal_rvs_scalar
    procedure, pass(self) :: get_normal_rvs_vector
    generic, public :: normal_rvs => get_normal_rvs_scalar, get_normal_rvs_vector

    procedure, pass(self) :: get_t_rvs_scalar
    procedure, pass(self) :: get_t_rvs_vector
    generic, public :: t_rvs => get_t_rvs_scalar, get_t_rvs_vector

    procedure, pass(self) :: weighted_choose_i4i4, weighted_choose_i4i8, weighted_choose_i8i4, weighted_choose_i8i8
    generic, public :: weighted_choose => weighted_choose_i4i4, weighted_choose_i4i8, weighted_choose_i8i4, weighted_choose_i8i8

    final :: prob_dists_destructor
  end type prob_dists

  interface prob_dists
    module procedure init_prob_dists
  end interface prob_dists

  type(prob_dists), public :: pdists

contains

  function init_prob_dists(seed) result(self)
    implicit none
    integer(kind=c_int64_t), intent(in), optional :: seed
    type(prob_dists) :: self

    integer(kind=c_int64_t) :: seed_

    if (present(seed)) then
      seed_ = seed
    else
      call system_clock(count=seed_)
    end if

    self%t = fgsl_rng_env_setup()
    self%t = fgsl_rng_default
    self%r = fgsl_rng_alloc(self%t)
    call fgsl_rng_set(self%r, seed_)

    self%initialized = .true.
    return
  end function init_prob_dists

  subroutine prob_dists_destructor(self)
    implicit none
    type(prob_dists) :: self

    call fgsl_rng_free(self%r)
    return
  end subroutine prob_dists_destructor

  subroutine get_unif_rvs_scalar_double(self, f, a, b)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(out) :: f

    real(kind=c_double), optional, intent(in) :: a
    real(kind=c_double), optional, intent(in) :: b
    real(kind=c_double) :: a_
    real(kind=c_double) :: b_
    
    a_ = 0.0_dp
    b_ = 1.0_dp
    if (present(a)) a_ = a
    if (present(b)) b_ = b

    f = a_ + (b_-a_)*fgsl_rng_uniform(self%r)
    return
  end subroutine get_unif_rvs_scalar_double

  subroutine get_unif_rvs_vector_double(self, N, rvs, a, b)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int64_t), intent(in) :: N
    real(kind=c_double), dimension(:), allocatable, intent(out) :: rvs

    real(kind=c_double), optional, intent(in) :: a
    real(kind=c_double), optional, intent(in) :: b
    real(kind=c_double) :: a_
    real(kind=c_double) :: b_

    integer(kind=c_int64_t) :: j

    a_ = 0.0_dp
    b_ = 1.0_dp
    if (present(a)) a_ = a
    if (present(b)) b_ = b

    if (allocated(rvs)) deallocate(rvs)
    allocate(rvs(N))

    do j=1,N
      rvs(j) = a_ + (b_-a_)*fgsl_rng_uniform(self%r)
    end do
    return
  end subroutine get_unif_rvs_vector_double

  subroutine get_unif_rvs_scalar_i4(self, f, a, b)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int), intent(out) :: f
    integer(kind=c_int), optional, intent(in) :: a
    integer(kind=c_int), optional, intent(in) :: b
    real(kind=c_double) :: f_r

    integer(kind=c_int) :: a_
    integer(kind=c_int) :: b_

    a_ = 0
    b_ = 1
    if (present(a)) a_ = a
    if (present(b)) b_ = b
    
    f_r = fgsl_rng_uniform(self%r)

    f = a_ + (b_-a_)*nint(f_r, kind=c_int)

    return
  end subroutine get_unif_rvs_scalar_i4

  subroutine get_unif_rvs_vector_i4(self, N, rvs, a, b)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int), intent(in) :: N
    integer(kind=c_int), dimension(:), allocatable, intent(out) :: rvs
    
    integer(kind=c_int), intent(in), optional :: a
    integer(kind=c_int), intent(in), optional :: b
    
    real(kind=c_double), dimension(:) :: f_r(N)

    integer(kind=c_int) :: a_
    integer(kind=c_int) :: b_

    integer :: i

    a_ = 0
    b_ = 1
    if (present(a)) a_ = a
    if (present(b)) b_ = b

    if (allocated(rvs)) deallocate(rvs)
    allocate(rvs(N))

    do i=1,N
      f_r(i) = fgsl_rng_uniform(self%r)
    end do

    rvs(:) = a_ + (b_-a_)*nint(f_r(:), kind=c_int)

    return
  end subroutine get_unif_rvs_vector_i4

  subroutine get_unif_rvs_scalar_i8(self, f, a, b)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int64_t), intent(out) :: f
    integer(kind=c_int64_t), intent(in), optional :: a
    integer(kind=c_int64_t), intent(in), optional :: b
    
    real(kind=c_double) :: f_r

    integer(kind=c_int64_t) :: a_
    integer(kind=c_int64_t) :: b_

    a_ = 0
    b_ = 1
    if (present(a)) a_ = a
    if (present(b)) b_ = b

    f_r = fgsl_rng_uniform(self%r)

    f = a_ + (b_-a_)*nint(f_r, kind=c_int64_t)

    return
  end subroutine get_unif_rvs_scalar_i8

  subroutine get_unif_rvs_vector_i8(self, N, rvs, a, b)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int64_t), intent(in) :: N
    integer(kind=c_int64_t), dimension(:), allocatable, intent(out) :: rvs
    integer(kind=c_int64_t), intent(in), optional :: a
    integer(kind=c_int64_t), intent(in), optional :: b
    
    real(kind=c_double), dimension(:) :: f_r(N)

    integer(kind=c_int64_t) :: a_
    integer(kind=c_int64_t) :: b_

    integer(kind=c_int64_t) :: i

    a_ = 0
    b_ = 1
    if (present(a)) a_ = a
    if (present(b)) b_ = b

    if (allocated(rvs)) deallocate(rvs)
    allocate(rvs(N))

    do i=1,N
      f_r(i) = fgsl_rng_uniform(self%r)
    end do

    rvs(:) = a_ + (b_-a_)*nint(f_r(:), kind=c_int64_t)

    return
  end subroutine get_unif_rvs_vector_i8

  function get_poisson_rvs_scalar_int4(self, mu) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int), intent(in) :: mu
    integer(kind=c_int64_t) :: rvs
    rvs = fgsl_ran_poisson(self%r, real(mu, kind=c_double))
    return
  end function get_poisson_rvs_scalar_int4

  function get_poisson_rvs_vector_int4(self, mu, N) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int), intent(in) :: mu
    integer(kind=c_int64_t), intent(in) :: N
    integer(kind=c_int64_t), dimension(:), allocatable :: rvs
    integer(kind=c_int64_t) :: j

    if (allocated(rvs)) deallocate(rvs)
    allocate(rvs(N))

    do j=1,N
      rvs(j) = fgsl_ran_poisson(self%r, real(mu, kind=c_double))
    end do
    return
  end function get_poisson_rvs_vector_int4

  function get_poisson_rvs_scalar_int8(self, mu) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int64_t), intent(in) :: mu
    integer(kind=c_int64_t) :: rvs
    rvs = fgsl_ran_poisson(self%r, real(mu, kind=c_double))
    return
  end function get_poisson_rvs_scalar_int8

  function get_poisson_rvs_vector_int8(self, mu, N) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int64_t), intent(in) :: mu
    integer(kind=c_int64_t), intent(in) :: N
    integer(kind=c_int64_t), dimension(:), allocatable :: rvs
    integer(kind=c_int64_t) :: j

    if (allocated(rvs)) deallocate(rvs)
    allocate(rvs(N))

    do j=1,N
      rvs(j) = fgsl_ran_poisson(self%r, real(mu, kind=c_double))
    end do
    return
  end function get_poisson_rvs_vector_int8

  function get_poisson_rvs_scalar_double(self, mu) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: mu
    integer(kind=c_int64_t) :: rvs
    rvs = fgsl_ran_poisson(self%r, mu)
    return
  end function get_poisson_rvs_scalar_double

  function get_poisson_rvs_vector_double(self, mu, N) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: mu
    integer(kind=c_int64_t), intent(in) :: N
    integer(kind=c_int64_t), dimension(:), allocatable :: rvs
    integer(kind=c_int64_t) :: j

    if (allocated(rvs)) deallocate(rvs)
    allocate(rvs(N))

    do j=1,N
      rvs(j) = fgsl_ran_poisson(self%r, mu)
    end do
    return
  end function get_poisson_rvs_vector_double

  function get_exponential_rvs_scalar(self, mu) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: mu
    real(kind=c_double) :: rvs
    rvs = fgsl_ran_exponential(self%r, mu)
    return
  end function get_exponential_rvs_scalar

  function get_exponential_rvs_vector(self, mu, N) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: mu
    integer(kind=c_int64_t), intent(in) :: N
    real(kind=c_double), dimension(:), allocatable :: rvs
    integer(kind=c_int64_t) :: j

    if (allocated(rvs)) deallocate(rvs)
    allocate(rvs(N))

    do j=1,N
      rvs(j) = fgsl_ran_exponential(self%r, mu)
    end do
    return
  end function get_exponential_rvs_vector

  function get_lognormal_rvs_scalar(self, mu, sigma) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: mu
    real(kind=c_double), intent(in) :: sigma
    real(kind=c_double) :: rvs
    rvs = fgsl_ran_lognormal(self%r, mu, sigma)
    return
  end function get_lognormal_rvs_scalar

  function get_lognormal_rvs_vector(self, mu, sigma, N) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: mu
    real(kind=c_double), intent(in) :: sigma
    integer(kind=c_int64_t), intent(in) :: N
    real(kind=c_double), dimension(:), allocatable :: rvs
    integer(kind=c_int64_t) :: j

    if (allocated(rvs)) deallocate(rvs)
    allocate(rvs(N))

    do j=1,N
      rvs(j) = fgsl_ran_lognormal(self%r, mu, sigma)
    end do
    return
  end function get_lognormal_rvs_vector

  function get_normal_rvs_scalar(self, mu, sigma) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: mu
    real(kind=c_double), intent(in) :: sigma
    real(kind=c_double) :: rvs
    rvs = mu + fgsl_ran_gaussian(self%r, sigma)
    return
  end function get_normal_rvs_scalar

  function get_normal_rvs_vector(self, mu, sigma, N) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: mu
    real(kind=c_double), intent(in) :: sigma
    integer(kind=c_int64_t), intent(in) :: N
    real(kind=c_double), dimension(:), allocatable :: rvs
    integer(kind=c_int64_t) :: j

    if (allocated(rvs)) deallocate(rvs)
    allocate(rvs(N))

    do j=1,N
      rvs(j) = mu + fgsl_ran_gaussian(self%r, sigma)
    end do
    return
  end function get_normal_rvs_vector

  function get_t_rvs_scalar(self, nu) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: nu
    real(kind=c_double) :: rvs
    rvs = fgsl_ran_tdist(self%r, nu)
    return
  end function get_t_rvs_scalar

  function get_t_rvs_vector(self, nu, N) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: nu
    integer(kind=c_int64_t), intent(in) :: N
    real(kind=c_double), dimension(:), allocatable :: rvs
    integer(kind=c_int64_t) :: j

    if (allocated(rvs)) deallocate(rvs)
    allocate(rvs(N))

    do j=1,N
      rvs(j) = fgsl_ran_tdist(self%r, nu)
    end do
    return
  end function get_t_rvs_vector

  function get_bernoulli_rvs_scalar(self, p) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: p
    integer(kind=c_int64_t) :: rvs
    rvs = fgsl_ran_bernoulli(self%r, p)
    return
  end function get_bernoulli_rvs_scalar

  function get_bernoulli_rvs_vector(self, p, N) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: p
    integer(kind=c_int64_t), intent(in) :: N
    integer(kind=c_int64_t), dimension(:), allocatable :: rvs
    integer(kind=c_int64_t) :: j

    if (allocated(rvs)) deallocate(rvs)
    allocate(rvs(N))

    do j=1,N
      rvs(j) = fgsl_ran_bernoulli(self%r, p)
    end do
    return
  end function get_bernoulli_rvs_vector

  function get_bernoulli_rvs_scalar_logical(self, p) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: p
    integer(kind=c_int64_t) :: draw
    logical :: rvs
    draw = fgsl_ran_bernoulli(self%r, p)
    rvs = merge(.True.,.False., draw==1)
    return
  end function get_bernoulli_rvs_scalar_logical

  function get_bernoulli_rvs_vector_logical(self, p, N) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: p
    integer(kind=c_int64_t), intent(in) :: N
    integer(kind=c_int64_t) :: draw
    logical, dimension(:), allocatable :: rvs
    integer(kind=c_int64_t) :: j

    if (allocated(rvs)) deallocate(rvs)
    allocate(rvs(N))

    do j=1,N
      draw = fgsl_ran_bernoulli(self%r, p)
      rvs(j) = merge(.True.,.False., draw==1)
    end do
    return
  end function get_bernoulli_rvs_vector_logical

  function ran_sample(self,N) result(sample)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int64_t), intent(in) :: N
    integer(kind=c_int64_t), dimension(:) :: sample(N)

    integer(kind=c_int64_t) :: j

    do j=1,N
      sample(j) = fgsl_rng_uniform_int(self%r, int(N, kind=c_long))
    end do
    return
  end function ran_sample

  function ran_choose_scalar(self,array) result(sample)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int64_t), target, dimension(:), intent(in) :: array
    integer(kind=c_int64_t), target, dimension(:), allocatable :: sample

    integer(kind=c_int64_t), target, dimension(:), allocatable :: sample_tmp
    type(c_ptr) :: dest, src
    integer(kind=c_int64_t) :: f
    integer(kind=c_int64_t) :: k

    allocate(sample_tmp(1))

    src = c_loc(array(1))
    dest = c_loc(sample_tmp(1))

    k = int(size(array), kind=c_int64_t)

    f = fgsl_ran_choose(self%r, dest, int(1, kind=c_int64_t), src, k, sizeof(array(1)))

    sample = sample_tmp(1)

    return

  end function ran_choose_scalar

  function ran_choose_i4(self,array,N) result(sample)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int64_t), target, dimension(:), intent(in) :: array
    integer(kind=c_int), intent(in) :: N
    integer(kind=c_int64_t), target, dimension(:), allocatable :: sample

    type(c_ptr) :: dest, src
    integer(kind=c_int64_t) :: f
    integer(kind=c_int64_t) :: k

    if (allocated(sample)) deallocate(sample)
    allocate(sample(N))

    src = c_loc(array(1))
    dest = c_loc(sample(1))

    k = int(size(array), kind=c_int64_t)

    f = fgsl_ran_choose(self%r, dest, int(N, kind=c_int64_t), src, k, sizeof(array(1)))

    return

  end function ran_choose_i4

  function ran_choose_i8(self,array,N) result(sample)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int64_t), target, dimension(:), intent(in) :: array
    integer(kind=c_int64_t), intent(in) :: N
    integer(kind=c_int64_t), target, dimension(:), allocatable :: sample

    type(c_ptr) :: dest, src
    integer(kind=c_int64_t) :: f
    integer(kind=c_int64_t) :: k

    if (allocated(sample)) deallocate(sample)
    allocate(sample(N))

    k = int(size(array), kind=c_int64_t)

    if (N==k) then
        sample = array
    else
      src = c_loc(array(1))
      dest = c_loc(sample(1))
      f = fgsl_ran_choose(self%r, dest, N, src, k, sizeof(array(1)))
    end if

    return

  end function ran_choose_i8

  function ran_shuffle_i4(self,array) result(array_out)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int), dimension(:), intent(in) :: array
    integer(kind=c_int), dimension(:), allocatable :: array_out

    integer(kind=c_int64_t), target, dimension(:), allocatable :: a

    integer(kind=c_int) :: N

    type(c_ptr) :: base

    N = size(array)

    if (allocated(array_out)) deallocate(array_out)
    allocate(array_out(N))

    allocate(a(N))
    a = int(array, kind=c_int64_t)

    base = c_loc(a(1))

    call fgsl_ran_shuffle(self%r, base, int(N, kind=c_int64_t), sizeof(a(1)))

    array_out = int(a, kind=c_int)

    return
  end function ran_shuffle_i4

  function ran_shuffle_i8(self,array) result(array_out)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int64_t), dimension(:), intent(in) :: array
    integer(kind=c_int64_t), dimension(:), allocatable :: array_out

    integer(kind=c_int64_t) :: N

    integer(kind=c_int64_t), target, dimension(:), allocatable :: a

    type(c_ptr) :: base

    N = size(array)

    if (allocated(array_out)) deallocate(array_out)
    allocate(array_out(N))

    allocate(a(N))
    a = array

    base = c_loc(a(1))

    call fgsl_ran_shuffle(self%r, base, N, sizeof(a(1)))

    array_out = a

    return
  end function ran_shuffle_i8

  function ran_shuffle_double(self,array) result(array_out)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), dimension(:), intent(in) :: array
    real(kind=c_double), dimension(:), allocatable :: array_out

    integer(kind=c_int64_t) :: N

    type(c_ptr) :: base

    real(kind=c_double), target, dimension(:), allocatable :: a

    N = size(array)

    if (allocated(array_out)) deallocate(array_out)
    allocate(array_out(N))

    allocate(a(N))
    a = array

    base = c_loc(a(1))

    call fgsl_ran_shuffle(self%r, base, N, sizeof(a(1)))

    array_out = a

    return
  end function ran_shuffle_double

  function get_ran_int_rvs_i4(self,a,b) result(val)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int), intent(in) :: a
    integer(kind=c_int), intent(in) :: b
    integer(kind=c_int64_t) :: val

    val = floor(fgsl_ran_flat(self%r,real(a, kind=c_double),real(b, kind=c_double)))

    return

  end function get_ran_int_rvs_i4

  function get_ran_int_rvs_i8(self,a,b) result(val)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int64_t), intent(in) :: a
    integer(kind=c_int64_t), intent(in) :: b
    integer(kind=c_int64_t) :: val

    val = floor(fgsl_ran_flat(self%r,real(a, kind=c_double),real(b, kind=c_double)))

    return

  end function get_ran_int_rvs_i8

  function get_loglogistic_rvs_scalar(Self,scale,shape) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: scale
    real(kind=c_double), intent(in) :: shape
    real(kind=c_double) :: rvs

    real(kind=c_double) :: logistic_rvs

    logistic_rvs = fgsl_ran_logistic(self%r, 1.0_dp/shape)

    rvs = scale*exp(logistic_rvs)

    return

  end function get_loglogistic_rvs_scalar

  function get_loglogistic_rvs_vector(Self,scale,shape,N) result(rvs)
    implicit none
    class(prob_dists), intent(in) :: self
    real(kind=c_double), intent(in) :: scale
    real(kind=c_double), intent(in) :: shape
    integer(kind=c_int64_t), intent(in) :: N
    real(kind=c_double), dimension(:), allocatable :: rvs

    integer(kind=c_int64_t) :: j

    real(kind=c_double) :: logistic_rvs

    allocate(rvs(N))

    do j=1,N
      logistic_rvs = fgsl_ran_logistic(self%r, 1.0_dp/shape)
      rvs(j) = scale*exp(logistic_rvs)
    end do

    return

  end function get_loglogistic_rvs_vector

  function weighted_choose_i4i4(self, n, items, weights) result(choice)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int), intent(in) :: n
    integer(kind=c_int), dimension(:), intent(in) :: items
    real(kind=c_double), dimension(:), intent(in) :: weights

    integer(kind=c_int), dimension(:) :: choice(n)
    real(kind=c_double), dimension(:) :: cum(size(weights))
    real(kind=c_double) :: rnum

    integer :: i, j
    integer :: M
    M = size(weights)
    
    cum(1) = weights(1)
    do j=2, M
      cum(j) = cum(j-1) + weights(j)
    end do

    do i=1, n
      call self%unif_rvs(rnum)
      do j=1, M
        if (rnum < cum(j)) then
          choice(i) = items(j)
          exit
        end if
      end do
    end do
    return
  end function weighted_choose_i4i4

  function weighted_choose_i4i8(self, n, items, weights) result(choice)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int), intent(in) :: n
    integer(kind=c_int64_t), dimension(:), intent(in) :: items
    real(kind=c_double), dimension(:), intent(in) :: weights
    integer(kind=c_int64_t), dimension(:) :: choice(n)

    real(kind=c_double), dimension(:) :: cum(size(weights))
    real(kind=c_double) :: rnum

    integer :: i, j
    integer :: M
    M = size(weights)
    
    cum(1) = weights(1)
    do j=2, M
      cum(j) = cum(j-1) + weights(j)
    end do

    do i=1, n
      call self%unif_rvs(rnum)
      do j=1, M
        if (rnum < cum(j)) then
          choice(i) = items(j)
          exit
        end if
      end do
    end do
    return
  end function weighted_choose_i4i8

  function weighted_choose_i8i4(self, n, items, weights) result(choice)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int64_t), intent(in) :: n
    integer(kind=c_int), dimension(:), intent(in) :: items
    real(kind=c_double), dimension(:), intent(in) :: weights
    integer(kind=c_int64_t), dimension(:) :: choice(n)

    real(kind=c_double), dimension(:) :: cum(size(weights))
    real(kind=c_double) :: rnum

    integer :: i, j
    integer :: M
    M = size(weights)
    
    cum(1) = weights(1)
    do j=2, M
      cum(j) = cum(j-1) + weights(j)
    end do

    do i=1, int(n, kind=c_int)
      call self%unif_rvs(rnum)
      do j=1, M
        if (rnum < cum(j)) then
          choice(i) = items(j)
          exit
        end if
      end do
    end do
    return
  end function weighted_choose_i8i4

  function weighted_choose_i8i8(self, n, items, weights) result(choice)
    implicit none
    class(prob_dists), intent(in) :: self
    integer(kind=c_int64_t), intent(in) :: n
    integer(kind=c_int64_t), dimension(:), intent(in) :: items
    real(kind=c_double), dimension(:), intent(in) :: weights
    integer(kind=c_int64_t), dimension(:) :: choice(n)

    real(kind=c_double), dimension(:) :: cum(size(weights))
    real(kind=c_double) :: rnum

    integer :: i, j
    integer :: M
    M = size(weights)
    
    cum(1) = weights(1)
    do j=2, M
      cum(j) = cum(j-1) + weights(j)
    end do

    do i=1, int(n, kind=c_int)
      call self%unif_rvs(rnum)
      do j=1, M
        if (rnum < cum(j)) then
          choice(i) = items(j)
          exit
        end if
      end do
    end do
    return
  end function weighted_choose_i8i8

end module prob