Skip to content

Commit 2f45d19

Browse files
Jim-215-FisherJim-215-Fisher
authored
chg. complex number with kinds
1 parent 1a5ba93 commit 2f45d19

File tree

3 files changed

+38
-31
lines changed

3 files changed

+38
-31
lines changed

‎src/stdlib_stats_distribution_PRNG.fypp

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,9 @@
11
#:include "common.fypp"
22
module stdlib_stats_distribution_PRNG
3-
use stdlib_kinds
3+
use stdlib_kinds, only: int8, int16, int32, int64
44
implicit none
55
private
6-
integer, parameter :: MAX_INT_BIT_SIZE = 64
6+
integer, parameter :: MAX_INT_BIT_SIZE = bit_size(1_int64)
77
integer(int64), save :: st(4), si = 614872703977525537_int64
88
logical, save :: seed_initialized = .false.
99

@@ -43,7 +43,8 @@ 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
@@ -58,6 +59,7 @@ module stdlib_stats_distribution_PRNG
5859
function xoshiro256ss( ) result (res)
5960
! Generate random 64-bit integers
6061
! Written in 2018 by David Blackman and Sebastiano Vigna ([email protected])
62+
! http://prng.di.unimi.it/xoshiro256starstar.c
6163
!
6264
! This is xoshiro256** 1.0, one of our all-purpose, rock-solid
6365
! generators. It has excellent (sub-ns) speed, a state (256 bits) that is
@@ -81,7 +83,6 @@ module stdlib_stats_distribution_PRNG
8183
st(1) = ieor(st(1), st(4))
8284
st(3) = ieor(st(3), t)
8385
st(4) = rol64(st(4), 45)
84-
return
8586
end function xoshiro256ss
8687

8788
function rol64(x, k) result(res)
@@ -92,7 +93,6 @@ module stdlib_stats_distribution_PRNG
9293
t1 = shiftr(x, (64 - k))
9394
t2 = shiftl(x, k)
9495
res = ior(t1, t2)
95-
return
9696
end function rol64
9797

9898

