Skip to content

Commit 0beff2a

Browse files
Merge pull request fortran-lang#20 from fortran-lang/master
merge
2 parents ba70104 + a1c911c commit 0beff2a

9 files changed

+364
-3
lines changed

README.md

+6
Original file line numberDiff line numberDiff line change
@@ -82,6 +82,12 @@ GCC Fortran (MSYS) | 10 | Windows Server 2019 | x86_64
8282
GCC Fortran (MinGW) | 10 | Windows Server 2019 | x86_64, i686
8383
Intel oneAPI classic | 2021.1 | Ubuntu 20.04 | x86_64
8484

85+
The following combinations are known to work, but they are not tested in the CI:
86+
87+
Name | Version | Platform | Architecture
88+
--- | --- | --- | ---
89+
GCC Fortran (MinGW) | 8.4.0, 9.3.0, 10.2.0 | Windows 10 | x86_64, i686
90+
8591
We try to test as many available compilers and platforms as possible.
8692
A list of tested compilers which are currently not working and the respective issue are listed below.
8793

doc/specs/index.md

+1
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@ This is and index/directory of the specifications (specs) for each new module/fe
1919
- [optval](./stdlib_optval.html) - Fallback value for optional arguments
2020
- [quadrature](./stdlib_quadrature.html) - Numerical integration
2121
- [stats](./stdlib_stats.html) - Descriptive Statistics
22+
- [stats_distribution_PRNG](./stdlib_stats_distribution_PRNG.html) - Probability Distributions random number generator
2223

2324
## Missing specs
2425

Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
---
2+
title: stats_distribution
3+
---
4+
5+
# Statistical Distributions -- Pseudorandom Number Generator Module
6+
7+
[TOC]
8+
9+
## `random_seed` - set or get a value of seed to the probability distribution pseudorandom number generator
10+
11+
### Status
12+
13+
Experimental
14+
15+
### Description
16+
17+
Set or get the seed value before calling the probability distribution pseudorandom number generator for variates.
18+
19+
### Syntax
20+
21+
`call [[stdlib_stats_distribution_PRNG(module):random_seed(interface)]](put, get)`
22+
23+
### Arguments
24+
25+
`put`: argument has `intent(in)` and may be a scalar of type `integer`.
26+
27+
`get`: argument has `intent(out)` and is a scalar of type `integer`.
28+
29+
### Return value
30+
31+
Return a scalar of type `integer`.
32+
33+
### Example
34+
35+
```fortran
36+
program demo_random_seed
37+
use stdlib_stats_distribution_PRNG, only : random_seed
38+
implicit none
39+
integer :: seed_put, seed_get
40+
41+
seed_put = 1234567
42+
call random_seed(seed_put, seed_get) ! set and get current value of seed
43+
end program demo_random_seed
44+
```
45+
46+
## `dist_rand` - Get a random integer with specified kind
47+
48+
### Status
49+
50+
Experimental
51+
52+
### Description
53+
54+
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.
55+
56+
### Syntax
57+
58+
`result = [[stdlib_stats_distribution_PRNG(module):dist_rand(interface)]](n)`
59+
60+
### Arguments
61+
62+
`n`: argument has `intent(in)` is a scalar of type `integer`.
63+
64+
### Return value
65+
66+
Return a scalar of type `integer`.
67+
68+
### Example
69+
70+
```fortran
71+
program demo_dist_rand
72+
use stdlib_stats_distribution_PRNG, only : dist_rand, random_seed
73+
implicit none
74+
integer :: put, get
75+
76+
put = 135792468
77+
call random_seed(put, get) ! set and get current value of seed
78+
print *, dist_rand(1_int8) ! random integer in [-2^7, 2^7 - 1]
79+
! -90
80+
print *, dist_rand(1_int16) ! random integer in [-2^15, 2^15 - 1]
81+
! -32725
82+
print *, dist_rand(1_int32) ! random integer in [-2^31, 2^31 - 1]
83+
! -1601563881
84+
print *, dist_rand(1_int64) ! random integer in [-2^63, 2^63 - 1]
85+
! 180977695517992208
86+
end program demo_dist_rand
87+
```

