Probability Distribution and Statistical Functions -- PRNG Module

doc/specs/index.md

@@ -19,6 +19,7 @@ This is and index/directory of the specifications (specs) for each new module/fe
  - [optval](./stdlib_optval.html) - Fallback value for optional arguments
  - [quadrature](./stdlib_quadrature.html) - Numerical integration
  - [stats](./stdlib_stats.html) - Descriptive Statistics
+ - [stats_distribution_PRNG](./stdlib_stats_distribution_PRNG.html) - Probability Distributions random number generator
 
 ## Missing specs
 

doc/specs/stdlib_stats_distribution_PRNG.md

@@ -0,0 +1,87 @@
+---
+title: stats_distribution
+---
+
+# Statistical Distributions -- Pseudorandom Number Generator Module
+
+[TOC]
+
+## `random_seed` - set or get a value of seed to the probability distribution pseudorandom number generator
+
+### Status
+
+Experimental
+
+### Description
+
+Set or get the seed value before calling the probability distribution pseudorandom number generator for variates.
+
+### Syntax
+
+`call [[stdlib_stats_distribution_PRNG(module):random_seed(interface)]](put, get)`
+
+### Arguments
+
+`put`: argument has `intent(in)` and may be a scalar of type `integer`.
+
+`get`: argument has `intent(out)` and is a scalar of type `integer`.
+
+### Return value
+
+Return a scalar of type `integer`.
+
+### Example
+
+```fortran
+program demo_random_seed
+    use stdlib_stats_distribution_PRNG, only : random_seed
+    implicit none
+    integer :: seed_put, seed_get
+
+    seed_put = 1234567
+    call random_seed(seed_put, seed_get)     ! set and get current value of seed
+end program demo_random_seed
+```
+
+## `dist_rand` - Get a random integer with specified kind
+
+### Status
+
+Experimental
+
+### Description
+
+Generate an integer pseudorandom number in a specific range [-2^k, 2^k - 1] according to the input integer kind n. This pseudorandom number will be operated by bit opeartors instead of normal arithmetic operators.
+
+### Syntax
+
+`result = [[stdlib_stats_distribution_PRNG(module):dist_rand(interface)]](n)`
+
+### Arguments
+
+`n`: argument has `intent(in)` is a scalar of type `integer`.
+
+### Return value
+
+Return a scalar of type `integer`.
+
+### Example
+
+```fortran
+program demo_dist_rand
+    use stdlib_stats_distribution_PRNG, only : dist_rand, random_seed
+    implicit none
+    integer :: put, get
+
+    put = 135792468
+    call random_seed(put, get)     ! set and get current value of seed
+    print *, dist_rand(1_int8)     ! random integer in [-2^7, 2^7 - 1]
+! -90
+    print *, dist_rand(1_int16)    ! random integer in [-2^15, 2^15 - 1]
+! -32725
+    print *, dist_rand(1_int32)    ! random integer in [-2^31, 2^31 - 1]
+! -1601563881
+    print *, dist_rand(1_int64)    ! random integer in [-2^63, 2^63 - 1]
+! 180977695517992208
+end program demo_dist_rand
+```

src/CMakeLists.txt

@@ -21,6 +21,7 @@ set(fppFiles
     stdlib_quadrature.fypp
     stdlib_quadrature_trapz.fypp
     stdlib_quadrature_simps.fypp
+    stdlib_stats_distribution_PRNG.fypp
 )
 
 

src/Makefile.manual

@@ -17,7 +17,8 @@ SRCFYPP =\
         stdlib_stats_moment_all.fypp \
         stdlib_stats_moment_mask.fypp \
         stdlib_stats_moment_scalar.fypp \
-        stdlib_stats_var.fypp	
+        stdlib_stats_var.fypp \
+	stdlib_stats_distribution_PRNG.fypp
 
 SRC = f18estop.f90 \
       stdlib_ascii.f90 \
@@ -104,3 +105,6 @@ stdlib_stats_var.o: \
 	stdlib_optval.o \
 	stdlib_kinds.o \
 	stdlib_stats.o
+stdlib_stats_distribution_PRNG.o: \
+        stdlib_kinds.o \
+	stdlib_error.o

src/stdlib_stats_distribution_PRNG.fypp

