stdlib_random.fypp Source File

#:include "common.fypp"
module stdlib_random
    use stdlib_kinds, only: int8, int16, int32, int64
    use stdlib_optval, only: optval
    use stdlib_error, only: error_stop
    implicit none
    private
    integer, parameter :: MAX_INT_BIT_SIZE = bit_size(1_int64)
    integer(int64) :: st(4) ! internal states for xoshiro256ss function
    integer(int64) :: si = 614872703977525537_int64 ! default seed value
    logical :: seed_initialized = .false.

    public :: random_seed
    public :: dist_rand


    interface dist_rand
    !! Version experimental
    !!
    !! Generation of random integers with different kinds
    !! ([Specification](../page/specs/stdlib_random.html#
    !! dist_rand-get-a-random-integer-with-specified-kind))
        #:for k1, t1 in INT_KINDS_TYPES
        module procedure dist_rand_${t1[0]}$${k1}$
        #:endfor
    end interface dist_rand

    interface random_seed
    !! Version experimental
    !!
    !! Set seed value for random number generator
    !! ([Specification](../page/specs/stdlib_random.html#
    !! random_seed-set-or-get-a-value-of-seed-to-the-probability-distribution-pseudorandom-number-generator))
    !!
        #:for k1, t1 in INT_KINDS_TYPES
        module procedure random_distribution_seed_${t1[0]}$${k1}$
        #:endfor
    end interface random_seed


    contains

    #:for k1, t1 in INT_KINDS_TYPES
    function dist_rand_${t1[0]}$${k1}$(n) result(res)
    !! Random integer generation for various kinds
    !! result = [-2^k, 2^k - 1], k = 7, 15, 31, 63, depending on input kind
    !! Result will be operated by bitwise operators to generate desired integer
    !! and real pseudorandom numbers
    !!
        ${t1}$, intent(in) :: n
        ${t1}$ :: res
        integer :: k

        k = MAX_INT_BIT_SIZE - bit_size(n)
        if(k < 0) call error_stop("Error(dist_rand): Integer bit size is"       &
                  //" greater than 64bit")
        res = shiftr(xoshiro256ss( ), k)
    end function dist_rand_${t1[0]}$${k1}$

    #:endfor

    function xoshiro256ss( ) result (res)
    ! Generate random 64-bit integers
    ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
    ! http://prng.di.unimi.it/xoshiro256starstar.c
    !
    ! This is xoshiro256** 1.0, one of our all-purpose, rock-solid
    ! generators. It has excellent (sub-ns) speed, a state (256 bits) that is
    ! large enough for any parallel application, and it passes all tests we
    ! are aware of.
    !
    ! The state must be seeded so that it is not everywhere zero. If you have
    ! a 64-bit seed, we suggest to seed a splitmix64 generator and use its
    ! output to fill st.
    !
    ! Fortran 90 version translated from C by Jim-215-Fisher
    !
        integer(int64) :: res, t

        if(.not. seed_initialized) call random_distribution_seed_iint64(si,t)
        res = ishftc(st(2) * 5, 7) * 9
        t = shiftl(st(2), 17)
        st(3) = ieor(st(3), st(1))
        st(4) = ieor(st(4), st(2))
        st(2) = ieor(st(2), st(3))
        st(1) = ieor(st(1), st(4))
        st(3) = ieor(st(3), t)
        st(4) = ishftc(st(4), 45)
    end function xoshiro256ss

    function splitmix64(s) result(res)
    ! Written in 2015 by Sebastiano Vigna (vigna@acm.org)
    ! This is a fixed-increment version of Java 8's SplittableRandom
    ! generator.
    ! See http://dx.doi.org/10.1145/2714064.2660195 and
    ! http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html
    !
    ! It is a very fast generator passing BigCrush, and it can be useful if
    ! for some reason you absolutely want 64 bits of state.
    !
    ! Fortran 90 translated from C by Jim-215-Fisher
    !
        integer(int64) :: res
        integer(int64), intent(in), optional :: s
        integer(int64) :: int01 = -7046029254386353131_int64,                  &
                          int02 = -4658895280553007687_int64,                  &
                          int03 = -7723592293110705685_int64
    ! Values are converted from C unsigned integer of 0x9e3779b97f4a7c15,
    ! 0xbf58476d1ce4e5b9, 0x94d049bb133111eb

        res = optval(s, si)
        si = res + int01
        res = ieor(res, shiftr(res, 30)) * int02
        res = ieor(res, shiftr(res, 27)) * int03
        res = ieor(res, shiftr(res, 31))
    end function splitmix64

    #:for k1, t1 in INT_KINDS_TYPES
    subroutine random_distribution_seed_${t1[0]}$${k1}$(put, get)
    !! Set seed value for random number generator
    !!
        ${t1}$, intent(in) :: put
        ${t1}$, intent(out) :: get
        integer(int64) :: tmp
        integer :: i

        tmp = splitmix64(int(put, kind = int64))
        do i = 1, 10
            tmp = splitmix64( )
        end do
        do i = 1, 4
            tmp = splitmix64( )
            st(i) = tmp
        end do
        get = int(tmp, kind = ${k1}$)
        seed_initialized = .true.
    end subroutine random_distribution_seed_${t1[0]}$${k1}$

    #:endfor
end module stdlib_random