src/CMakeLists.txt

+1
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@ set(fppFiles
2121
stdlib_quadrature.fypp
2222
stdlib_quadrature_trapz.fypp
2323
stdlib_quadrature_simps.fypp
24+
stdlib_stats_distribution_PRNG.fypp
2425
)
2526

2627

src/Makefile.manual

+5-1
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,8 @@ SRCFYPP =\
1717
stdlib_stats_moment_all.fypp \
1818
stdlib_stats_moment_mask.fypp \
1919
stdlib_stats_moment_scalar.fypp \
20-
stdlib_stats_var.fypp
20+
stdlib_stats_var.fypp \
21+
stdlib_stats_distribution_PRNG.fypp
2122

2223
SRC = f18estop.f90 \
2324
stdlib_ascii.f90 \
@@ -104,3 +105,6 @@ stdlib_stats_var.o: \
104105
stdlib_optval.o \
105106
stdlib_kinds.o \
106107
stdlib_stats.o
108+
stdlib_stats_distribution_PRNG.o: \
109+
stdlib_kinds.o \
110+
stdlib_error.o
+151
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,151 @@
1+
#:include "common.fypp"
2+
module stdlib_stats_distribution_PRNG
3+
use stdlib_kinds, only: int8, int16, int32, int64
4+
use stdlib_error, only: error_stop
5+
implicit none
6+
private
7+
integer, parameter :: MAX_INT_BIT_SIZE = bit_size(1_int64)
8+
integer(int64) :: st(4) ! internal states for xoshiro256ss function
9+
integer(int64) :: si = 614872703977525537_int64 ! default seed value
10+
logical :: seed_initialized = .false.
11+
12+
public :: random_seed
13+
public :: dist_rand
14+
15+
16+
interface dist_rand
17+
!! Version experimental
18+
!!
19+
!! Generation of random integers with different kinds
20+
!! ([Specification](../page/specs/stdlib_stats_distribution_PRNG.html#
21+
!! description))
22+
#:for k1, t1 in INT_KINDS_TYPES
23+
module procedure dist_rand_${t1[0]}$${k1}$
24+
#:endfor
25+
end interface dist_rand
26+
27+
interface random_seed
28+
!! Version experimental
29+
!!
30+
!! Set seed value for random number generator
31+
!! ([Specification](../page/specs/stdlib_stats_distribution_PRNG.html#
32+
!! description))
33+
!!
34+
#:for k1, t1 in INT_KINDS_TYPES
35+
module procedure random_distribution_seed_${t1[0]}$${k1}$
36+
#:endfor
37+
end interface random_seed
38+
39+
40+
contains
41+
42+
#:for k1, t1 in INT_KINDS_TYPES
43+
function dist_rand_${t1[0]}$${k1}$(n) result(res)
44+
!! Random integer generation for various kinds
45+
!! result = [-2^k, 2^k - 1], k = 7, 15, 31, 63, depending on input kind
46+
!! Result will be operated by bitwise operators to generate desired integer
47+
!! and real pseudorandom numbers
48+
!!
49+
${t1}$, intent(in) :: n
50+
${t1}$ :: res
51+
integer :: k
52+
53+
k = MAX_INT_BIT_SIZE - bit_size(n)
54+
if(k < 0) call error_stop("Error(dist_rand): Integer bit size is" &
55+
//" greater than 64bit")
56+
res = shiftr(xoshiro256ss( ), k)
57+
end function dist_rand_${t1[0]}$${k1}$
58+
59+
#:endfor
60+
61+
function xoshiro256ss( ) result (res)
62+
! Generate random 64-bit integers
63+
! Written in 2018 by David Blackman and Sebastiano Vigna ([email protected])
64+
! http://prng.di.unimi.it/xoshiro256starstar.c
65+
!
66+
! This is xoshiro256** 1.0, one of our all-purpose, rock-solid
67+
! generators. It has excellent (sub-ns) speed, a state (256 bits) that is
68+
! large enough for any parallel application, and it passes all tests we
69+
! are aware of.
70+
!
71+
! The state must be seeded so that it is not everywhere zero. If you have
72+
! a 64-bit seed, we suggest to seed a splitmix64 generator and use its
73+
! output to fill st.
74+
!
75+
! Fortran 90 version translated from C by Jim-215-Fisher
76+
!
77+
integer(int64) :: res, t
78+
79+
if(.not. seed_initialized) call random_distribution_seed_iint64(si,t)
80+
res = rol64(st(2) * 5, 7) * 9
81+
t = shiftl(st(2), 17)
82+
st(3) = ieor(st(3), st(1))
83+
st(4) = ieor(st(4), st(2))
84+
st(2) = ieor(st(2), st(3))
85+
st(1) = ieor(st(1), st(4))
86+
st(3) = ieor(st(3), t)
87+
st(4) = rol64(st(4), 45)
88+
end function xoshiro256ss
89+
90+
pure function rol64(x, k) result(res)
91+
integer(int64), intent(in) :: x
92+
integer, intent(in) :: k
93+
integer(int64) :: t1, t2, res
94+
95+
t1 = shiftr(x, (64 - k))
96+
t2 = shiftl(x, k)
97+
res = ior(t1, t2)
98+
end function rol64
99+
100+
101+
function splitmix64(s) result(res)
102+
! Written in 2015 by Sebastiano Vigna ([email protected])
103+
! This is a fixed-increment version of Java 8's SplittableRandom
104+
! generator.
105+
! See http://dx.doi.org/10.1145/2714064.2660195 and
106+
! http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html
107+
!
108+
! It is a very fast generator passing BigCrush, and it can be useful if
109+
! for some reason you absolutely want 64 bits of state.
110+
!
111+
! Fortran 90 translated from C by Jim-215-Fisher
112+
!
113+
integer(int64) :: res
114+
integer(int64), intent(in), optional :: s
115+
integer(int64) :: int01 = -7046029254386353131_int64, &
116+
int02 = -4658895280553007687_int64, &
117+
int03 = -7723592293110705685_int64
118+
! Values are converted from C unsigned integer of 0x9e3779b97f4a7c15,
119+
! 0xbf58476d1ce4e5b9, 0x94d049bb133111eb
120+
121+
if(present(s)) si = s
122+
res = si
123+
si = res + int01
124+
res = ieor(res, shiftr(res, 30)) * int02
125+
res = ieor(res, shiftr(res, 27)) * int03
126+
res = ieor(res, shiftr(res, 31))
127+
end function splitmix64
128+
129+
#:for k1, t1 in INT_KINDS_TYPES
130+
subroutine random_distribution_seed_${t1[0]}$${k1}$(put, get)
131+
!! Set seed value for random number generator
132+
!!
133+
${t1}$, intent(in) :: put
134+
${t1}$, intent(out) :: get
135+
integer(int64) :: tmp
136+
integer :: i
137+
138+
tmp = splitmix64(int(put, kind = int64))
139+
do i = 1, 10
140+
tmp = splitmix64( )
141+
end do
142+
do i = 1, 4
143+
tmp = splitmix64( )
144+
st(i) = tmp
145+
end do
146+
get = int(tmp, kind = ${k1}$)
147+
seed_initialized = .true.
148+
end subroutine random_distribution_seed_${t1[0]}$${k1}$
149+
150+
#:endfor
151+
end module stdlib_stats_distribution_PRNG

src/tests/stats/CMakeLists.txt

+1
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ ADDTEST(moment)
55
ADDTEST(rawmoment)
66
ADDTEST(var)
77
ADDTEST(varn)
8+
ADDTEST(distribution_PRNG)
89

910
if(DEFINED CMAKE_MAXIMUM_RANK)
1011
if(${CMAKE_MAXIMUM_RANK} GREATER 7)

src/tests/stats/Makefile.manual

+4-2
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,5 @@
1-
PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90
21

3-
include ../Makefile.manual.test.mk
2+
PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90 \
3+
test_distribution_PRNG.f90
4+
5+
include ../Makefile.manual.test.mk

0 commit comments

Comments
 (0)