Skip to content

Commit fe9aea4

Browse files
Jim-215-FisherJim-215-Fisher
authored
Add files via upload
1 parent b53a270 commit fe9aea4

File tree

3 files changed

+48
-103
lines changed

3 files changed

+48
-103
lines changed

‎src/stdlib_stats_distribution_PRNG.fypp

Lines changed: 7 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -1,16 +1,16 @@
11
#:include "common.fypp"
22
module stdlib_stats_distribution_PRNG
33
use stdlib_kinds, only: int8, int16, int32, int64
4+
use stdlib_error
45
implicit none
56
private
67
integer, parameter :: MAX_INT_BIT_SIZE = bit_size(1_int64)
7-
integer(int64), save :: st(4), si = 614872703977525537_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
@@ -51,6 +51,8 @@ module stdlib_stats_distribution_PRNG
5151
integer :: k
5252

5353
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")
5456
res = shiftr(xoshiro256ss( ), k)
5557
end function dist_rand_${t1[0]}$${k1}$
5658

@@ -96,71 +98,6 @@ module stdlib_stats_distribution_PRNG
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-
! 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 ([email protected])
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-
164101
function splitmix64(s) result(res)
165102
! Written in 2015 by Sebastiano Vigna ([email protected])
166103
! This is a fixed-increment version of Java 8's SplittableRandom
@@ -178,6 +115,8 @@ module stdlib_stats_distribution_PRNG
178115
data int01, int02, int03/-7046029254386353131_int64, &
179116
-4658895280553007687_int64, &
180117
-7723592293110705685_int64/
118+
! Values are converted from C unsigned integer of 0x9e3779b97f4a7c15,
119+
! 0xbf58476d1ce4e5b9, 0x94d049bb133111eb
181120

182121
if(present(s)) si = s
183122
res = si

‎src/stdlib_stats_distribution_exponential.fypp

Lines changed: 16 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
#:include "common.fypp"
22
#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
3-
Module stdlib_stats_distribution_expon
3+
Module stdlib_stats_distribution_exponential
44
use stdlib_kinds
55
use stdlib_error, only : error_stop
66
use stdlib_stats_distribution_PRNG, only : dist_rand
@@ -157,8 +157,9 @@ Module stdlib_stats_distribution_expon
157157
${t1}$ :: r = 7.69711747013104972_${k1}$
158158
integer(int32) :: jz, iz
159159

160-
if(lamda <= 0.0_${k1}$) call error_stop("Error: Exponential" &
161-
//" distribution lamda parameter must be greaeter than zero")
160+
if(lamda <= 0.0_${k1}$) call error_stop("Error(exp_dist_rvs):" &
161+
//" Exponential distribution lamda parameter must be greater than zero")
162+
162163
if( .not. zig_exp_initialized ) call zigset
163164

164165
! Original algorithm use 32bit
@@ -218,8 +219,9 @@ Module stdlib_stats_distribution_expon
218219
${t1}$ :: r = 7.69711747013104972_${k1}$
219220
integer :: jz, iz, i
220221

221-
if(lamda <= 0.0_${k1}$) call error_stop("Error: Exponential" &
222-
//" distribution lamda parameter must be greaeter than zero")
222+
if(lamda <= 0.0_${k1}$) call error_stop("Error(exp_dist_rvs_array):" &
223+
//" Exponential distribution lamda parameter must be greater than zero")
224+
223225
if( .not. zig_exp_initialized ) call zigset
224226
allocate(res(array_size))
225227
do i =1, array_size
@@ -284,10 +286,10 @@ Module stdlib_stats_distribution_expon
284286
${t1}$, intent(in) :: x, lamda
285287
real :: res
286288

287-
if(lamda <= 0.0_${k1}$) call error_stop("Error: Exponential" &
288-
//" distribution lamda parameter must be greaeter than zero")
289-
if(x < 0.0_${k1}$) call error_stop("Error: Exponential distribution" &
290-
//" variate x must be non-negative")
289+
if(lamda <= 0.0_${k1}$) call error_stop("Error(exp_dist_pdf):" &
290+
//" Exponential distribution lamda parameter must be greater than zero")
291+
if(x < 0.0_${k1}$) call error_stop("Error(exp_dist_pdf): Exponential" &
292+
//" distribution variate x must be non-negative")
291293
res = exp(- x * lamda) * lamda
292294
return
293295
end function exp_dist_pdf_${t1[0]}$${k1}$
@@ -313,10 +315,10 @@ Module stdlib_stats_distribution_expon
313315
${t1}$, intent(in) :: x, lamda
314316
real :: res
315317

316-
if(lamda <= 0.0_${k1}$) call error_stop("Error: Exponential" &
317-
//" distribution lamda parameter must be greaeter than zero")
318-
if(x < 0.0_${k1}$) call error_stop("Error: Exponential distribution" &
319-
//" variate x must be non-negative")
318+
if(lamda <= 0.0_${k1}$) call error_stop("Error(exp_dist_cdf):" &
319+
//" Exponential distribution lamda parameter must be greater than zero")
320+
if(x < 0.0_${k1}$) call error_stop("Error(exp_dist_cdf): Exponential" &
321+
//" distribution variate x must be non-negative")
320322
res = (1.0 - exp(- x * lamda))
321323
return
322324
end function exp_dist_cdf_${t1[0]}$${k1}$
@@ -334,4 +336,4 @@ Module stdlib_stats_distribution_expon
334336
end function exp_dist_cdf_${t1[0]}$${k1}$
335337

336338
#:endfor
337-
end module stdlib_stats_distribution_expon
339+
end module stdlib_stats_distribution_exponential

‎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)