Skip to content

Commit da63c88

Browse files
Jim-215-FisherJim-215-Fisher
authored
Add files via upload
1 parent 0c54889 commit da63c88

File tree

3 files changed

+45
-101
lines changed

3 files changed

+45
-101
lines changed

‎src/stdlib_stats_distribution_PRNG.fypp

Lines changed: 12 additions & 72 deletions
Original file line numberDiff line numberDiff line change
@@ -1,16 +1,16 @@
11
#:include "common.fypp"
22
module stdlib_stats_distribution_PRNG
3-
use stdlib_kinds
3+
use stdlib_kinds, only: int8, int16, int32, int64
4+
use stdlib_error
45
implicit none
56
private
6-
integer, parameter :: MAX_INT_BIT_SIZE = 64
7-
integer(int64), save :: st(4), si = 614872703977525537_int64
7+
integer, parameter :: MAX_INT_BIT_SIZE = bit_size(1_int64)
8+
integer(int64), save :: st(4) ! internal states for xoshiro256ss function
9+
integer(int64), save :: si = 614872703977525537_int64 ! default seed value
810
logical, save :: seed_initialized = .false.
911

1012
public :: random_seed
1113
public :: dist_rand
12-
public :: jump
13-
public :: long_jump
1414

1515

1616
interface dist_rand
@@ -43,13 +43,16 @@ module stdlib_stats_distribution_PRNG
4343
function dist_rand_${t1[0]}$${k1}$(n) result(res)
4444
!! Random integer generation for various kinds
4545
!! result = [-2^k, 2^k - 1], k = 7, 15, 31, 63, depending on input kind
46-
!! Result is used as bit model number instead of regular arithmetic number
46+
!! Result will be operated by bitwise operators to generate desired integer
47+
!! and real pseudorandom numbers
4748
!!
4849
${t1}$, intent(in) :: n
4950
${t1}$ :: res
5051
integer :: k
5152

5253
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")
5356
res = shiftr(xoshiro256ss( ), k)
5457
end function dist_rand_${t1[0]}$${k1}$
5558

@@ -58,6 +61,7 @@ module stdlib_stats_distribution_PRNG
5861
function xoshiro256ss( ) result (res)
5962
! Generate random 64-bit integers
6063
! Written in 2018 by David Blackman and Sebastiano Vigna ([email protected])
64+
! http://prng.di.unimi.it/xoshiro256starstar.c
6165
!
6266
! This is xoshiro256** 1.0, one of our all-purpose, rock-solid
6367
! generators. It has excellent (sub-ns) speed, a state (256 bits) that is
@@ -81,7 +85,6 @@ module stdlib_stats_distribution_PRNG
8185
st(1) = ieor(st(1), st(4))
8286
st(3) = ieor(st(3), t)
8387
st(4) = rol64(st(4), 45)
84-
return
8588
end function xoshiro256ss
8689

8790
function rol64(x, k) result(res)
@@ -92,73 +95,9 @@ module stdlib_stats_distribution_PRNG
9295
t1 = shiftr(x, (64 - k))
9396
t2 = shiftl(x, k)
9497
res = ior(t1, t2)
95-
return
9698
end function rol64
9799

