Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This include file contains the implementation of procedures in [pm_distUnif](@ref pm_distUnif).
19 : !>
20 : !> \author
21 : !> \FatemehBagheri, Wednesday 12:20 PM, September 22, 2021, Dallas, TX
22 :
23 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
24 :
25 : ! Define array elements.
26 : #if setUnifCDF_ENABLED && D0_ENABLED
27 : #define GET_ELEMENT(ARRAY) ARRAY
28 : #elif setUnifCDF_ENABLED && D1_ENABLED
29 : #define GET_ELEMENT(ARRAY) ARRAY(i)
30 : #elif setUnifCDF_ENABLED
31 : #error "Unrecognized interface."
32 : #endif
33 : ! Define the optional RNG argument for the RNG procedures.
34 : #if getUnifRand_ENABLED || setUnifRand_ENABLED
35 : #if RNGD_ENABLED || RNGF_ENABLED
36 : #define RNG
37 : #elif SM64_ENABLED || X256SSG_ENABLED || X256SSW_ENABLED
38 : #define RNG rng,
39 : #else
40 : #error "Unrecognized interface."
41 : #endif
42 : #endif
43 : !%%%%%%%%%%%%%%%%%
44 : #if getUnifCDF_ENABLED
45 : !%%%%%%%%%%%%%%%%%
46 :
47 : #if DD_ENABLED
48 440 : call setUnifCDF(cdf, x)
49 : #elif LU_ENABLED
50 2017 : call setUnifCDF(cdf, x, lower, upper)
51 : #else
52 : #error "Unrecognized interface."
53 : #endif
54 :
55 : !%%%%%%%%%%%%%%%%%
56 : #elif setUnifCDF_ENABLED
57 : !%%%%%%%%%%%%%%%%%
58 :
59 : ! Set default bounds.
60 : #if IK_ENABLED && DD_ENABLED
61 : integer(IKC), parameter :: lower = 0_IKC, upper = 1_IKC
62 : #elif CK_ENABLED && DD_ENABLED
63 : complex(CKC), parameter :: lower = 0._CKC, upper = 1._CKC
64 : #elif RK_ENABLED && DD_ENABLED
65 : real(RKC), parameter :: lower = 0._RKC, upper = 1._RKC
66 : #elif !LU_ENABLED
67 : #error "Unrecognized interface."
68 : #endif
69 : ! Define the normalization constant for vector computations.
70 : #if D1_ENABLED
71 : integer(IK) :: i
72 : #if IK_ENABLED && LU_ENABLED
73 : real(RKC) :: inverseUpperMinusLowerPlusOne
74 1 : inverseUpperMinusLowerPlusOne = 1._RKC / real(upper - lower + 1_IKC, RKC)
75 : #elif CK_ENABLED && LU_ENABLED
76 : complex(CKC) :: inverseUpperMinusLowerPlusOne
77 1 : inverseUpperMinusLowerPlusOne%re = 1._CKC / (upper%re - lower%re)
78 1 : inverseUpperMinusLowerPlusOne%im = 1._CKC / (upper%im - lower%im)
79 : #elif RK_ENABLED && LU_ENABLED
80 : real(RKC) :: inverseUpperMinusLowerPlusOne
81 1 : inverseUpperMinusLowerPlusOne = 1._RKC / (upper - lower)
82 : #elif !DD_ENABLED
83 : #error "Unrecognized interface."
84 : #endif
85 9 : CHECK_ASSERTION(__LINE__, size(cdf, kind = IK) == size(X, kind = IK), SK_"setUnifCDF(): The condition `size(cdf) == size(X)` must hold. size(cdf), size(X) = "//getStr([size(cdf, kind = IK), size(X, kind = IK)])) ! fpp
86 : #elif !D0_ENABLED
87 : #error "Unrecognized interface."
88 : #endif
89 : #if LU_ENABLED && CK_ENABLED
90 3009 : CHECK_ASSERTION(__LINE__, lower%re <= upper%re .and. lower%im <= upper%im, SK_"setUnifCDF(): The conditions `lower%re <= upper%re .and. lower%im <= upper%im` must hold. lower, upper = "//getStr([lower, upper])) ! fpp
91 : #elif LU_ENABLED
92 3060 : CHECK_ASSERTION(__LINE__, lower <= upper, SK_"setUnifCDF(): The condition `lower <= upper` must hold. lower, upper = "//getStr([lower, upper])) ! fpp
93 : #elif !DD_ENABLED
94 : #error "Unrecognized interface."
95 : #endif
96 : ! Begin the computation.
97 : #if D1_ENABLED
98 : do concurrent(i = 1 : size(X, kind = IK))
99 : #endif
100 : ! integer.
101 : #if IK_ENABLED
102 33 : if (GET_ELEMENT(X) < lower) then
103 6 : GET_ELEMENT(cdf) = 0._RKC
104 26 : elseif (GET_ELEMENT(X) < upper) then
105 : #if DD_ENABLED
106 2 : GET_ELEMENT(cdf) = real(GET_ELEMENT(X) + 1_IKC, RKC) * 0.5_RKC
107 : #elif LU_ENABLED && D0_ENABLED
108 9 : GET_ELEMENT(cdf) = real(GET_ELEMENT(X) + 1_IKC - lower, RKC) / real(upper - lower + 1_IKC, RKC)
109 : #elif LU_ENABLED && D1_ENABLED
110 7 : GET_ELEMENT(cdf) = real(GET_ELEMENT(X) + 1_IKC - lower, RKC) * inverseUpperMinusLowerPlusOne
111 : #else
112 : #error "Unrecognized interface."
113 : #endif
114 : else
115 8 : GET_ELEMENT(cdf) = 1._RKC
116 : end if
117 : ! real.
118 : #elif RK_ENABLED
119 2442 : if (GET_ELEMENT(X) < lower) then
120 500 : GET_ELEMENT(cdf) = 0._RKC
121 1941 : elseif (GET_ELEMENT(X) < upper) then
122 : #if DD_ENABLED
123 439 : GET_ELEMENT(cdf) = GET_ELEMENT(X)
124 : #elif LU_ENABLED && D0_ENABLED
125 502 : GET_ELEMENT(cdf) = (GET_ELEMENT(X) - lower) / (upper - lower)
126 : #elif LU_ENABLED && D1_ENABLED
127 500 : GET_ELEMENT(cdf) = (GET_ELEMENT(X) - lower) * inverseUpperMinusLowerPlusOne
128 : #else
129 : #error "Unrecognized interface."
130 : #endif
131 : else
132 500 : GET_ELEMENT(cdf) = 1._RKC
133 : end if
134 : ! complex.
135 : #elif CK_ENABLED
136 : ! real part.
137 2004 : if (GET_ELEMENT(X)%re < real(lower,CKC)) then
138 500 : GET_ELEMENT(cdf)%re = 0._CKC ! fpp
139 1504 : elseif (GET_ELEMENT(X)%re < real(upper,CKC)) then
140 : #if DD_ENABLED
141 2 : GET_ELEMENT(cdf)%re = GET_ELEMENT(X)%re
142 : #elif LU_ENABLED && D0_ENABLED
143 502 : GET_ELEMENT(cdf)%re = (GET_ELEMENT(X)%re - real(lower,CKC)) / (real(upper,CKC) - real(lower,CKC))
144 : #elif LU_ENABLED && D1_ENABLED
145 500 : GET_ELEMENT(cdf)%re = (GET_ELEMENT(X)%re - real(lower,CKC)) * inverseUpperMinusLowerPlusOne%re
146 : #else
147 : #error "Unrecognized interface."
148 : #endif
149 : else
150 500 : GET_ELEMENT(cdf)%re = 1._CKC ! fpp
151 : end if
152 : ! imaginary part.
153 3005 : if (GET_ELEMENT(X)%im < aimag(lower)) then
154 500 : GET_ELEMENT(cdf)%im = 0._CKC ! fpp
155 1502 : elseif (GET_ELEMENT(X)%im < aimag(upper)) then
156 : #if DD_ENABLED
157 : GET_ELEMENT(cdf)%im = GET_ELEMENT(X)%im
158 : #elif LU_ENABLED && D0_ENABLED
159 502 : GET_ELEMENT(cdf)%im = (GET_ELEMENT(X)%im - aimag(lower)) / (aimag(upper) - aimag(lower))
160 : #elif LU_ENABLED && D1_ENABLED
161 500 : GET_ELEMENT(cdf)%im = (GET_ELEMENT(X)%im - aimag(lower)) * inverseUpperMinusLowerPlusOne%im
162 : #else
163 : #error "Unrecognized interface."
164 : #endif
165 : else
166 502 : GET_ELEMENT(cdf)%im = 1._CKC ! fpp
167 : end if
168 : #else
169 : #error "Unrecognized interface."
170 : #endif
171 : #if D1_ENABLED
172 : end do
173 : #endif
174 :
175 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
176 : #elif constructSplitmix64_ENABLED
177 : !%%%%%%%%%%%%%%%%%%%%%%%%%%
178 :
179 : integer(IK64) :: count
180 387381 : if (present(seed)) then
181 13 : rng%state = seed
182 : else ! By default, the seed is randomly initialized for every new instance of the RNG.
183 387368 : call system_clock(count)
184 : rng%state = 324108011427370141_IK64 ! This must be present, otherwise GNU 10.3 uninitliazation warning bug.
185 387368 : rng%state = ieor(rng%state, count)
186 : end if
187 387381 : if (present(imageID)) then
188 0 : CHECK_ASSERTION(__LINE__, 0_IK < imageID, \
189 : SK_"@constructSplitmix64(): The condition `0 < imageID` must hold. imageID = "//getStr(imageID))
190 0 : rng%state = ieor(rng%state, int(imageID, IK64))
191 : else
192 387381 : rng%state = ieor(rng%state, 1_IK64)
193 : end if
194 :
195 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
196 : #elif setStateNext_ENABLED && SM64_ENABLED
197 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
198 :
199 : ! Equivalent to unsigned hexadecimal integers
200 : ! Z"9e3779b97f4a7c15", Z"bf58476d1ce4e5b9", Z"94d049bb133111eb"
201 : integer(IK64), parameter :: TRIPLE(3) = [ -7046029254386353131_IK64 &
202 : , -4658895280553007687_IK64 &
203 : , -7723592293110705685_IK64 ]
204 1569568 : if (rng%state < 0_IK64) then
205 784777 : if (rng%state < -huge(0_IK64) - 1_IK64 - TRIPLE(1)) then
206 782407 : rng%state = rng%state + huge(0_IK64) + 1_IK64
207 782407 : rng%state = rng%state + TRIPLE(1)
208 782407 : rng%state = rng%state + huge(0_IK64) + 1_IK64
209 : else
210 2370 : rng%state = rng%state + TRIPLE(1)
211 : end if
212 : else
213 784791 : rng%state = rng%state + TRIPLE(1)
214 : end if
215 : !rng%state = rng%state + TRIPLE(1)
216 1569568 : rng%stream = rng%state
217 1569568 : rng%stream = ieor(rng%stream, shiftr(rng%stream, 30)) * TRIPLE(2)
218 1569568 : rng%stream = ieor(rng%stream, shiftr(rng%stream, 27)) * TRIPLE(3)
219 1569568 : rng%stream = ieor(rng%stream, shiftr(rng%stream, 31))
220 :
221 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
222 : #elif constructXoshiro256ssg_ENABLED || constructXoshiro256ssw_ENABLED
223 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
224 :
225 : use pm_kind, only: IKC => IK64, RKC => RK64
226 : integer(IKC) :: ijump
227 : type(splitmix64_type) :: rngsm
228 387367 : rngsm = splitmix64_type(seed = seed)
229 387367 : call setUnifRand(rngsm, rng%state)
230 387367 : call setStateJump(rng)
231 387367 : if (present(imageID)) then
232 387179 : CHECK_ASSERTION(__LINE__, 0_IK < imageID, SK_"@constructXoshiro256ssw(): The condition `0 < imageID` must hold. imageID = "//getStr(imageID))
233 387179 : do ijump = 2_IKC, imageID
234 387179 : call setStateJump(rng)
235 : end do
236 : end if
237 :
238 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
239 : #elif setStateNext_ENABLED && (X256SSG_ENABLED || X256SSW_ENABLED)
240 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
241 :
242 : integer(IK64) :: dum
243 115810840 : rng%stream = ishftc(rng%state(2) * 5_IK64, 7) * 9_IK64
244 115810840 : dum = shiftl(rng%state(2), 17)
245 115810840 : rng%state(3) = ieor(rng%state(3), rng%state(1))
246 115810840 : rng%state(4) = ieor(rng%state(4), rng%state(2))
247 115810840 : rng%state(2) = ieor(rng%state(2), rng%state(3))
248 115810840 : rng%state(1) = ieor(rng%state(1), rng%state(4))
249 115810840 : rng%state(3) = ieor(rng%state(3), dum)
250 115810840 : rng%state(4) = ishftc(rng%state(4), 45)
251 : #if X256SSG_ENABLED
252 16554 : rng%pos = 0_IK
253 : #endif
254 :
255 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
256 : #elif setStateJump_ENABLED && (X256SSG_ENABLED || X256SSW_ENABLED)
257 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
258 :
259 : integer(IK64) :: state(4)
260 : integer(IK) :: istate, ibit
261 : #if DJ_ENABLED
262 : #define JUMP xoshiro256ssJump128
263 : #elif AJ_ENABLED
264 0 : CHECK_ASSERTION(__LINE__, size(jump, 1, IK) == size(rng%state, 1, IK), \
265 : SK_": The condition `size(jump, 1) == size(rng%state, 1)` must hold. size(jump), size(rng%state) = "\
266 : //getStr([size(jump, 1, IK), size(rng%state, 1, IK)]))
267 : #else
268 : #error "Unreocgnized interface."
269 : #endif
270 387367 : state = 0_IK64
271 : ! Jump 2^128 (or 2^64) steps ahead from the specified `jump`.
272 1936835 : do istate = 1_IK, size(rng%state, 1, IK)
273 101102787 : do ibit = 0_IK, int(bit_size(rng%stream), IK) - 1_IK ! 63_IK
274 335459822 : if (btest(JUMP(istate), ibit)) state = ieor(state, rng%state) ! fpp
275 100715420 : call setStateNext(rng)
276 : end do
277 : end do
278 1936835 : rng%state = state
279 : #undef JUMP
280 :
281 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
282 : #elif getUnifRand_ENABLED && D0_ENABLED && LK_ENABLED && (RNGD_ENABLED || RNGF_ENABLED) && DD_ENABLED
283 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
284 :
285 : real :: dummy
286 362308 : call random_number(dummy)
287 362308 : rand = logical(dummy < .5, LKC)
288 :
289 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
290 : #elif getUnifRand_ENABLED && D0_ENABLED && LK_ENABLED && (SM64_ENABLED || X256SSW_ENABLED) && DD_ENABLED
291 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
292 :
293 0 : rand = logical(int(0, kind(rng%stream)) < rng%stream, LKC)
294 0 : call setStateNext(rng)
295 :
296 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
297 : #elif getUnifRand_ENABLED && D0_ENABLED && LK_ENABLED && X256SSG_ENABLED && DD_ENABLED
298 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
299 :
300 0 : rand = logical(btest(rng%stream, rng%pos), LKC)
301 0 : rng%pos = rng%pos + 1_IK
302 0 : if (rng%pos == xoshiro256ssStreamBitSize) call setStateNext(rng)
303 :
304 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
305 : #elif getUnifRand_ENABLED && (RNGD_ENABLED || RNGF_ENABLED || SM64_ENABLED || X256SSG_ENABLED || X256SSW_ENABLED) && LU_ENABLED
306 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
307 :
308 13347423 : call setUnifRand(RNG rand, lb, ub)
309 :
310 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
311 : #elif setUnifRand_ENABLED && D0_ENABLED && SK_ENABLED && DD_ENABLED
312 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
313 :
314 : integer(IK) :: i, randint
315 5098929 : do i = 1_IK, len(rand, IK)
316 3799559 : call setUnifRand(RNG randint, 1_IK, 127_IK)
317 5098929 : rand(i:i) = char(randint)
318 : end do
319 :
320 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
321 : #elif setUnifRand_ENABLED && D0_ENABLED && SK_ENABLED && LU_ENABLED
322 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
323 :
324 : integer(IK) :: lbLen, ubLen
325 : integer(IK) :: i, randint
326 : integer(IK) :: lbi, ubi
327 1237326 : lbLen = len(lb, IK)
328 1237326 : ubLen = len(ub, IK)
329 3711978 : CHECK_ASSERTION(__LINE__, lb <= ub, SK_"@setUnifRand(): The condition `lb <= ub` must hold. lb, ub = "//getStr([lb, ub]))
330 3711978 : CHECK_ASSERTION(__LINE__, lbLen == ubLen .or. lbLen == 1_IK .or. ubLen == 1_IK, SK_"@setUnifRand(): The condition `len(lb) == len(ub) .or. len(lb) == 1 .or. len(ub) == 1` must hold. len(lb), len(ub) = "//getStr([lbLen, ubLen]))
331 3711978 : CHECK_ASSERTION(__LINE__, lbLen == len(rand, IK) .or. lbLen == 1_IK, SK_"@setUnifRand(): The condition `len(lb) == len(rand) .or. len(lb) == 1` must hold. len(lb), len(rand) = "//getStr([lbLen, len(rand , IK)]))
332 3711978 : CHECK_ASSERTION(__LINE__, ubLen == len(rand, IK) .or. ubLen == 1_IK, SK_"@setUnifRand(): The condition `len(ub) == len(rand) .or. len(ub) == 1` must hold. len(ub), len(rand) = "//getStr([ubLen, len(rand , IK)]))
333 1237326 : if (1_IK < lbLen .and. 1_IK < ubLen) then
334 2610666 : do i = 1_IK, len(rand, IK)
335 1777015 : call setUnifRand(RNG randint, ichar(lb(i:i), IK), ichar(ub(i:i), IK))
336 2610666 : rand(i:i) = char(randint)
337 : end do
338 403675 : elseif (1_IK == lbLen .and. 1_IK == ubLen) then
339 403333 : lbi = ichar(lb, IK)
340 403333 : ubi = ichar(ub, IK)
341 8357461 : do i = 1_IK, len(rand, IK)
342 7954128 : call setUnifRand(RNG randint, lbi, ubi)
343 8357461 : rand(i:i) = char(randint)
344 : end do
345 342 : elseif (1_IK == lbLen) then
346 0 : lbi = ichar(lb, IK)
347 0 : do i = 1_IK, len(rand, IK)
348 0 : call setUnifRand(RNG randint, lbi, ichar(ub(i:i), IK))
349 0 : rand(i:i) = char(randint)
350 : end do
351 342 : elseif (1_IK == ubLen) then
352 0 : ubi = ichar(ub, IK)
353 0 : do i = 1_IK, len(rand, IK)
354 0 : call setUnifRand(RNG randint, ichar(lb(i:i), IK), ubi)
355 0 : rand(i:i) = char(randint)
356 : end do
357 342 : elseif (len(rand, IK) /= 0_IK) then
358 : error stop "@setUnifRand(): Invalid user-specified input arguments. "& ! LCOV_EXCL_LINE
359 : //"Recompile & rerun with macro CHECK_ENABLED=1 for more information." ! LCOV_EXCL_LINE
360 : end if
361 :
362 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
363 : #elif setUnifRand_ENABLED && D0_ENABLED && IK_ENABLED && (RNGD_ENABLED || RNGF_ENABLED)
364 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
365 :
366 : ! Convert significant bit count to decimal precision.
367 : ! RKT: target real kind with at least `digits(rand)` bits.
368 : ! This is to ensure full coverage of the range of the specific-kind integer.
369 : integer, parameter :: RKT = selected_real_kind(floor(digits(rand) * log10(2._RKB)))
370 : integer, parameter :: RKC = merge(RKB, RKT, -6 < RKT .and. RKT < 0)
371 : real(RKC) :: temp
372 : #if DD_ENABLED
373 : !real(RKD) :: temp
374 : !call random_number(temp)
375 : ! \bug GNU Fortran compiler 10.3
376 : ! The following comment is a bug with `IKC = integer_kinds(5)`.
377 : !rand = floor(.5_RKD + temp, kind = IKC) ! rand = nint(temp, kind = IKC)
378 : integer(IKC), parameter :: lb = -huge(0_IKC), ub = +huge(0_IKC)
379 : #elif LU_ENABLED
380 41704167 : CHECK_ASSERTION(__LINE__, lb <= ub, SK_"@setUnifRand(): The condition `lb <= ub` must hold. lb, ub = "//getStr([lb, ub]))
381 : #else
382 : #error "Unrecognized interface."
383 : #endif
384 13901389 : if (lb + 1_IKC < ub) then
385 19862439 : call random_number(temp)
386 : ! early conversion to `real` avoids possible overflow with `huge` limits.
387 : ! rand = lb + int(temp * real(ub - lb + 1_IKC, kind(temp)), kind = IKC)
388 19862439 : rand = floor((1._RKC - temp) * real(lb, RKC) + temp * real(ub, RKC) + temp, kind = IKC)
389 384448 : elseif (lb == ub) then
390 189306 : rand = lb
391 : else
392 195142 : call random_number(temp)
393 195142 : if (temp < 0.5_RKC) then
394 97717 : rand = lb
395 : else
396 97425 : rand = ub
397 : end if
398 : end if
399 :
400 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
401 : #elif setUnifRand_ENABLED && D0_ENABLED && IK_ENABLED && (SM64_ENABLED || X256SSG_ENABLED || X256SSW_ENABLED)
402 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
403 :
404 : integer , parameter :: IKS = kind(rng%stream)
405 : integer(IK) , parameter :: streamBitSize = bit_size(rng%stream) ! this has to remain generic, supportive of all RNGs.
406 : #if DD_ENABLED
407 : integer(IK) , parameter :: randBitSize = bit_size(rand)
408 : integer(IK) , parameter :: streamBitExcess = streamBitSize - randBitSize
409 : integer(IK) , parameter :: randStreamBitSizeRatio = int(real(randBitSize) / real(streamBitSize), IKC) + 1_IK
410 : integer(IKS) :: buffer(randStreamBitSizeRatio)
411 : integer(IK) :: ibuf
412 : #if X256SSG_ENABLED
413 10061 : ibuf = rng%pos + randBitSize
414 10061 : if (ibuf < streamBitSize) then
415 8796 : rand = int(ibits(rng%stream, rng%pos, randBitSize), IKC)
416 8796 : rng%pos = ibuf
417 10057 : return
418 1265 : elseif (ibuf == streamBitSize) then
419 1261 : rand = int(ibits(rng%stream, rng%pos, randBitSize), IKC)
420 1261 : call setStateNext(rng)
421 1261 : return
422 : end if
423 : ! Here, the `rand` kind is higher than kind(rng%stream).
424 : ! There is currently no neat solution in Fortran for greedy storage of the remaining bits.
425 : ! For now, accept the wasteful approach below for the greedy algorithm too.
426 4 : if (0 < rng%pos) call setStateNext(rng)
427 : #else
428 : if (0_IK < streamBitExcess) then
429 : ! Both shifting and transfer below are possible solutions.
430 : ! The `abs` below, though redundant, bypasses the standard constraint.
431 : !rand = int(shiftr(rng%stream, abs(streamBitExcess)), IKC)
432 11718880 : rand = transfer(source = rng%stream, mold = rand)
433 11718880 : call setStateNext(rng)
434 : return
435 : elseif (0_IK == streamBitExcess) then
436 1549468 : rand = int(rng%stream, IKC)
437 1549468 : call setStateNext(rng)
438 : return
439 : end if
440 : #endif
441 8 : do ibuf = 1_IK, randStreamBitSizeRatio
442 4 : buffer(ibuf) = rng%stream
443 8 : call setStateNext(rng)
444 : end do
445 4 : rand = transfer(source = buffer, mold = rand)
446 : #elif LU_ENABLED
447 : integer(IKC), parameter :: HUGE_IKC = huge(0_IKC)
448 : integer(IKC) :: scale, lbb, ubb, nzeros, nsignif, imask, temp, diff
449 35150460 : CHECK_ASSERTION(__LINE__, lb <= ub, SK_"@setUnifRand(): The condition `lb <= ub` must hold. lb, ub = "//getStr([lb, ub]))
450 11716820 : if (lb /= ub) then
451 11716534 : if (lb + 1_IKC == ub) then ! impossible range overflow.
452 : #if X256SSG_ENABLED
453 5 : if (btest(rng%stream, rng%pos)) then
454 1 : rand = lb
455 : else
456 4 : rand = ub
457 : end if
458 5 : rng%pos = rng%pos + 1_IK
459 5 : if (rng%pos == streamBitSize) call setStateNext(rng)
460 : #else
461 643 : if (0_IKS < rng%stream) then
462 312 : rand = lb
463 : else
464 331 : rand = ub
465 : end if
466 643 : call setStateNext(rng)
467 : #endif
468 648 : return
469 : end if
470 11715886 : if (lb < 0_IKC .and. 0_IKC < ub) then ! possible range overflow.
471 10045 : if (HUGE_IKC + lb < ub) then ! overflowed.
472 15 : lbb = -lb - HUGE_IKC - 1_IKC
473 15 : ubb = +ub - HUGE_IKC - 1_IKC
474 15 : scale = ubb + lbb
475 0 : nzeros = 0_IKC
476 : else
477 10030 : scale = ub - lb
478 10030 : nzeros = int(leadz(scale), IKC)
479 : end if
480 : !scale = ub .uadd. -lb
481 : else ! impossible range overflow (assuming `lb < ub`).
482 11705841 : scale = ub - lb
483 11705841 : nzeros = int(leadz(scale), IKC)
484 : end if
485 11715886 : nsignif = bit_size(scale) - nzeros
486 11715886 : imask = shiftr(not(0_IKC), nzeros)
487 : loopTry: do
488 11718931 : call setUnifRand(rng, temp)
489 11718931 : rand = iand(temp, imask)
490 11718931 : if(rand <= scale) exit loopTry
491 0 : diff = nzeros
492 11715886 : loopReject: do
493 1848807 : if(diff < nsignif) exit loopReject
494 1845762 : temp = shiftr(temp, nsignif)
495 1845762 : rand = iand(temp, imask)
496 1845762 : if(rand <= scale) exit loopTry
497 351913 : diff = diff - nsignif
498 : end do loopReject
499 : end do loopTry
500 11715886 : rand = rand + lb
501 : else
502 286 : rand = lb
503 : end if
504 : #else
505 : #error "Unrecognized interface."
506 : #endif
507 :
508 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
509 : #elif setUnifRand_ENABLED && D0_ENABLED && LK_ENABLED && (RNGD_ENABLED || RNGF_ENABLED)
510 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
511 :
512 : #if DD_ENABLED
513 : real :: dummy
514 6338811 : call random_number(dummy)
515 6338811 : rand = logical(dummy < .5, kind = LKC)
516 : #elif LU_ENABLED
517 : real :: dummy
518 3987920 : if (lb .neqv. ub) then
519 3987920 : call random_number(dummy)
520 3987920 : rand = logical(dummy < .5, kind = LKC)
521 : else
522 0 : rand = lb
523 : end if
524 : #else
525 : #error "Unrecognized interface."
526 : #endif
527 :
528 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
529 : #elif setUnifRand_ENABLED && D0_ENABLED && LK_ENABLED && (SM64_ENABLED || X256SSG_ENABLED || X256SSW_ENABLED)
530 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
531 :
532 : #if DD_ENABLED && X256SSG_ENABLED
533 : integer(IK) , parameter :: streamBitSize = bit_size(rng%stream)
534 10 : rand = logical(btest(rng%stream, rng%pos), LKC)
535 10 : rng%pos = rng%pos + 1_IK
536 10 : if (rng%pos == streamBitSize) call setStateNext(rng)
537 : #elif DD_ENABLED && (SM64_ENABLED || X256SSW_ENABLED)
538 20 : rand = logical(int(0, kind(rng%stream)) < rng%stream, kind = LKC)
539 20 : call setStateNext(rng)
540 : #elif LU_ENABLED
541 30 : if (lb .neqv. ub) then
542 30 : call setUnifRand(rng, rand)
543 : else
544 0 : rand = lb
545 : end if
546 : #else
547 : #error "Unrecognized interface."
548 : #endif
549 :
550 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
551 : #elif setUnifRand_ENABLED && D0_ENABLED && CK_ENABLED
552 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
553 :
554 : ! Curse you Intel.
555 : #if __INTEL_COMPILER
556 : real(CKC) :: temp(2)
557 : #if DD_ENABLED
558 : call setUnifRand(RNG temp)
559 : #elif LU_ENABLED
560 : if (lb /= ub) then
561 : call setUnifRand(RNG temp(1), lb%re, ub%re)
562 : call setUnifRand(RNG temp(2), lb%im, ub%im)
563 : else
564 : rand = lb
565 : end if
566 : #else
567 : #error "Unrecognized interface."
568 : #endif
569 : rand = cmplx(temp(1), temp(2), CKC)
570 : #else
571 : #if DD_ENABLED
572 5082527 : call setUnifRand(RNG rand%re)
573 5082527 : call setUnifRand(RNG rand%im)
574 : #elif LU_ENABLED
575 4374602 : if (lb /= ub) then
576 4373708 : call setUnifRand(RNG rand%re, lb%re, ub%re)
577 4373708 : call setUnifRand(RNG rand%im, lb%im, ub%im)
578 : else
579 894 : rand = lb
580 : end if
581 : #else
582 : #error "Unrecognized interface."
583 : #endif
584 : #endif
585 :
586 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
587 : #elif setUnifRand_ENABLED && D0_ENABLED && RK_ENABLED && (RNGD_ENABLED || RNGF_ENABLED)
588 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
589 :
590 : #if DD_ENABLED
591 16298499 : call random_number(rand)
592 : #elif LU_ENABLED
593 19728941 : if (lb /= ub) then
594 59171328 : CHECK_ASSERTION(__LINE__, lb < ub, SK_"@setUnifRand(): The condition `lb <= ub` must hold. lb, ub = "//getStr([lb, ub]))
595 : ! The following loop ensures the random number falls in the half-open range `[lb, ub)`,
596 : ! even in the extreme case when `ub = lb + spacing(lb)`.
597 : do
598 19723776 : call random_number(rand)
599 : !rand = lb + rand * (ub - lb)
600 : ! The product expansion, although more expensive, avoids possible overflow with `huge` limits.
601 19723776 : rand = (1._RKC - rand) * lb + rand * ub
602 : ! The equality (instead of <) ensures NAN cases are handled gracefully without infinite loop.
603 : ! Hard lesson learned.
604 : !if (rand < ub) return
605 19723776 : if (ub <= rand) cycle
606 0 : exit
607 : end do
608 : else
609 5165 : rand = lb
610 : end if
611 : #else
612 : #error "Unrecognized interface."
613 : #endif
614 :
615 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
616 : #elif setUnifRand_ENABLED && D0_ENABLED && RK_ENABLED && (SM64_ENABLED || X256SSG_ENABLED || X256SSW_ENABLED)
617 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
618 :
619 : #if DD_ENABLED
620 : integer , parameter :: IKS = kind(rng%stream)
621 : integer , parameter :: DIGITS_RKC = digits(rand)
622 : integer , parameter :: DIGITS_RNG = digits(rng%stream)
623 : !#if X256SSG_ENABLED
624 : ! The following original approach is too costly.
625 : ! We will instead use the wasteful approach for the greedy.
626 : !real(RKC) , parameter :: INVPDIGRKC = 1._RKC / 2._RKC**DIGITS_RKC
627 : !integer :: remaining, remainders
628 : !remainders = DIGITS_RNG - DIGITS_RKC - rng%pos
629 : !if (0 < remainders) then
630 : ! rand = real(ibits(rng%stream, rng%pos, DIGITS_RKC), RKC) * INVPDIGRKC
631 : ! rng%pos = rng%pos + DIGITS_RKC
632 : !else
633 : ! remainders = DIGITS_RNG - rng%pos
634 : ! rand = real(ibits(rng%stream, rng%pos, remainders), RKC) * 0.5_RKC**remainders
635 : ! call setStateNext(rng)
636 : ! remaining = DIGITS_RKC - remainders
637 : ! do
638 : ! remainders = remaining - DIGITS_RNG
639 : ! if (0 < remainders) then ! use full stream and cycle
640 : ! rand = (rand + real(abs(rng%stream), RKC)) * 0.5_RKC**DIGITS_RNG
641 : ! remaining = remainders
642 : ! call setStateNext(rng)
643 : ! else ! use part of stream and exit.
644 : ! rng%pos = -remainders
645 : ! rand = (rand + real(ibits(rng%stream, 0_IK, rng%pos), RKC)) * 0.5_RKC**rng%pos
646 : ! return
647 : ! end if
648 : ! end do
649 : !end if
650 : !#else
651 : integer :: i
652 : integer , parameter :: NUMBIT_RNG = bit_size(rng%stream) ! DIGITS_RNG + 1 ! The `1` makes up for the missing sign bit in the counting.
653 : integer , parameter :: REPEAT_RNG = int(real(DIGITS_RKC) / real(DIGITS_RNG))
654 : integer , parameter :: REMAINDERS = DIGITS_RKC - DIGITS_RNG * REPEAT_RNG
655 : ! bit-shifting or integer exponentiation overflows for exponent 63. Use real exponentation.
656 : real(RKC) , parameter :: INV_POWREM = 1._RKC / 2._RKC**REMAINDERS ! real(shiftl(1_IKS, REMAINDERS), RKC) ! 1. / 2_IKS**REMAINDERS
657 : real(RKC) , parameter :: INV_POWDIG = 1._RKC / 2._RKC**DIGITS_RNG ! real(shiftl(1_IKS, DIGITS_RNG), RKC) ! 1. / 2_IKS**DIGITS_RNG
658 : #if X256SSG_ENABLED
659 15022 : if (DIGITS_RNG < REMAINDERS + rng%pos) call setStateNext(rng)
660 15022 : rand = real(shiftr(rng%stream, NUMBIT_RNG - REMAINDERS + rng%pos), RKC) * INV_POWREM
661 : #else
662 4929130 : rand = real(shiftr(rng%stream, NUMBIT_RNG - REMAINDERS), RKC) * INV_POWREM
663 : #endif
664 4944152 : call setStateNext(rng)
665 48 : do i = 1, REPEAT_RNG ! exists only for higher-than-double precision real kinds.
666 24 : rand = (rand + real(shiftr(rng%stream, NUMBIT_RNG - DIGITS_RNG), RKC)) * INV_POWDIG
667 48 : call setStateNext(rng)
668 : end do
669 : !#endif
670 : #elif LU_ENABLED
671 3779256 : if (lb /= ub) then
672 11337768 : CHECK_ASSERTION(__LINE__, lb < ub, SK_"@setUnifRand(): The condition `lb <= ub` must hold. lb, ub = "//getStr([lb, ub]))
673 : ! The following loop ensures the random number falls in the half-open range `[lb, ub)`,
674 : ! even in the extreme case when `ub = lb + spacing(lb)`.
675 : do
676 3779265 : call setUnifRand(rng, rand)
677 : ! The product expansion, although more expensive, avoids possible overflow with `huge` limits.
678 3779265 : rand = (1._RKC - rand) * lb + rand * ub
679 : !rand = lb + rand * (ub - lb)
680 3779265 : if (ub <= rand) cycle
681 9 : exit
682 : end do
683 : else
684 0 : rand = lb
685 : end if
686 : #else
687 : #error "Unrecognized interface."
688 : #endif
689 :
690 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
691 : #elif setUnifRand_ENABLED && !D0_ENABLED && (RNGD_ENABLED || RNGF_ENABLED || SM64_ENABLED || X256SSG_ENABLED || X256SSW_ENABLED)
692 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
693 :
694 : ! Define the bounds.
695 : #if DD_ENABLED
696 : #define BOUNDS
697 : #elif LU_ENABLED
698 : #define BOUNDS, lb, ub
699 : #else
700 : #error "Unrecognized interface."
701 : #endif
702 : #if D1_ENABLED
703 : integer(IK) :: irow, nrow
704 387385 : nrow = size(rand, 1, IK)
705 2004357 : do irow = 1, nrow
706 2004357 : call setUnifRand(RNG rand(irow) BOUNDS)
707 : end do
708 : #elif D2_ENABLED
709 : integer(IK) :: irow, nrow, icol, ncol
710 0 : nrow = size(rand, 1, IK)
711 0 : ncol = size(rand, 2, IK)
712 0 : do icol = 1, size(rand, 2, IK)
713 0 : do irow = 1, nrow
714 0 : call setUnifRand(RNG rand(irow, icol) BOUNDS)
715 : end do
716 : end do
717 : #elif D3_ENABLED
718 : integer(IK) :: irow, nrow, icol, ncol, idim, ndim
719 0 : nrow = size(rand, 1, IK)
720 0 : ncol = size(rand, 2, IK)
721 0 : ndim = size(rand, 3, IK)
722 0 : do idim = 1, ndim
723 0 : do icol = 1, ncol
724 0 : do irow = 1, nrow
725 0 : call setUnifRand(RNG rand(irow, icol, idim) BOUNDS)
726 : end do
727 : end do
728 : end do
729 : #else
730 : #error "Unrecognized interface."
731 : #endif
732 : #undef BOUNDS
733 :
734 : #else
735 : !%%%%%%%%%%%%%%%%%%%%%%%%
736 : #error "Unrecognized interface."
737 : !%%%%%%%%%%%%%%%%%%%%%%%%
738 : #endif
739 : #undef GET_ELEMENT
740 : #undef RNG
|