@@ -0,0 +1,151 @@
+#:include "common.fypp"
+module stdlib_stats_distribution_PRNG
+    use stdlib_kinds, only: int8, int16, int32, int64
+    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_stats_distribution_PRNG.html#
+    !! description))
+        #: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_stats_distribution_PRNG.html#
+    !! description))
+    !!
+        #: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 ([email protected])
+    ! 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 = rol64(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) = rol64(st(4), 45)
+    end function xoshiro256ss
+
+    pure function rol64(x, k) result(res)
+        integer(int64), intent(in) :: x
+        integer, intent(in) :: k
+        integer(int64) :: t1, t2, res
+
+        t1 = shiftr(x, (64 - k))
+        t2 = shiftl(x, k)
+        res = ior(t1, t2)
+    end function rol64
+
+
+    function splitmix64(s) result(res)
+    ! Written in 2015 by Sebastiano Vigna ([email protected])
+    ! 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
+
+        if(present(s)) si = s
+        res = 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_stats_distribution_PRNG

src/tests/stats/CMakeLists.txt

@@ -5,6 +5,7 @@ ADDTEST(moment)
 ADDTEST(rawmoment)
 ADDTEST(var)
 ADDTEST(varn)
+ADDTEST(distribution_PRNG)
 
 if(DEFINED CMAKE_MAXIMUM_RANK)
     if(${CMAKE_MAXIMUM_RANK} GREATER 7)

src/tests/stats/Makefile.manual

@@ -1,3 +1,5 @@
-PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90
 
-include ../Makefile.manual.test.mk
+PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90 \
+            test_distribution_PRNG.f90
+
+include ../Makefile.manual.test.mk

src/tests/stats/test_distribution_PRNG.f90

@@ -0,0 +1,108 @@
+program test_distribution_PRNG
+    use stdlib_error, only : check
+    use stdlib_kinds, only: int8, int16, int32, int64
+    use stdlib_stats_distribution_PRNG, only : random_seed, dist_rand
+
+    implicit none
+    logical ::  warn = .true.
+
+    call test_random_seed
+    call test_random_rand_iint8
+    call test_random_rand_iint16
+    call test_random_rand_iint32
+    call test_random_rand_iint64
+
+
+    contains
+
+    subroutine test_random_seed
+        integer :: put, get, res(5)
+        integer :: ans(5) = [-1859553078, -1933696596, -642834430, &
+                             1711399314,  1548311463]
+        integer :: i
+
+        print *, ""
+        print *, "Test random_seed"
+        put = 135792468
+        do i = 1, 5
+            call random_seed(put,get)
+            res(i) = get
+            put = get
+        end do
+        call check(all(res == ans), msg="random seed test failed.",warn=warn)
+    end subroutine test_random_seed
+
+    subroutine test_random_rand_iint8
+        integer :: put, get, i
+
+        integer(int8) :: res(5), ans(5)=[118, -15, -72, 101, 70]
+
+
+        print *, ""
+        print *, "Test random_rand with kind int8"
+        put = 12345678
+        call random_seed(put, get)
+        do i = 1, 5
+            res(i) = dist_rand(1_int8)
+        end do
+        call check(all(res == ans), msg="random_rand with kind int8 test" &
+                                     //" failed.", warn=warn)
+    end subroutine test_random_rand_iint8
+
+    subroutine test_random_rand_iint16
+        integer :: put, get, i
+
+        integer(int16) :: res(5), ans(5)=[30286, -3799, -18204, 25947, 18148]
+
+
+        print *, ""
+        print *, "Test random_rand with kind int16"
+        put = 12345678
+        call random_seed(put, get)
+        do i = 1, 5
+            res(i) = dist_rand(1_int16)
+        end do
+        call check(all(res == ans), msg="random_rand with kind int16 test" &
+                                     //" failed.", warn=warn)
+    end subroutine test_random_rand_iint16
+
+    subroutine test_random_rand_iint32
+        integer :: put, get, i
+
+        integer(int32) :: res(5), ans(5)=[1984865646, -248954393, -1192993267, &
+                                  1700514835, 1189401802]
+
+
+        print *, ""
+        print *, "Test random_rand with kind int32"
+        put = 12345678
+        call random_seed(put, get)
+        do i = 1, 5
+            res(i) = dist_rand(1_int32)
+        end do
+        call check(all(res == ans), msg="random_rand with kind int32 test" &
+                                     //" failed.", warn=warn)
+    end subroutine test_random_rand_iint32
+
+    subroutine test_random_rand_iint64
+        integer :: put, get, i
+
+        integer(int64) :: res(5), ans(5)=[8524933037632333570_int64, &
+                                  -1069250973542918798_int64,        &
+                                  -5123867065024149335_int64,        &
+                                  7303655603304982073_int64,         &
+                                  5108441843522503546_int64]
+
+
+        print *, ""
+        print *, "Test random_rand with kind int64"
+        put = 12345678
+        call random_seed(put, get)
+        do i = 1, 5
+            res(i) = dist_rand(1_int64)
+        end do
+        call check(all(res == ans), msg="random_rand with kind int64 test" &
+                                     //" failed.", warn=warn)
+    end subroutine test_random_rand_iint64
+
+end program test_distribution_PRNG