98100

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 ([email protected])
104-
!
105-
! Fortran 90 version translated from C by Jim-215-Fisher
106-
integer(int64) :: jp(4) = [1733541517147835066_int64, &
107-
-3051731464161248980_int64, &
108-
-6244198995065845334_int64, &
109-
4155657270789760540_int64]
110-
integer(int64) :: s1 = 0, s2 = 0, s3 = 0, s4 = 0, c = 1_int64
111-
integer :: i, j, k
112-
113-
do i = 1, 4
114-
do j = 1, 64
115-
if(iand(jp(i), shiftl(c, j - 1)) /= 0) then
116-
s1 = ieor(s1, st(1))
117-
s2 = ieor(s2, st(2))
118-
s3 = ieor(s3, st(3))
119-
s4 = ieor(s4, st(4))
120-
end if
121-
k = xoshiro256ss( )
122-
end do
123-
end do
124-
st(1) = s1
125-
st(2) = s2
126-
st(3) = s3
127-
st(4) = s4
128-
end subroutine jump
129-
130-
subroutine long_jump
131-
! This is the long-jump function for the xoshiro256ss generator. It is
132-
! equivalent to 2^192 calls to xoshiro256ss(); it can be used to generate
133-
! 2^64 starting points, from each of which jump() will generate 2^64
134-
! non-overlapping subsequences for parallel distributed computations
135-
! Written in 2018 by David Blackman and Sebastiano Vigna ([email protected])
136-
!
137-
! Fortran 90 version translated from C by Jim-215-Fisher
138-
integer(int64) :: jp(4) = [8566230491382795199_int64, &
139-
-4251311993797857357_int64, &
140-
8606660816089834049_int64, &
141-
4111957640723818037_int64]
142-
integer(int64) :: s1 = 0, s2 = 0, s3 = 0, s4 = 0, c = 1_int64
143-
integer(int32) :: i, j, k
144-
145-
do i = 1, 4
146-
do j = 1, 64
147-
if(iand(jp(i), shiftl(c, j - 1)) /= 0) then
148-
s1 = ieor(s1, st(1))
149-
s2 = ieor(s2, st(2))
150-
s3 = ieor(s3, st(3))
151-
s4 = ieor(s4, st(4))
152-
end if
153-
k = xoshiro256ss()
154-
end do
155-
end do
156-
st(1) = s1
157-
st(2) = s2
158-
st(3) = s3
159-
st(4) = s4
160-
end subroutine long_jump
161-
162101
function splitmix64(s) result(res)
163102
! Written in 2015 by Sebastiano Vigna ([email protected])
164103
! This is a fixed-increment version of Java 8's SplittableRandom
@@ -176,14 +115,15 @@ module stdlib_stats_distribution_PRNG
176115
data int01, int02, int03/-7046029254386353131_int64, &
177116
-4658895280553007687_int64, &
178117
-7723592293110705685_int64/
118+
! Values are converted from C unsigned integer of 0x9e3779b97f4a7c15,
119+
! 0xbf58476d1ce4e5b9, 0x94d049bb133111eb
179120

180121
if(present(s)) si = s
181122
res = si
182123
si = res + int01
183124
res = ieor(res, shiftr(res, 30)) * int02
184125
res = ieor(res, shiftr(res, 27)) * int03
185126
res = ieor(res, shiftr(res, 31))
186-
return
187127
end function splitmix64
188128

189129
#:for k1, t1 in INT_KINDS_TYPES

‎src/stdlib_stats_distribution_normal.fypp

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -163,8 +163,8 @@ Module stdlib_stats_distribution_normal
163163
${t1}$ :: x, y
164164
integer :: hz, iz
165165

