Skip to content

Commit e026ab7

Browse files
Jim-215-FisherJim-215-Fisher
authored
initial commit
1 parent e4f3de1 commit e026ab7

File tree

4 files changed

+1658
-0
lines changed

4 files changed

+1658
-0
lines changed

‎src/stdlib_stats_distribution_PRNG.fypp‎

Lines changed: 212 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,212 @@
1+
#:include "common.fypp"
2+
module stdlib_stats_distribution_PRNG
3+
use stdlib_kinds, only: int8, int16, int32, int64
4+
implicit none
5+
private
6+
integer, parameter :: MAX_INT_BIT_SIZE = bit_size(1_int64)
7+
integer(int64), save :: st(4), si = 614872703977525537_int64
8+
logical, save :: seed_initialized = .false.
9+
10+
public :: random_seed
11+
public :: dist_rand
12+
public :: jump
13+
public :: long_jump
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+
res = shiftr(xoshiro256ss( ), k)
55+
end function dist_rand_${t1[0]}$${k1}$
56+
57+
#:endfor
58+
59+
function xoshiro256ss( ) result (res)
60+
! Generate random 64-bit integers
61+
! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
62+
! http://prng.di.unimi.it/xoshiro256starstar.c
63+
!
64+
! This is xoshiro256** 1.0, one of our all-purpose, rock-solid
65+
! generators. It has excellent (sub-ns) speed, a state (256 bits) that is
66+
! large enough for any parallel application, and it passes all tests we
67+
! are aware of.
68+
!
69+
! The state must be seeded so that it is not everywhere zero. If you have
70+
! a 64-bit seed, we suggest to seed a splitmix64 generator and use its
71+
! output to fill st.
72+
!
73+
! Fortran 90 version translated from C by Jim-215-Fisher
74+
!
75+
integer(int64) :: res, t
76+
77+
if(.not. seed_initialized) call random_distribution_seed_iint64(si,t)
78+
res = rol64(st(2) * 5 , 7) * 9
79+
t = shiftl(st(2), 17)
80+
st(3) = ieor(st(3), st(1))
81+
st(4) = ieor(st(4), st(2))
82+
st(2) = ieor(st(2), st(3))
83+
st(1) = ieor(st(1), st(4))
84+
st(3) = ieor(st(3), t)
85+
st(4) = rol64(st(4), 45)
86+
end function xoshiro256ss
87+
88+
function rol64(x, k) result(res)
89+
integer(int64), intent(in) :: x
90+
integer, intent(in) :: k
91+
integer(int64) :: t1, t2, res
92+
93+
t1 = shiftr(x, (64 - k))
94+
t2 = shiftl(x, k)
95+
res = ior(t1, t2)
96+
end function rol64
97+
98+
99+
subroutine jump
100+
! This is the jump function for the xoshiro256ss generator. It is equivalent
101+
! to 2^128 calls to xoshiro256ss(); it can be used to generate 2^128
102+
! non-overlapping subsequences for parallel computations.
103+
! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
104+
! http://prng.di.unimi.it/xoshiro256starstar.c
105+
!
106+
! Fortran 90 version translated from C by Jim-215-Fisher
107+
integer(int64) :: jp(4) = [1733541517147835066_int64, &
108+
-3051731464161248980_int64, &
109+
-6244198995065845334_int64, &
110+
4155657270789760540_int64]
111+
integer(int64) :: s1 = 0, s2 = 0, s3 = 0, s4 = 0, c = 1_int64
112+
integer :: i, j, k
113+
114+
do i = 1, 4
115+
do j = 1, 64
116+
if(iand(jp(i), shiftl(c, j - 1)) /= 0) then
117+
s1 = ieor(s1, st(1))
118+
s2 = ieor(s2, st(2))
119+
s3 = ieor(s3, st(3))
120+
s4 = ieor(s4, st(4))
121+
end if
122+
k = xoshiro256ss( )
123+
end do
124+
end do
125+
st(1) = s1
126+
st(2) = s2
127+
st(3) = s3
128+
st(4) = s4
129+
end subroutine jump
130+
131+
subroutine long_jump
132+
! This is the long-jump function for the xoshiro256ss generator. It is
133+
! equivalent to 2^192 calls to xoshiro256ss(); it can be used to generate
134+
! 2^64 starting points, from each of which jump() will generate 2^64
135+
! non-overlapping subsequences for parallel distributed computations
136+
! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
137+
! http://prng.di.unimi.it/xoshiro256starstar.c
138+
!
139+
! Fortran 90 version translated from C by Jim-215-Fisher
140+
integer(int64) :: jp(4) = [8566230491382795199_int64, &
141+
-4251311993797857357_int64, &
142+
8606660816089834049_int64, &
143+
4111957640723818037_int64]
144+
integer(int64) :: s1 = 0, s2 = 0, s3 = 0, s4 = 0, c = 1_int64
145+
integer(int32) :: i, j, k
146+
147+
do i = 1, 4
148+
do j = 1, 64
149+
if(iand(jp(i), shiftl(c, j - 1)) /= 0) then
150+
s1 = ieor(s1, st(1))
151+
s2 = ieor(s2, st(2))
152+
s3 = ieor(s3, st(3))
153+
s4 = ieor(s4, st(4))
154+
end if
155+
k = xoshiro256ss()
156+
end do
157+
end do
158+
st(1) = s1
159+
st(2) = s2
160+
st(3) = s3
161+
st(4) = s4
162+
end subroutine long_jump
163+
164+
function splitmix64(s) result(res)
165+
! Written in 2015 by Sebastiano Vigna (vigna@acm.org)
166+
! This is a fixed-increment version of Java 8's SplittableRandom
167+
! generator.
168+
! See http://dx.doi.org/10.1145/2714064.2660195 and
169+
! http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html
170+
!
171+
! It is a very fast generator passing BigCrush, and it can be useful if
172+
! for some reason you absolutely want 64 bits of state.
173+
!
174+
! Fortran 90 translated from C by Jim-215-Fisher
175+
!
176+
integer(int64) :: res, int01, int02, int03
177+
integer(int64), intent(in), optional :: s
178+
data int01, int02, int03/-7046029254386353131_int64, &
179+
-4658895280553007687_int64, &
180+
-7723592293110705685_int64/
181+
182+
if(present(s)) si = s
183+
res = si
184+
si = res + int01
185+
res = ieor(res, shiftr(res, 30)) * int02
186+
res = ieor(res, shiftr(res, 27)) * int03
187+
res = ieor(res, shiftr(res, 31))
188+
end function splitmix64
189+
190+
#:for k1, t1 in INT_KINDS_TYPES
191+
subroutine random_distribution_seed_${t1[0]}$${k1}$(put, get)
192+
!! Set seed value for random number generator
193+
!!
194+
${t1}$, intent(in) :: put
195+
${t1}$, intent(out) :: get
196+
integer(int64) :: tmp
197+
integer :: i
198+
199+
tmp = splitmix64(int(put, kind = int64))
200+
do i = 1, 10
201+
tmp = splitmix64( )
202+
end do
203+
do i = 1, 4
204+
tmp = splitmix64( )
205+
st(i) = tmp
206+
end do
207+
get = int(tmp, kind = ${k1}$)
208+
seed_initialized = .true.
209+
end subroutine random_distribution_seed_${t1[0]}$${k1}$
210+
211+
#:endfor
212+
end module stdlib_stats_distribution_PRNG

0 commit comments

Comments
 (0)