stdlib_stats_distribution_uniform.fypp Source File


Contents


Source Code

#: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