@@ -101,6 +101,7 @@ module stdlib_stats_distribution_PRNG
101101
! to 2^128 calls to xoshiro256ss(); it can be used to generate 2^128
102102
! non-overlapping subsequences for parallel computations.
103103
! Written in 2018 by David Blackman and Sebastiano Vigna ([email protected])
104+
! http://prng.di.unimi.it/xoshiro256starstar.c
104105
!
105106
! Fortran 90 version translated from C by Jim-215-Fisher
106107
integer(int64) :: jp(4) = [1733541517147835066_int64, &
@@ -133,6 +134,7 @@ module stdlib_stats_distribution_PRNG
133134
! 2^64 starting points, from each of which jump() will generate 2^64
134135
! non-overlapping subsequences for parallel distributed computations
135136
! Written in 2018 by David Blackman and Sebastiano Vigna ([email protected])
137+
! http://prng.di.unimi.it/xoshiro256starstar.c
136138
!
137139
! Fortran 90 version translated from C by Jim-215-Fisher
138140
integer(int64) :: jp(4) = [8566230491382795199_int64, &
@@ -183,7 +185,6 @@ module stdlib_stats_distribution_PRNG
183185
res = ieor(res, shiftr(res, 30)) * int02
184186
res = ieor(res, shiftr(res, 27)) * int03
185187
res = ieor(res, shiftr(res, 31))
186-
return
187188
end function splitmix64
188189

189190
#:for k1, t1 in INT_KINDS_TYPES

‎src/stdlib_stats_distribution_exponential.fypp

Lines changed: 23 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@ Module stdlib_stats_distribution_expon
2121
!! Version experimental
2222
!!
2323
!! Exponential Distribution Random Variates
24-
!!([Specification](../page/specs/stdlib_stats_distribution_expon.html#
24+
!!([Specification](../page/specs/stdlib_stats_distribution_exponential.html#
2525
!! description))
2626
!!
2727
module procedure exp_dist_rvs_0_rsp !0 dummy variable
@@ -39,7 +39,7 @@ Module stdlib_stats_distribution_expon
3939
!! Version experimental
4040
!!
4141
!! Exponential Distribution Probability Density Function
42-
!!([Specification](../page/specs/stdlib_stats_distribution_expon.html#
42+
!!([Specification](../page/specs/stdlib_stats_distribution_exponential.html#
4343
!! description))
4444
!!
4545
#:for k1, t1 in RC_KINDS_TYPES
@@ -51,7 +51,7 @@ Module stdlib_stats_distribution_expon
5151
!! Version experimental
5252
!!
5353
!! Exponential Distribution Cumulative Distribution Function
54-
!! ([Specification](../page/specs/stdlib_stats_distribution_expon.html#
54+
!! ([Specification](../page/specs/stdlib_stats_distribution_exponential.html#
5555
!! description))
5656
!!
5757
#:for k1, t1 in RC_KINDS_TYPES
@@ -77,6 +77,8 @@ Module stdlib_stats_distribution_expon
7777
! unsigned integers in Fortran.
7878
!
7979
! Latest version - 1 January 2001
80+
!
81+
! Fotran 90 program translated from C by Jim-215-Fisher
8082
!
8183
real(dp), parameter :: M2 = 2147483648.0_dp
8284
real(dp) :: de = 7.697117470131487_dp, te, &
@@ -115,25 +117,25 @@ Module stdlib_stats_distribution_expon
115117

116118
! Original algorithm use 32bit
117119
iz = 0
118-
jz = dist_rand(iz)
120+
jz = dist_rand(1_int32)
119121

120122
iz = iand( jz, 255 )
121123
if( abs( jz ) < ke(iz) ) then
122124
res = abs(jz) * we(iz)
123125
else
124126
L1: do
125127
if( iz == 0 ) then
126-
res = r - log( uni( ) )
128+
res = r - log( uni(1.0_${k1}$) )
127129
exit L1
128130
end if
129131
x = abs( jz ) * we(iz)
130-
if( fe(iz) + uni( ) * (fe(iz-1) - fe(iz)) < exp( -x ) ) then
132+
if(fe(iz) + uni(1.0_${k1}$) * (fe(iz-1) - fe(iz)) < exp(-x)) then
131133
res = x
132134
exit L1
133135
end if
134136

135137
!original algorithm use 32bit
136-
jz = dist_rand(iz)
138+
jz = dist_rand(1_int32)
137139
iz = iand( jz, 255 )
138140
if( abs( jz ) < ke(iz) ) then
139141
res = abs( jz ) * we(iz)
@@ -161,25 +163,25 @@ Module stdlib_stats_distribution_expon
161163

162164
! Original algorithm use 32bit
163165
iz = 0
164-
jz = dist_rand(iz)
166+
jz = dist_rand(1_int32)
165167

166168
iz = iand( jz, 255 )
167169
if( abs( jz ) < ke(iz) ) then
168170
res = abs(jz) * we(iz)
169171
else
170172
L1: do
171173
if( iz == 0 ) then
172-
res = r - log( uni( ) )
174+
res = r - log( uni(1.0_${k1}$) )
173175
exit L1
174176
end if
175177
x = abs( jz ) * we(iz)
176-
if( fe(iz) + uni( ) * (fe(iz-1) - fe(iz)) < exp( -x ) ) then
178+
if(fe(iz) + uni(1.0_${k1}$) * (fe(iz-1) - fe(iz)) < exp(-x)) then
177179
res = x
178180
exit L1
179181
end if
180182

181183
!original algorithm use 32bit
182-
jz = dist_rand(iz)
184+
jz = dist_rand(1_int32)
183185
iz = iand( jz, 255 )
184186
if( abs( jz ) < ke(iz) ) then
185187
res = abs( jz ) * we(iz)
@@ -201,7 +203,7 @@ Module stdlib_stats_distribution_expon
201203

202204
tr = exp_dist_rvs_r${k1}$(real(lamda))
203205
ti = exp_dist_rvs_r${k1}$(aimag(lamda))
204-
res = cmplx(tr, ti)
206+
res = cmplx(tr, ti, kind=${k1}$)
205207
return
206208
end function exp_dist_rvs_${t1[0]}$${k1}$
207209

@@ -223,25 +225,25 @@ Module stdlib_stats_distribution_expon
223225
do i =1, array_size
224226
! Original algorithm use 32bit
225227
iz = 0
226-
jz = dist_rand(iz)
228+
jz = dist_rand(1_int32)
227229

228230
iz = iand( jz, 255 )
229231
if( abs( jz ) < ke(iz) ) then
230232
re = abs(jz) * we(iz)
231233
else
232234
L1: do
233235
if( iz == 0 ) then
234-
re = r - log( uni( ) )
236+
re = r - log( uni(1.0_${k1}$) )
235237
exit L1
236238
end if
237239
x = abs( jz ) * we(iz)
238-
if( fe(iz) + uni( ) * (fe(iz-1) - fe(iz)) < exp( -x ) ) then
240+
if(fe(iz) + uni(1.0_${k1}$)*(fe(iz-1)-fe(iz)) < exp(-x)) then
239241
re = x
240242
exit L1
241243
end if
242244

243245
!original algorithm use 32bit
244-
jz = dist_rand(iz)
246+
jz = dist_rand(1_int32)
245247
iz = iand( jz, 255 )
246248
if( abs( jz ) < ke(iz) ) then
247249
re = abs( jz ) * we(iz)
@@ -268,7 +270,7 @@ Module stdlib_stats_distribution_expon
268270
do i = 1, array_size
269271
tr = exp_dist_rvs_r${k1}$(real(lamda))
270272
ti = exp_dist_rvs_r${k1}$(aimag(lamda))
271-
res(i) = cmplx(tr, ti)
273+
res(i) = cmplx(tr, ti, kind=${k1}$)
272274
end do
273275
return
274276
end function exp_dist_rvs_array_${t1[0]}$${k1}$
@@ -284,6 +286,8 @@ Module stdlib_stats_distribution_expon
284286

285287
if(lamda <= 0.0_${k1}$) call error_stop("Error: Exponential" &
286288
//" 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")
287291
res = exp(- x * lamda) * lamda
288292
return
289293
end function exp_dist_pdf_${t1[0]}$${k1}$
@@ -311,6 +315,8 @@ Module stdlib_stats_distribution_expon
311315

312316
if(lamda <= 0.0_${k1}$) call error_stop("Error: Exponential" &
313317
//" 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")
314320
res = (1.0 - exp(- x * lamda))
315321
return
316322
end function exp_dist_cdf_${t1[0]}$${k1}$

‎src/stdlib_stats_distribution_uniform.fypp

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -134,7 +134,7 @@ Module stdlib_stats_distribution_uniform
134134
! Uniformly distributed float in [0,1]
135135
! Based on the paper by Frederic Goualard, "Generating Random Floating-
136136
! Point Numbers By Dividing Integers: a Case Study", Proceedings of
137-
! ICCS 2020, June 20202, Amsterdam, Netherlands
137+
! ICCS 2020, June 2020, Amsterdam, Netherlands
138138
!
139139
${t1}$ :: res
140140
integer(int64) :: tmp
@@ -201,7 +201,7 @@ Module stdlib_stats_distribution_uniform
201201
tr = real(scale) * r1
202202
ti = aimag(scale) * r2
203203
endif
204-
res = cmplx(tr, ti)
204+
res = cmplx(tr, ti, kind=${k1}$)
205205
return
206206
end function unif_dist_rvs_1_${t1[0]}$${k1}$
207207

@@ -233,7 +233,7 @@ Module stdlib_stats_distribution_uniform
233233
tr = real(loc) + real(scale) * r1
234234
ti = aimag(loc) + aimag(scale) * r2
235235
endif
236-
res = cmplx(tr, ti)
236+
res = cmplx(tr, ti, kind=${k1}$)
237237
return
238238
end function unif_dist_rvs_${t1[0]}$${k1}$
239239

@@ -329,7 +329,7 @@ Module stdlib_stats_distribution_uniform
329329
tr = real(loc) + real(scale) * r1
330330
ti = aimag(loc) + aimag(scale) * r2
331331
endif
332-
res(i) = cmplx(tr, ti)
332+
res(i) = cmplx(tr, ti, kind=${k1}$)
333333
enddo
334334
return
335335
end function unif_dist_rvs_array_${t1[0]}$${k1}$
@@ -343,10 +343,10 @@ Module stdlib_stats_distribution_uniform
343343

344344
if(scale == 0) then
345345
res = 0.0
346-
elseif(x < loc .or. x >loc + scale) then
346+
elseif(x < loc .or. x > (loc + scale)) then
347347
res = 0.0
348348
else
349-
res = 1. / (scale + 1)
349+
res = 1. / (scale + 1_${k1}$)
350350
end if
351351
return
352352
end function unif_dist_pdf_${t1[0]}$${k1}$
@@ -400,7 +400,7 @@ Module stdlib_stats_distribution_uniform
400400
elseif(x < loc) then
401401
res = 0.0
402402
elseif(x >= loc .and. x <= (loc + scale)) then
403-
res = real((x - loc + 1)) / real((scale + 1))
403+
res = real((x - loc + 1_${k1}$)) / real((scale + 1_${k1}$))
404404
else
405405
res = 1.0
406406
end if
@@ -479,4 +479,4 @@ Module stdlib_stats_distribution_uniform
479479
end function shuffle_${t1[0]}$${k1}$
480480

481481
#:endfor
482-
end module stdlib_stats_distribution_uniform
482+
end module stdlib_stats_distribution_uniform

0 commit comments

Comments
 (0)