diff --git a/doc/specs/index.md b/doc/specs/index.md index 6ea78b52e..edd00ed95 100644 --- a/doc/specs/index.md +++ b/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 diff --git a/doc/specs/stdlib_stats_distribution_PRNG.md b/doc/specs/stdlib_stats_distribution_PRNG.md new file mode 100644 index 000000000..60067bde6 --- /dev/null +++ b/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 +``` diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 02413b415..2429f555a 100644 --- a/src/CMakeLists.txt +++ b/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 ) diff --git a/src/Makefile.manual b/src/Makefile.manual index 9351a374a..a57253e2b 100644 --- a/src/Makefile.manual +++ b/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 diff --git a/src/stdlib_stats_distribution_PRNG.fypp b/src/stdlib_stats_distribution_PRNG.fypp new file mode 100644 index 000000000..60448c778 --- /dev/null +++ b/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 (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 = 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 (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 + + 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 diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index 36ffc7aeb..38f9bb84b 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/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) diff --git a/src/tests/stats/Makefile.manual b/src/tests/stats/Makefile.manual index aacaf98ab..a0731e8f4 100644 --- a/src/tests/stats/Makefile.manual +++ b/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 \ No newline at end of file diff --git a/src/tests/stats/test_distribution_PRNG.f90 b/src/tests/stats/test_distribution_PRNG.f90 new file mode 100644 index 000000000..923091dd3 --- /dev/null +++ b/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