#:include "common.fypp"
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
#:set ALL_KINDS_TYPES = INT_KINDS_TYPES + RC_KINDS_TYPES
module stdlib_stats_distribution_uniform
use stdlib_kinds, only : sp, dp, xdp, qp, int8, int16, int32, int64
use stdlib_error, only : error_stop
use stdlib_random, only : dist_rand
implicit none
private
real(dp), parameter :: MESENNE_NUMBER = 1.0_dp / (2.0_dp ** 53 - 1.0_dp)
integer(int64), parameter :: INT_ONE = 1_int64
public :: rvs_uniform
public :: pdf_uniform
public :: cdf_uniform
public :: shuffle
interface rvs_uniform
!! version: experimental
!!
!! Get uniformly distributed random variate for integer, real and complex
!! variables.
!! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html#
!! rvs_uniform-uniform-distribution-random-variates))
module procedure rvs_unif_0_rsp ! 0 dummy variable
#:for k1, t1 in ALL_KINDS_TYPES
module procedure rvs_unif_1_${t1[0]}$${k1}$ ! 1 dummy variable
#:endfor
#:for k1, t1 in ALL_KINDS_TYPES
module procedure rvs_unif_${t1[0]}$${k1}$ ! 2 dummy variables
#:endfor
#:for k1, t1 in ALL_KINDS_TYPES
module procedure rvs_unif_array_${t1[0]}$${k1}$ ! 3 dummy variables
#:endfor
end interface rvs_uniform
interface pdf_uniform
!! version: experimental
!!
!! Get uniform distribution probability density (pdf) for integer, real and
!! complex variables.
!! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html#
!! pdf_uniform-uniform-probability-density-function))
#:for k1, t1 in ALL_KINDS_TYPES
module procedure pdf_unif_${t1[0]}$${k1}$
#:endfor
end interface pdf_uniform
interface cdf_uniform
!! version: experimental
!!
!! Get uniform distribution cumulative distribution function (cdf) for integer,
!! real and complex variables.
!! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html#
!! cdf_uniform-uniform-cumulative-distribution-function))
!!
#:for k1, t1 in ALL_KINDS_TYPES
module procedure cdf_unif_${t1[0]}$${k1}$
#:endfor
end interface cdf_uniform
interface shuffle
!! version: experimental
!!
!! Fisher-Yates shuffle algorithm for a rank one array of integer, real and
!! complex variables.
!! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html#
!! shuffle-using-fisher-yates-algorithm-to-generate-a-random-permutation-of-a-list))
!!
#:for k1, t1 in ALL_KINDS_TYPES
module procedure shuffle_${t1[0]}$${k1}$
#:endfor
end interface shuffle
contains
#:for k1, t1 in INT_KINDS_TYPES
impure elemental function rvs_unif_1_${t1[0]}$${k1}$(scale) result(res)
!
! Uniformly distributed integer in [0, scale]
! Bitmask with rejection
! https://www.pcg-random.org/posts/bounded-rands.html
!
! Fortran 90 translated from c by Jim-215-fisher
!
${t1}$, intent(in) :: scale
${t1}$ :: res, u, mask
integer :: zeros, bits_left, bits
if(scale <= 0_${k1}$) call error_stop("Error(rvs_unif_1): Uniform" &
//" distribution scale parameter must be positive")
zeros = leadz(scale)
bits = bit_size(scale) - zeros
mask = shiftr(not(0_${k1}$), zeros)
L1 : do
u = dist_rand(scale)
res = iand(u, mask)
if(res <= scale) exit L1
bits_left = zeros
L2 : do
if(bits_left < bits) exit L2
u = shiftr(u, bits)
res = iand(u, mask)
if(res <= scale) exit L1
bits_left = bits_left - bits
end do L2
end do L1
end function rvs_unif_1_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
impure elemental function rvs_unif_${t1[0]}$${k1}$(loc, scale) result(res)
!
! Uniformly distributed integer in [loc, loc + scale]
!
${t1}$, intent(in) :: loc, scale
${t1}$ :: res
if(scale <= 0_${k1}$) call error_stop("Error(rvs_unif): Uniform" &
//" distribution scale parameter must be positive")
res = loc + rvs_unif_1_${t1[0]}$${k1}$(scale)
end function rvs_unif_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in REAL_KINDS_TYPES
impure elemental function rvs_unif_0_${t1[0]}$${k1}$( ) result(res)
!
! Uniformly distributed float in [0,1]
! Based on the paper by Frederic Goualard, "Generating Random Floating-
! Point Numbers By Dividing Integers: a Case Study", Proceedings of
! ICCS 2020, June 2020, Amsterdam, Netherlands
!
${t1}$ :: res
integer(int64) :: tmp
tmp = shiftr(dist_rand(INT_ONE), 11) ! Get random from [0,2^53-1]
res = real(tmp * MESENNE_NUMBER, kind = ${k1}$) ! convert to [0,1]
end function rvs_unif_0_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in REAL_KINDS_TYPES
impure elemental function rvs_unif_1_${t1[0]}$${k1}$(scale) result(res)
!
! Uniformly distributed float in [0, scale]
!
${t1}$, intent(in) :: scale
${t1}$ :: res
if(scale == 0._${k1}$) call error_stop("Error(rvs_unif_1): " &
//"Uniform distribution scale parameter must be non-zero")
res = scale * rvs_unif_0_${t1[0]}$${k1}$( )
end function rvs_unif_1_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in REAL_KINDS_TYPES
impure elemental function rvs_unif_${t1[0]}$${k1}$(loc, scale) result(res)
!
! Uniformly distributed float in [loc, loc + scale]
!
${t1}$, intent(in) :: loc, scale
${t1}$ :: res
if(scale == 0._${k1}$) call error_stop("Error(rvs_unif): " &
//"Uniform distribution scale parameter must be non-zero")
res = loc + scale * rvs_unif_0_${t1[0]}$${k1}$( )
end function rvs_unif_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in CMPLX_KINDS_TYPES
impure elemental function rvs_unif_1_${t1[0]}$${k1}$(scale) result(res)
!
! Uniformly distributed complex in [(0,0i), (scale, i(scale))]
! The real part and imaginary part are independent of each other, so that
! the joint distribution is on an unit square [(0,0i), (scale,i(scale))]
!
${t1}$, intent(in) :: scale
${t1}$ :: res
real(${k1}$) :: r1, tr, ti
if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(rvs_uni_" &
//"1): Uniform distribution scale parameter must be non-zero")
r1 = rvs_unif_0_r${k1}$( )
if(scale % re == 0.0_${k1}$) then
ti = scale % im * r1
tr = 0.0_${k1}$
else if(scale % im == 0.0_${k1}$) then
tr = scale % re * r1
ti = 0.0_${k1}$
else
tr = scale % re * r1
r1 = rvs_unif_0_r${k1}$( )
ti = scale % im * r1
end if
res = cmplx(tr, ti, kind=${k1}$)
end function rvs_unif_1_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in CMPLX_KINDS_TYPES
impure elemental function rvs_unif_${t1[0]}$${k1}$(loc, scale) result(res)
!
! Uniformly distributed complex in [(loc,iloc), (loc + scale, i(loc +
! scale))].
! The real part and imaginary part are independent of each other, so that
! the joint distribution is on an unit square [(loc,iloc), (loc + scale,
! i(loc + scale))]
!
${t1}$, intent(in) :: loc, scale
${t1}$ :: res
real(${k1}$) :: r1, tr, ti
if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(rvs_uni_" &
//"): Uniform distribution scale parameter must be non-zero")
r1 = rvs_unif_0_r${k1}$( )
if(scale % re == 0.0_${k1}$) then
tr = loc % re
ti = loc % im + scale % im * r1
else if(scale % im == 0.0_${k1}$) then
tr = loc % re + scale % re * r1
ti = loc % im
else
tr = loc % re + scale % re * r1
r1 = rvs_unif_0_r${k1}$( )
ti = loc % im + scale % im * r1
end if
res = cmplx(tr, ti, kind=${k1}$)
end function rvs_unif_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
function rvs_unif_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res)
integer, intent(in) :: array_size
${t1}$, intent(in) :: loc, scale
${t1}$ :: res(array_size)
${t1}$ :: u, mask, nn
integer :: i, zeros, bits_left, bits
if(scale == 0_${k1}$) call error_stop("Error(rvs_unif_array): " &
//"Uniform distribution scale parameter must be non-zero")
zeros = leadz(scale)
bits = bit_size(scale) - zeros
mask = shiftr(not(0_${k1}$), zeros)
do i = 1, array_size
L1 : do
u = dist_rand(scale)
nn = iand(u, mask)
if(nn <= scale) exit L1
bits_left = zeros
L2 : do
if(bits_left < bits) exit L2
u = shiftr(u, bits)
nn = iand(u, mask)
if(nn <= scale) exit L1
bits_left = bits_left - bits
end do L2
end do L1
res(i) = loc + nn
end do
end function rvs_unif_array_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in REAL_KINDS_TYPES
function rvs_unif_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res)
integer, intent(in) :: array_size
${t1}$, intent(in) :: loc, scale
${t1}$ :: res(array_size)
${t1}$ :: t
integer(int64) :: tmp
integer :: i
if(scale == 0._${k1}$) call error_stop("Error(rvs_unif_array):" &
//" Uniform distribution scale parameter must be non-zero")
do i = 1, array_size
tmp = shiftr(dist_rand(INT_ONE), 11)
t = real(tmp * MESENNE_NUMBER, kind = ${k1}$)
res(i) = loc + scale * t
end do
end function rvs_unif_array_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in CMPLX_KINDS_TYPES
function rvs_unif_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res)
integer, intent(in) :: array_size
${t1}$, intent(in) :: loc, scale
${t1}$ :: res(array_size)
real(${k1}$) :: r1, tr, ti
integer(int64) :: tmp
integer :: i
if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(rvs_unif" &
//"_array): Uniform distribution scale parameter must be non-zero")
do i = 1, array_size
tmp = shiftr(dist_rand(INT_ONE), 11)
r1 = real(tmp * MESENNE_NUMBER, kind = ${k1}$)
if(scale % re == 0.0_${k1}$) then
tr = loc % re
ti = loc % im + scale % im * r1
else if(scale % im == 0.0_${k1}$) then
tr = loc % re + scale % re * r1
ti = loc % im
else
tr = loc % re + scale % re * r1
tmp = shiftr(dist_rand(INT_ONE), 11)
r1 = real(tmp * MESENNE_NUMBER, kind = ${k1}$)
ti = loc % im + scale % im * r1
end if
res(i) = cmplx(tr, ti, kind=${k1}$)
end do
end function rvs_unif_array_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
elemental function pdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res)
${t1}$, intent(in) :: x, loc, scale
real :: res
if(scale == 0_${k1}$) then
res = 0.0
else if(x < loc .or. x > (loc + scale)) then
res = 0.0
else
res = 1. / (scale + 1_${k1}$)
end if
end function pdf_unif_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in REAL_KINDS_TYPES
elemental function pdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res)
${t1}$, intent(in) :: x, loc, scale
${t1}$ :: res
${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$
if(scale == zero) then
res = zero
else if(x < loc .or. x > (loc + scale)) then
res = zero
else
res = one / scale
end if
end function pdf_unif_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in CMPLX_KINDS_TYPES
elemental function pdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res)
${t1}$, intent(in) :: x, loc, scale
real(${k1}$) :: res, tr, ti
real(${k1}$), parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$
tr = loc % re + scale % re; ti = loc % im + scale % im
if(scale == (zero, zero)) then
res = zero
else if((x % re >= loc % re .and. x % re <= tr) .and. &
(x % im >= loc % im .and. x % im <= ti)) then
res = one / (scale % re * scale % im)
else
res = zero
end if
end function pdf_unif_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in INT_KINDS_TYPES
elemental function cdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res)
${t1}$, intent(in) :: x, loc, scale
real :: res
if(scale == 0_${k1}$) then
res = 0.0
else if(x < loc) then
res = 0.0
else if(x >= loc .and. x <= (loc + scale)) then
res = real((x - loc + 1_${k1}$)) / real((scale + 1_${k1}$))
else
res = 1.0
end if
end function cdf_unif_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in REAL_KINDS_TYPES
elemental function cdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res)
${t1}$, intent(in) :: x, loc, scale
${t1}$ :: res
${t1}$, parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$
if(scale == zero) then
res = zero
else if(x < loc) then
res = zero
else if(x >= loc .and. x <= (loc + scale)) then
res = (x - loc) / scale
else
res = one
end if
end function cdf_unif_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in CMPLX_KINDS_TYPES
elemental function cdf_unif_${t1[0]}$${k1}$(x, loc, scale) result(res)
${t1}$, intent(in) :: x, loc, scale
real(${k1}$) :: res
real(${k1}$), parameter :: zero = 0.0_${k1}$, one = 1.0_${k1}$
logical :: r1, r2, i1, i2
if(scale == (zero, zero)) then
res = zero
return
end if
r1 = x % re < loc % re
r2 = x % re > (loc % re + scale % re)
i1 = x % im < loc % im
i2 = x % im > (loc % im + scale % im)
if(r1 .or. i1) then
res = zero
else if((.not. r1) .and. (.not. r2) .and. i2) then
res = (x % re - loc % re) / scale % re
else if((.not. i1) .and. (.not. i2) .and. r2) then
res = (x % im - loc % im) / scale % im
else if((.not. r1) .and. (.not. r2) .and. (.not. i1) .and. (.not. i2)) &
then
res = ((x % re - loc % re) / scale % re) * ((x % im - loc % im) / &
scale % im)
else if(r2 .and. i2)then
res = one
end if
end function cdf_unif_${t1[0]}$${k1}$
#:endfor
#:for k1, t1 in ALL_KINDS_TYPES
function shuffle_${t1[0]}$${k1}$( list ) result(res)
${t1}$, intent(in) :: list(:)
${t1}$ :: res(size(list))
${t1}$ :: tmp
integer :: n, i, j
n = size(list)
res = list
do i = 1, n - 1
j = rvs_uniform(n - i) + i
tmp = res(i)
res(i) = res(j)
res(j) = tmp
end do
end function shuffle_${t1[0]}$${k1}$
#:endfor
end module stdlib_stats_distribution_uniform