166-
if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" &
167-
//" parameter must be non-zero")
166+
if(scale==0._${k1}$) call error_stop("Error(norm_dist_rvs): Normal" &
167+
//" distribution scale parameter must be non-zero")
168168
if( .not. zig_norm_initialized ) call zigset
169169
iz = 0
170170
! original algorithm use 32bit
@@ -235,8 +235,8 @@ Module stdlib_stats_distribution_normal
235235
${t1}$ :: x, y, re
236236
integer :: hz, iz, i
237237

238-
if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" &
239-
//" parameter must be non-zero")
238+
if(scale==0._${k1}$) call error_stop("Error(norm_dist_rvs_array):" &
239+
//" Normal distribution scale parameter must be non-zero")
240240
if( .not. zig_norm_initialized ) call zigset
241241
allocate(res(array_size))
242242
do i = 1, array_size
@@ -311,8 +311,8 @@ Module stdlib_stats_distribution_normal
311311
real :: res
312312
${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$ * acos(-1.0_${k1}$))
313313

314-
if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" &
315-
//" parameter must be non-zero")
314+
if(scale==0._${k1}$) call error_stop("Error(norm_dist_pdf):" &
315+
//" Normal distribution scale parameter must be non-zero")
316316
res = exp(- 0.5_${k1}$ * (x - loc) * (x - loc) / (scale * scale)) / &
317317
(sqrt_2_Pi * scale)
318318
return
@@ -342,8 +342,8 @@ Module stdlib_stats_distribution_normal
342342
real :: res
343343
${t1}$, parameter :: sqrt_2 = sqrt(2.0_${k1}$)
344344

345-
if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" &
346-
//" parameter must be non-zero")
345+
if(scale==0._${k1}$) call error_stop("Error(norm_dist_cdf):" &
346+
//" Normal distribution scale parameter must be non-zero")
347347
res = (1.0_${k1}$ + erf((x - loc) / (scale * sqrt_2))) / 2.0_${k1}$
348348
return
349349
end function norm_dist_cdf_${t1[0]}$${k1}$

‎src/stdlib_stats_distribution_uniform.fypp

Lines changed: 25 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -42,7 +42,8 @@ Module stdlib_stats_distribution_uniform
4242
interface uniform_distribution_pdf
4343
!! Version experiment
4444
!!
45-
!! Get uniform distribution probability density (pdf) for integer, real and complex variables
45+
!! Get uniform distribution probability density (pdf) for integer, real and
46+
!! complex variables
4647
!! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html#
4748
!! description))
4849

@@ -54,7 +55,8 @@ Module stdlib_stats_distribution_uniform
5455
interface uniform_distribution_cdf
5556
!! Version experimental
5657
!!
57-
!! Get uniform distribution cumulative distribution function (cdf) for integer, real and complex variables
58+
!! Get uniform distribution cumulative distribution function (cdf) for
59+
!! integer, real and complex variables
5860
!! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html#
5961
!! description))
6062
!!
@@ -66,7 +68,8 @@ Module stdlib_stats_distribution_uniform
6668
interface shuffle
6769
!! Version experimental
6870
!!
69-
!! Fisher-Yates shuffle algorithm for a rank one array of integer, real and complex variables
71+
!! Fisher-Yates shuffle algorithm for a rank one array of integer, real and
72+
!! complex variables
7073
!! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html#
7174
!! description))
7275
!!
@@ -85,13 +88,14 @@ Module stdlib_stats_distribution_uniform
8588
! https://www.pcg-random.org/posts/bounded-rands.html
8689
!
8790
! Fortran 90 translated from c by Jim-215-fisher
91+
!
8892
${t1}$, intent(in) :: scale
8993
${t1}$ :: res, u, mask, n
9094
integer :: zeros, bits_left, bits
9195

9296
n = scale
93-
if(n <= 0_${k1}$) call error_stop("Error: Uniform distribution scale" &
94-
//" parameter must be positive")
97+
if(n <= 0_${k1}$) call error_stop("Error(unif_dist_rvs_1): Uniform" &
98+
//" distribution scale parameter must be positive")
9599
zeros = leadz(n)
96100
bits = bit_size(n) - zeros
97101
mask = shiftr(not(0_${k1}$), zeros)
@@ -121,8 +125,8 @@ Module stdlib_stats_distribution_uniform
121125
${t1}$, intent(in) :: loc, scale
122126
${t1}$ :: res
123127

124-
if(scale == 0_${k1}$) call error_stop("Error: Uniform distribution" &
125-
//" scale parameter must be non-zero")
128+
if(scale <= 0_${k1}$) call error_stop("Error(unif_dist_rvs): Uniform" &
129+
//" distribution scale parameter must be positive")
126130
res = loc + unif_dist_rvs_1_${t1[0]}$${k1}$(scale)
127131
return
128132
end function unif_dist_rvs_${t1[0]}$${k1}$
@@ -153,8 +157,8 @@ Module stdlib_stats_distribution_uniform
153157
${t1}$, intent(in) :: scale
154158
${t1}$ :: res
155159

156-
if(scale == 0._${k1}$) call error_stop("Error: Uniform distribution" &
157-
//" scale parameter must be non-zero")
160+
if(scale == 0._${k1}$) call error_stop("Error(unif_dist_rvs_1): " &
161+
//"Uniform distribution scale parameter must be non-zero")
158162
res = scale * unif_dist_rvs_0_${t1[0]}$${k1}$( )
159163
return
160164
end function unif_dist_rvs_1_${t1[0]}$${k1}$
@@ -169,8 +173,8 @@ Module stdlib_stats_distribution_uniform
169173
${t1}$, intent(in) :: loc, scale
170174
${t1}$ :: res
171175

172-
if(scale == 0._${k1}$) call error_stop("Error: Uniform distribution" &
173-
//" scale parameter must be non-zero")
176+
if(scale == 0._${k1}$) call error_stop("Error(unif_dist_rvs): " &
177+
//"Uniform distribution scale parameter must be non-zero")
174178
res = loc + scale * unif_dist_rvs_0_${t1[0]}$${k1}$( )
175179
return
176180
end function unif_dist_rvs_${t1[0]}$${k1}$
@@ -187,8 +191,8 @@ Module stdlib_stats_distribution_uniform
187191
${t1}$ :: res
188192
real(${k1}$) :: r1, r2, tr, ti
189193

190-
if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" &
191-
//" distribution scale parameter must be non-zero")
194+
if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(uni_dist_" &
195+
//"rvs_1): Uniform distribution scale parameter must be non-zero")
192196
r1 = unif_dist_rvs_0_r${k1}$( )
193197
if(real(scale) == 0.0_${k1}$) then
194198
ti = aimag(scale) * r1
@@ -219,8 +223,8 @@ Module stdlib_stats_distribution_uniform
219223
${t1}$ :: res
220224
real(${k1}$) :: r1, r2, tr, ti
221225

222-
if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" &
223-
//" distribution scale parameter must be non-zero")
226+
if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(uni_dist_" &
227+
//"rvs): Uniform distribution scale parameter must be non-zero")
224228
r1 = unif_dist_rvs_0_r${k1}$( )
225229
if(real(scale) == 0.0_${k1}$) then
226230
tr = real(loc)
@@ -249,8 +253,8 @@ Module stdlib_stats_distribution_uniform
249253
integer :: i, zeros, bits_left, bits
250254

251255
n = scale
252-
if(n == 0_${k1}$) call error_stop("Error: Uniform distribution" &
253-
//" scale parameter must be non-zero")
256+
if(n == 0_${k1}$) call error_stop("Error(unif_dist_rvs_array): Uniform" &
257+
//" distribution scale parameter must be non-zero")
254258
allocate(res(array_size))
255259
zeros = leadz(n)
256260
bits = bit_size(n) - zeros
@@ -287,8 +291,8 @@ Module stdlib_stats_distribution_uniform
287291
integer :: i
288292

289293

290-
if(scale == 0._${k1}$) call error_stop("Error: Uniform distribution" &
291-
//" scale parameter must be non-zero")
294+
if(scale == 0._${k1}$) call error_stop("Error(unif_dist_rvs_array):" &
295+
//" Uniform distribution scale parameter must be non-zero")
292296
allocate(res(array_size))
293297
do i = 1, array_size
294298
tmp = shiftr(dist_rand(INT_ONE), 11)
@@ -311,8 +315,8 @@ Module stdlib_stats_distribution_uniform
311315
integer :: i
312316

313317

314-
if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" &
315-
//" distribution scale parameter must be non-zero")
318+
if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(unif_dist_"&
319+
//"rvs_array): Uniform distribution scale parameter must be non-zero")
316320
allocate(res(array_size))
317321
do i = 1, array_size
318322
tmp = shiftr(dist_rand(INT_ONE), 11)

0 commit comments

Comments
 (0)