Line data Source code
1 : !> \file mo_xor4096.f90
2 : !> \brief \copybrief mo_xor4096
3 : !> \details \copydetails mo_xor4096
4 :
5 : !> \brief XOR4096-based random number generator
6 : !> \details This module provides random number generator based on xor4096 algorithm.
7 : !> \authors Juliane Mai
8 : !> \date Nov 2011
9 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
10 : !! FORCES is released under the LGPLv3+ license \license_note
11 : module mo_xor4096
12 :
13 : ! -----------------------------------------------------------------------------
14 : ! The original version of this source code (without multiple streams, optional
15 : ! arguments and gaussian distributed RN) is under GNU General Public Licence
16 : ! xorgens.c
17 : ! Copyright (C) 2004 R. P. Brent.
18 : ! This program is free software; you can redistribute it and/or
19 : ! modify it under the terms of the GNU General Public License,
20 : ! version 2, June 1991, as published by the Free Software Foundation.
21 : ! For details see http://www.gnu.org/copyleft/gpl.html .
22 :
23 : ! Author: Richard P. Brent (random@rpbrent.co.uk)
24 : ! -----------------------------------------------------------------------------
25 :
26 : use mo_kind, only : i4, i8, sp, dp
27 :
28 : Implicit NONE
29 :
30 : PUBLIC :: n_save_state ! dimension of vector keeping the state of a stream
31 :
32 : PUBLIC :: get_timeseed ! Returns a seed dependend on time
33 : PUBLIC :: xor4096 ! Generates uniform distributed random number
34 : PUBLIC :: xor4096g ! Generates gaussian distributed random number
35 :
36 : ! ------------------------------------------------------------------
37 :
38 : !> \brief Generate seed for xor4096
39 :
40 : !> \details
41 : !! Function which returns a scalar or a vector of integers
42 : !! dependent on time.
43 : !! This returned values can be used a seed for initializing
44 : !! random number generators like xor4096.\n
45 : !!
46 : !! If the return variable is a vector, only the first entry
47 : !! depends on time while the others are just the entry before
48 : !! increased by 1000.
49 : !!
50 : !! \b Example
51 : !!
52 : !! Results in a scalar seed
53 : !! \code{.f90}
54 : !! integer(i4) :: seed
55 : !! call get_timeseed(seed)
56 : !! \endcode
57 : !!
58 : !! Results in a array of seeds
59 : !! \code{.f90}
60 : !! integer(i8), dimension(3) :: seed
61 : !! call get_timeseed(seed)
62 : !! \endcode
63 : !! e.g. `seed = (/ 327_i8, 1327_i8, 2327_i8 /)`
64 :
65 :
66 : !> \param[in] "integer(i4/i8), dimension(:) :: seed" Seed to be generated
67 :
68 :
69 : !> \author Juliane Mai
70 : !> \date Aug 2012
71 :
72 : INTERFACE get_timeseed
73 : MODULE PROCEDURE get_timeseed_i4_0d, get_timeseed_i4_1d, &
74 : get_timeseed_i8_0d, get_timeseed_i8_1d
75 : END INTERFACE get_timeseed
76 :
77 : ! ------------------------------------------------------------------
78 :
79 : !> \brief Uniform XOR4096-based RNG
80 :
81 : !> \details
82 : !! Generates a uniform distributed random number based on xor4096 algorithm proposed by
83 : !! Brent et.al (2006).\n
84 : !! ****************************************************************************************
85 : !! The original version of this source code (without multiple streams and
86 : !! optional arguments) is under GNU General Public Licence
87 : !! xorgens.c
88 : !! Copyright (C) 2004 R. P. Brent.
89 : !! This program is free software; you can redistribute it and/or
90 : !! modify it under the terms of the GNU General Public License,
91 : !! version 2, June 1991, as published by the Free Software Foundation.
92 : !! For details see http://www.gnu.org/copyleft/gpl.html .
93 : !! Author: Richard P. Brent (random@rpbrent.co.uk)
94 : !! ****************************************************************************************
95 : !! The period of the generator is\n
96 : !! (2^4096 - 1)*2^32 for single precision version and
97 : !! (2^4096 - 1)*2^64 for double precision version.\n
98 : !!
99 : !! The generator is based on bitwise XOR compositions of left- and right-shifted numbers.\n
100 : !!
101 : !! The generator has to be called once with a non-zero seed value which initializes a new
102 : !! random number stream. The subsequent calls are with seed=0 which returns the following
103 : !! numbers within the beforehand initialized stream. Since the random numbers are based on
104 : !! the initial seed, the precison of the seed (sp/dp) determines the precision of the
105 : !! returned random number (sp/dp).\n
106 : !!
107 : !! If one initialize the generator with an array of seeds, one initializes n independent
108 : !! streams of random numbers.
109 : !! Lets assume that the streams with seed:
110 : !! 1_sp is 10, 20, 30 ...\n
111 : !! 2_sp is 40, 50, 60 ...\n
112 : !! 3_sp is 70, 80, 90 ...\n
113 : !!
114 : !! What you get is:\n
115 : !! 1st call
116 : !! \code{.f90}
117 : !! call xor( (/ 1_SP, 2_SP, 3_SP /), RN ) --> RN = (/ 10, 40, 70 /)
118 : !! \endcode
119 : !! 2nd call
120 : !! \code{.f90}
121 : !! call xor( (/ 0_SP, 0_SP, 0_SP /), RN ) --> RN = (/ 20, 50, 80 /)
122 : !! \endcode
123 : !! 3rd call
124 : !! \code{.f90}
125 : !! call xor( (/ 0_SP, 0_SP, 0_SP /), RN ) --> RN = (/ 30, 60, 90 /)
126 : !! \endcode
127 : !!
128 : !! Since after every initialization one looses the old stream of random numbers,
129 : !! it might be necessary to switch between different streams. Therefore, one needs
130 : !! the optional save_state argument.
131 : !! What you get is:
132 : !! 1st call of 1st stream
133 : !! \code{.f90}
134 : !! call xor( 1_SP, RN, save_state=save_stream_1 ) RN = 10
135 : !! \endcode
136 : !! 2nd call of 1st stream
137 : !! \code{.f90}
138 : !! call xor( 0_SP, RN, save_state=save_stream_1 ) RN = 20
139 : !! \endcode
140 : !! 1st call of 2nd stream
141 : !! \code{.f90}
142 : !! call xor( 2_SP, RN, save_state=save_stream_2 ) RN = 40
143 : !! \endcode
144 : !! 3rd call of 1st stream
145 : !! \code{.f90}
146 : !! call xor( 0_SP, RN, save_state=save_stream_1 ) RN = 30
147 : !! \endcode
148 : !! Note: If you would have called 4 times without optional argument,
149 : !! you would have get 50 in the 4th call.
150 : !!
151 : !! \b Example
152 : !!
153 : !! I4 seeds are generating I4 or SP random numbers,
154 : !! I8 seeds are generating I8 or DP random numbers.
155 : !!
156 : !! Initializing:
157 : !! \code{.f90}
158 : !! real(SP) :: RN(3)
159 : !! seed = (/ 1_I4, 100_I4, 2_I4 /)
160 : !! call xor4096(seed,RN)
161 : !! print*, RN --> (/ 0.1_sp, 0.2_sp, 0.6_sp /)
162 : !! \endcode
163 : !!
164 : !! After initializing
165 : !! \code{.f90}
166 : !! seed = (/ 0_I4, 0_I4, 0_I4 /)
167 : !! call xor4096(seed,RN)
168 : !! print*, RN --> (/ 0.3_sp, 0.1_sp, 0.5_sp /)
169 : !! \endcode
170 : !!
171 : !! \b Literature
172 : !!
173 : !! 1. Brent RP. _Some long-period random number generators using shifts and xors_, 2010
174 : !! 2. Brent RP. _From Mersenne Primes to Rndom Number Generators_, 2006
175 : !! 3. L''Ecuyer P & Simard R. _ACM: TestU01: A C Library for Empirical Testing of
176 : !! Random Number Generators_, 2007
177 : !!
178 : !> \param[in] "integer(i4/i8) :: seed/seed(:)" Value or 1D-array with non-zero seeds for
179 : !! initialization or zero for subsequent callsyy
180 : !> \param[out] "integer(i4/i8)/real(sp/dp) :: RN/RN(size(seed))"
181 : !! uniform distributed random number with
182 : !! interval:\n
183 : !! - i4: (-2^31,2^31-1)
184 : !! - i8: (-2^63,2^63-1)
185 : !! - sp: (0.0_sp, 1.0_sp)
186 : !! - dp: (0.0_dp, 1.0_dp)
187 : !> \param[inout] "integer(i4/i8), dimension(size(seed), n_save_state), opional :: save_state"
188 : !! Array carrying state of random number stream
189 : !! this should be used if several streams are used, i.e.
190 : !! each stream has its own save_state.
191 :
192 : !! \note
193 : !! - The dimension of the optional argument can be set using the public module parameter n_save_state.
194 : !! - If random numbers are in single precision (sp), one needs seed and save_state in (i4).
195 : !! - If random numbers are in double precision (dp), one needs seed and save_state in (i8).
196 :
197 : !> \author Juliane Mai
198 : !> \date Nov 2011
199 :
200 : INTERFACE xor4096
201 : MODULE PROCEDURE xor4096s_0d, xor4096s_1d, xor4096f_0d, xor4096f_1d, &
202 : xor4096l_0d, xor4096l_1d, xor4096d_0d, xor4096d_1d
203 : END INTERFACE xor4096
204 :
205 : ! ------------------------------------------------------------------
206 :
207 : !> \brief Gaussian XOR4096-based RNG
208 :
209 : !> \details
210 : !! Generates a gaussian distributed random number. First, a uniform distributed random
211 : !! number based on xor4096 algorithm is generated. Second, this number is transformed
212 : !! using the Polar method of Box-Mueller-transform to calculate the gaussian distributed
213 : !! numbers.\n
214 : !!
215 : !! ****************************************************************************************
216 : !! The original version of this source code (without multiple streams,
217 : !! optional arguments and gaussian distributed RN) is under GNU General Public Licence
218 : !! xorgens.c
219 : !! Copyright (C) 2004 R. P. Brent.
220 : !! This program is free software; you can redistribute it and/or
221 : !! modify it under the terms of the GNU General Public License,
222 : !! version 2, June 1991, as published by the Free Software Foundation.
223 : !! For details see http://www.gnu.org/copyleft/gpl.html .
224 : !!
225 : !! Author: Richard P. Brent (random@rpbrent.co.uk)
226 : !! ****************************************************************************************
227 : !!
228 : !! The generator has to be called once with a non-zero seed value which initializes a new
229 : !! random number stream. The subsequent calls are with seed=0 which returns the following
230 : !! numbers within the beforehand initialized stream. Since the random numbers are based on
231 : !! the initial seed, the precison of the seed (sp/dp) determines the precision of the
232 : !! returned random number (sp/dp).\n
233 : !!
234 : !! The returned values are gaussian distributed with mean 0 and variance 1.\n
235 : !!
236 : !! The polar method of the Box-Mueller transform transforms two uniform distributed
237 : !! random numbers \f$u_1\f$ and \f$u_2\f$ into two gaussian distributed random numbers \f$x_1\f$ and \f$x_2\f$.
238 : !! First, two uniform numbers a1, a2 distributed between (-1,1) are calculated:\n
239 : !! \f[a_1 = 2u_1 - 1\f]
240 : !! \f[a_2 = 2u_2 - 1\f]
241 : !! Second, \f$q = a_1^2 + a_2^2\f$ is calculated. If (\f$q = 0\f$) or (\f$q > 1\f$) one has to generate
242 : !! a new \f$x_1\f$ and \f$x_2\f$, since a divison by q is needed afterwards and q needs to be distributed
243 : !! within the unit circle.
244 : !! Third,
245 : !! \f[p = \sqrt{-2 \frac{\ln{q}}{q}}\f]
246 : !! is calculated.
247 : !! The numbers \f$z_1\f$ and \f$z_2\f$ with
248 : !! \f[z_1 = a_1 p\f]
249 : !! \f[z_2 = a_2 p\f]
250 : !! are then gaussian distributed with mean 0 and variance 1.
251 : !!
252 : !! \b Example
253 : !!
254 : !! See example for xor4096.
255 : !!
256 : !! \b Literature
257 : !!
258 : !! 1. Brent RP. _Some long-period random number generators using shifts and xors_, 2010
259 : !! 2. Brent RP. _From Mersenne Primes to Rndom Number Generators_, 2006
260 : !! 3. L''Ecuyer P & Simard R. _ACM: TestU01: A C Library for Empirical Testing of
261 : !! Random Number Generators_, 2007
262 : !! 4. http://en.wikipedia.org/wiki/Marsaglia_polar_method
263 : !! 5. http://de.wikipedia.org/wiki/Polar-Methode
264 : !!
265 : !> \param[in] "integer(i4/i8) :: seed/seed(:)" Value or 1D-array with non-zero seeds for
266 : !! initialization or zero for subsequent calls
267 : !> \param[out] "real(sp/dp) :: RN/RN(size(seed))"
268 : !! Gaussian distributed random number with
269 : !! interval:\n
270 : !! - sp: RN ~ N(0.0_sp, 1.0_sp)
271 : !! - dp: RN ~ N(0.0_dp, 1.0_dp)
272 : !> \param[inout] "integer(i4/i8), dimension(size(seed),n_save_state), optional :: save_state"
273 : !! Array carrying state of random number stream
274 : !! this should be used if several streams are used, i.e.
275 : !! each stream has its own save_state.
276 :
277 : !> \note
278 : !! - The dimension of the optional argument can be set using the public module parameter n_save_state.
279 : !! - If random numbers are in single precision (sp), one needs seed and save_state in (i4).
280 : !! - If random numbers are in double precision (dp), one needs seed and save_state in (i8).
281 :
282 : !> \author Juliane Mai
283 : !> \date Nov 2011
284 : !> \date Feb 2013
285 : !! - all optionals combined in save_state
286 :
287 : INTERFACE xor4096g
288 : MODULE PROCEDURE xor4096gf_0d, xor4096gf_1d, xor4096gd_0d, xor4096gd_1d
289 : END INTERFACE xor4096g
290 :
291 : ! ------------------------------------------------------------------
292 :
293 : PRIVATE
294 :
295 : !> Dimension of vector saving the state of a stream
296 : integer(i4), parameter :: n_save_state = 132_i4
297 :
298 : ! ------------------------------------------------------------------
299 :
300 : CONTAINS
301 :
302 : ! ------------------------------------------------------------------
303 :
304 0 : subroutine get_timeseed_i4_0d(seed)
305 :
306 : implicit none
307 : integer(i4), intent(inout) :: seed
308 :
309 : ! local variables
310 : integer(i4), dimension(8) :: time_array
311 :
312 0 : call date_and_time(values = time_array)
313 : seed = &
314 : time_array(5) * 3600000_i4 + & ! hour
315 : time_array(6) * 60000_i4 + & ! minutes
316 : time_array(7) * 1000_i4 + & ! seconds
317 0 : time_array(8) * 1_i4 ! milliseconds
318 :
319 0 : end subroutine get_timeseed_i4_0d
320 :
321 0 : subroutine get_timeseed_i4_1d(seed)
322 :
323 : implicit none
324 : integer(i4), dimension(:), intent(inout) :: seed
325 :
326 : ! local variables
327 : integer(i4), dimension(8) :: time_array
328 : integer(i4) :: i
329 :
330 0 : call date_and_time(values = time_array)
331 0 : seed(1) = &
332 : time_array(5) * 3600000_i4 + & ! hour
333 : time_array(6) * 60000_i4 + & ! minutes
334 : time_array(7) * 1000_i4 + & ! seconds
335 0 : time_array(8) * 1_i4 ! milliseconds
336 0 : do i = 2, size(seed)
337 0 : seed(i) = seed(i - 1) + 1000_i4
338 : end do
339 :
340 0 : end subroutine get_timeseed_i4_1d
341 :
342 0 : subroutine get_timeseed_i8_0d(seed)
343 :
344 : implicit none
345 : integer(i8), intent(inout) :: seed
346 :
347 : ! local variables
348 : integer(i4), dimension(8) :: time_array
349 :
350 0 : call date_and_time(values = time_array)
351 : seed = &
352 : int(time_array(5), i8) * 3600000_i8 + & ! hour
353 : int(time_array(6), i8) * 60000_i8 + & ! minutes
354 : int(time_array(7), i8) * 1000_i8 + & ! seconds
355 0 : int(time_array(8), i8) * 1_i8 ! milliseconds
356 :
357 0 : end subroutine get_timeseed_i8_0d
358 :
359 1 : subroutine get_timeseed_i8_1d(seed)
360 :
361 : implicit none
362 : integer(i8), dimension(:), intent(inout) :: seed
363 :
364 : ! local variables
365 : integer(i4), dimension(8) :: time_array
366 : integer(i4) :: i
367 :
368 1 : call date_and_time(values = time_array)
369 1 : seed(1) = &
370 : int(time_array(5), i8) * 3600000_i8 + & ! hour
371 : int(time_array(6), i8) * 60000_i8 + & ! minutes
372 : int(time_array(7), i8) * 1000_i8 + & ! seconds
373 1 : int(time_array(8), i8) * 1_i8 ! milliseconds
374 3 : do i = 2, size(seed)
375 3 : seed(i) = seed(i - 1) + 1000_i8
376 : end do
377 :
378 0 : end subroutine get_timeseed_i8_1d
379 :
380 : ! ------------------------------------------------------------------
381 :
382 11 : subroutine xor4096s_0d(seed, SingleIntegerRN, save_state)
383 : implicit none
384 :
385 : integer(i4), intent(in) :: seed
386 : integer(i4), intent(out) :: SingleIntegerRN
387 : integer(i4), optional, dimension(n_save_state), intent(inout) :: save_state
388 :
389 : integer(i4) :: wlen, r, s, a, b, c, d
390 :
391 : integer(i4), save :: w
392 : integer(i4), save :: x(0 : 127) ! x(0) ... x(r-1)
393 : integer(i4) :: weyl = 1640531527_i4 !Z'61C88647' ! Hexadecimal notation
394 : integer(i4) :: t, v, tmp
395 : integer(i4), save :: i = -1 ! i<0 indicates first call
396 : integer(i4) :: k
397 :
398 11 : wlen = 32
399 11 : r = 128
400 11 : s = 95
401 11 : a = 17
402 11 : b = 12
403 11 : c = 13
404 11 : d = 15
405 :
406 11 : if (present(save_state) .and. (seed .eq. 0)) then
407 516 : x(0 : r - 1) = save_state(1 : r)
408 4 : i = save_state(r + 1)
409 4 : w = save_state(r + 2)
410 : end if
411 :
412 11 : If ((i .lt. 0) .or. (seed .ne. 0)) then ! Initialization necessary
413 3 : If (seed .ne. 0) then ! v must be nonzero
414 : v = seed
415 : else
416 0 : v = NOT(seed)
417 : end if
418 :
419 99 : do k = wlen, 1, -1 ! Avoid correlations for close seeds
420 : ! This recurrence has period of 2^32-1
421 96 : v = IEOR(v, ISHFT(v, 13))
422 96 : v = IEOR(v, ISHFT(v, -17))
423 99 : v = IEOR(v, ISHFT(v, 5))
424 : end do
425 :
426 : ! Initialize circular array
427 3 : w = v
428 387 : do k = 0, r - 1
429 : ! w = w + weyl
430 384 : if (w < 0_i4) then
431 191 : w = w + weyl
432 193 : else if ((huge(1_i4) - w) > weyl) then
433 46 : w = w + weyl
434 : else
435 147 : tmp = -(huge(1_i4) - w - weyl)
436 147 : w = tmp - huge(1_i4) - 2_i4
437 : endif
438 384 : v = IEOR(v, ISHFT(v, 13))
439 384 : v = IEOR(v, ISHFT(v, -17))
440 384 : v = IEOR(v, ISHFT(v, 5))
441 387 : x(k) = v + w
442 : end do
443 :
444 : ! Discard first 4*r results (Gimeno)
445 3 : i = r - 1
446 1539 : do k = 4 * r, 1, -1
447 1536 : i = IAND(i + 1, r - 1)
448 1536 : t = x(i)
449 1536 : v = x(IAND(i + (r - s), r - 1))
450 1536 : t = IEOR(t, ISHFT(t, a))
451 1536 : t = IEOR(t, ISHFT(t, -b))
452 1536 : v = IEOR(v, ISHFT(v, c))
453 1536 : v = IEOR(v, IEOR(t, ISHFT(v, -d)))
454 1539 : x(i) = v
455 : end do
456 : end if ! end of initialization
457 :
458 : ! Apart from initialization (above), this is the generator
459 11 : i = IAND(i + 1, r - 1)
460 11 : t = x(i)
461 11 : v = x(IAND(i + (r - s), r - 1))
462 11 : t = IEOR(t, ISHFT(t, a))
463 11 : t = IEOR(t, ISHFT(t, -b))
464 11 : v = IEOR(v, ISHFT(v, c))
465 11 : v = IEOR(v, IEOR(t, ISHFT(v, -d)))
466 11 : x(i) = v
467 :
468 11 : w = w + weyl
469 :
470 11 : SingleIntegerRN = v + w
471 :
472 11 : if(present(save_state)) then
473 645 : save_state(1 : r) = x(0 : r - 1)
474 5 : save_state(r + 1) = i
475 5 : save_state(r + 2) = w
476 15 : if ((r + 3) <= n_save_state) save_state(r + 3 : n_save_state) = 0
477 : end if
478 :
479 1 : end subroutine xor4096s_0d
480 :
481 : ! -----------------------------------------------------------------------------
482 :
483 5 : subroutine xor4096s_1d(seed, SingleIntegerRN, save_state)
484 : implicit none
485 :
486 : integer(i4), dimension(:), intent(in) :: seed
487 : integer(i4), dimension(size(seed, 1)), intent(out) :: SingleIntegerRN
488 : integer(i4), optional, dimension(size(seed, 1), n_save_state), intent(inout) :: save_state
489 :
490 : integer(i4) :: m
491 : integer(i4) :: wlen, r, s, a, b, c, d
492 : integer(i4) :: weyl = 1640531527_i4 !Z'61C88647' ! Hexadecimal notation
493 : integer(i4) :: k, j, tmp
494 5 : integer(i4), dimension(size(seed, 1)) :: t, v
495 : integer(i4), dimension(:, :), allocatable, save :: x ! x(0) ... x(r-1)
496 : integer(i4), dimension(:), allocatable, save :: i, w ! i<0 indicates first call
497 :
498 5 : wlen = 32
499 5 : r = 128
500 5 : s = 95
501 5 : a = 17
502 5 : b = 12
503 5 : c = 13
504 5 : d = 15
505 :
506 5 : m = size(seed, 1)
507 :
508 25 : if (any(seed .eq. 0_i4) .and. any(seed .ne. 0_i4)) then
509 0 : stop 'xor4096: seeds have to be eigther all 0 or all larger than 0'
510 : end if
511 :
512 17 : if (present(save_state) .and. all(seed .eq. 0_i4)) then
513 0 : if (allocated(x)) then
514 0 : if (size(x, 1) .ne. m) then
515 0 : deallocate(x)
516 0 : deallocate(i)
517 0 : deallocate(w)
518 0 : allocate(i(m))
519 0 : allocate(x(m, 0 : r - 1))
520 0 : allocate(w(m))
521 : end if
522 : end if
523 0 : if (.not. allocated(x)) then
524 0 : allocate(i(m))
525 0 : allocate(x(m, 0 : r - 1))
526 0 : allocate(w(m))
527 : end if
528 0 : x(:, 0 : r - 1) = save_state(:, 1 : r)
529 0 : i(:) = save_state(:, r + 1)
530 0 : w(:) = save_state(:, r + 2)
531 : end if
532 :
533 8 : if(all(seed .ne. 0_i4)) then
534 1 : if (allocated(x)) then
535 0 : deallocate(x)
536 0 : deallocate(i)
537 0 : deallocate(w)
538 : end if
539 :
540 3 : allocate(i(m))
541 5 : i = -1
542 3 : allocate(x(m, 0 : r - 1))
543 2 : allocate(w(m))
544 : end if
545 :
546 20 : Do j = 1, m
547 20 : If ((i(j) .lt. 0) .or. (seed(j) .ne. 0)) then ! Initialization necessary
548 3 : If (seed(j) .ne. 0) then ! v must be nonzero
549 3 : v(j) = seed(j)
550 : else
551 0 : v(j) = NOT(seed(j))
552 : end if
553 :
554 99 : do k = wlen, 1, -1 ! Avoid correlations for close seeds
555 : ! This recurrence has period of 2^32-1
556 96 : v(j) = IEOR(v(j), ISHFT(v(j), 13))
557 96 : v(j) = IEOR(v(j), ISHFT(v(j), -17))
558 99 : v(j) = IEOR(v(j), ISHFT(v(j), 5))
559 : end do
560 :
561 : ! Initialize circular array
562 3 : w(j) = v(j)
563 387 : do k = 0, r - 1
564 : ! w(j) = w(j) + weyl
565 384 : if (w(j) < 0_i4) then
566 192 : w(j) = w(j) + weyl
567 192 : else if ((huge(1_i4) - w(j)) > weyl) then
568 45 : w(j) = w(j) + weyl
569 : else
570 147 : tmp = -(huge(1_i4) - w(j) - weyl)
571 147 : w(j) = tmp - huge(1_i4) - 2_i4
572 : endif
573 384 : v(j) = IEOR(v(j), ISHFT(v(j), 13))
574 384 : v(j) = IEOR(v(j), ISHFT(v(j), -17))
575 384 : v(j) = IEOR(v(j), ISHFT(v(j), 5))
576 387 : x(j, k) = v(j) + w(j)
577 : end do
578 :
579 : ! Discard first 4*r results (Gimeno)
580 3 : i(j) = r - 1
581 1539 : do k = 4 * r, 1, -1
582 1536 : i(j) = IAND(i(j) + 1, r - 1)
583 1536 : t(j) = x(j, i(j))
584 1536 : v(j) = x(j, IAND(i(j) + (r - s), r - 1))
585 1536 : t(j) = IEOR(t(j), ISHFT(t(j), a))
586 1536 : t(j) = IEOR(t(j), ISHFT(t(j), -b))
587 1536 : v(j) = IEOR(v(j), ISHFT(v(j), c))
588 1536 : v(j) = IEOR(v(j), IEOR(t(j), ISHFT(v(j), -d)))
589 1539 : x(j, i(j)) = v(j)
590 : end do
591 : end if ! end of initialization
592 : end do
593 :
594 : ! Apart from initialization (above), this is the generator
595 :
596 20 : do j = 1, m
597 15 : i(j) = IAND(i(j) + 1, r - 1)
598 15 : t(j) = x(j, i(j))
599 15 : v(j) = x(j, IAND(i(j) + (r - s), r - 1))
600 15 : t(j) = IEOR(t(j), ISHFT(t(j), a))
601 15 : t(j) = IEOR(t(j), ISHFT(t(j), -b))
602 15 : v(j) = IEOR(v(j), ISHFT(v(j), c))
603 15 : v(j) = IEOR(v(j), IEOR(t(j), ISHFT(v(j), -d)))
604 15 : x(j, i(j)) = v(j)
605 :
606 20 : w(j) = w(j) + weyl
607 : end do
608 :
609 20 : SingleIntegerRN = v + w
610 :
611 5 : if(present(save_state)) then
612 0 : save_state(:, 1 : r) = x(:, 0 : r - 1)
613 0 : save_state(:, r + 1) = i(:)
614 0 : save_state(:, r + 2) = w(:)
615 0 : if ((r + 3) <= n_save_state) save_state(:, r + 3 : n_save_state) = 0
616 : end if
617 :
618 16 : end subroutine xor4096s_1d
619 :
620 : ! -----------------------------------------------------------------------------
621 :
622 0 : subroutine xor4096f_0d(seed, SingleRealRN, save_state)
623 :
624 : implicit none
625 :
626 : integer(i4), intent(in) :: seed
627 : real(SP), intent(out) :: SingleRealRN
628 : integer(i4), optional, dimension(n_save_state), intent(inout) :: save_state
629 :
630 : integer(i4) :: wlen, r, s, a, b, c, d
631 : integer(i4), save :: w
632 : integer(i4), save :: x(0 : 127) ! x(0) ... x(r-1)
633 : integer(i4) :: weyl = 1640531527_i4 !Z'61C88647' ! Hexadecimal notation
634 : integer(i4) :: t, v, tmp
635 : integer(i4), save :: i = -1 ! i<0 indicates first call
636 : integer(i4) :: k
637 :
638 : real(SP) :: t24 = 1.0_SP / 16777216.0_SP ! = 0.5^24 = 1/2^24
639 :
640 : !$omp threadprivate(x,i,w)
641 : !
642 :
643 : ! produces a 24bit Integer Random Number (0...16777216) and
644 : ! scales it afterwards to (0.0,1.0)
645 :
646 0 : wlen = 32
647 0 : r = 128
648 0 : s = 95
649 0 : a = 17
650 0 : b = 12
651 0 : c = 13
652 0 : d = 15
653 :
654 0 : if (present(save_state) .and. (seed .eq. 0)) then
655 0 : x(0 : r - 1) = save_state(1 : r)
656 0 : i = save_state(r + 1)
657 0 : w = save_state(r + 2)
658 : end if
659 :
660 0 : If ((i .lt. 0) .or. (seed .ne. 0)) then ! Initialization necessary
661 0 : If (seed .ne. 0) then ! v must be nonzero
662 : v = seed
663 : else
664 0 : v = NOT(seed)
665 : end if
666 :
667 0 : do k = wlen, 1, -1 ! Avoid correlations for close seeds
668 : ! This recurrence has period of 2^32-1
669 0 : v = IEOR(v, ISHFT(v, 13))
670 0 : v = IEOR(v, ISHFT(v, -17))
671 0 : v = IEOR(v, ISHFT(v, 5))
672 : end do
673 :
674 : ! Initialize circular array
675 0 : w = v
676 0 : do k = 0, r - 1
677 : ! w = w + weyl
678 0 : if (w < 0_i4) then
679 0 : w = w + weyl
680 0 : else if ((huge(1_i4) - w) > weyl) then
681 0 : w = w + weyl
682 : else
683 0 : tmp = -(huge(1_i4) - w - weyl)
684 0 : w = tmp - huge(1_i4) - 2_i4
685 : endif
686 0 : v = IEOR(v, ISHFT(v, 13))
687 0 : v = IEOR(v, ISHFT(v, -17))
688 0 : v = IEOR(v, ISHFT(v, 5))
689 0 : x(k) = v + w
690 : end do
691 :
692 : ! Discard first 4*r results (Gimeno)
693 0 : i = r - 1
694 0 : do k = 4 * r, 1, -1
695 0 : i = IAND(i + 1, r - 1)
696 0 : t = x(i)
697 0 : v = x(IAND(i + (r - s), r - 1))
698 0 : t = IEOR(t, ISHFT(t, a))
699 0 : t = IEOR(t, ISHFT(t, -b))
700 0 : v = IEOR(v, ISHFT(v, c))
701 0 : v = IEOR(v, IEOR(t, ISHFT(v, -d)))
702 0 : x(i) = v
703 : end do
704 : end if ! end of initialization
705 :
706 : ! Apart from initialization (above), this is the generator
707 : v = 0_i4
708 0 : Do While (v .eq. 0_i4)
709 0 : i = IAND(i + 1, r - 1)
710 0 : t = x(i)
711 0 : v = x(IAND(i + (r - s), r - 1))
712 0 : t = IEOR(t, ISHFT(t, a))
713 0 : t = IEOR(t, ISHFT(t, -b))
714 0 : v = IEOR(v, ISHFT(v, c))
715 0 : v = IEOR(v, IEOR(t, ISHFT(v, -d)))
716 0 : x(i) = v
717 0 : w = w + weyl
718 0 : v = v + w
719 0 : v = ISHFT(v, -8)
720 : End Do
721 :
722 0 : SingleRealRN = t24 * v
723 :
724 0 : if(present(save_state)) then
725 0 : save_state(1 : r) = x(0 : r - 1)
726 0 : save_state(r + 1) = i
727 0 : save_state(r + 2) = w
728 0 : if ((r + 3) <= n_save_state) save_state(r + 3 : n_save_state) = 0
729 : end if
730 :
731 : !
732 :
733 5 : end subroutine xor4096f_0d
734 :
735 : ! -----------------------------------------------------------------------------
736 :
737 0 : subroutine xor4096f_1d(seed, SingleRealRN, save_state)
738 :
739 : implicit none
740 :
741 : integer(i4), dimension(:), intent(in) :: seed
742 : real(SP), dimension(size(seed)), intent(out) :: SingleRealRN
743 : integer(i4), optional, dimension(size(seed, 1), n_save_state), intent(inout) :: save_state
744 :
745 : integer(i4) :: m
746 : integer(i4) :: wlen, r, s, a, b, c, d
747 : integer(i4) :: weyl = 1640531527_i4 !Z'61C88647' = Hexadecimal notation
748 : integer(i4) :: k, j, tmp
749 : real(SP), save :: t24 = 1.0_SP / 16777216.0_SP ! = 0.5^24 = 1/2^24
750 0 : integer(i4), dimension(size(seed)) :: t, v
751 : integer(i4), dimension(:, :), allocatable, save :: x ! x(0) ... x(r-1)
752 : integer(i4), dimension(:), allocatable, save :: i, w ! i<0 indicates first call
753 :
754 : ! produces a 24bit Integer Random Number (0...16777216) and
755 : ! scales it afterwards to (0.0,1.0)
756 :
757 : !$omp threadprivate(x,i,w)
758 :
759 0 : wlen = 32
760 0 : r = 128
761 0 : s = 95
762 0 : a = 17
763 0 : b = 12
764 0 : c = 13
765 0 : d = 15
766 :
767 0 : m = size(seed, 1)
768 :
769 0 : if (any(seed .eq. 0_i4) .and. any(seed .ne. 0_i4)) then
770 0 : stop 'xor4096: seeds have to be eigther all 0 or all larger than 0'
771 : end if
772 :
773 0 : if (present(save_state) .and. all(seed .eq. 0_i4)) then
774 0 : if (allocated(x)) then
775 0 : if (size(x, 1) .ne. m) then
776 0 : deallocate(x)
777 0 : deallocate(i)
778 0 : deallocate(w)
779 0 : allocate(i(m))
780 0 : allocate(x(m, 0 : r - 1))
781 0 : allocate(w(m))
782 : end if
783 : end if
784 0 : if (.not. allocated(x)) then
785 0 : allocate(i(m))
786 0 : allocate(x(m, 0 : r - 1))
787 0 : allocate(w(m))
788 : end if
789 0 : x(:, 0 : r - 1) = save_state(:, 1 : r)
790 0 : i(:) = save_state(:, r + 1)
791 0 : w(:) = save_state(:, r + 2)
792 : end if
793 :
794 0 : if(all(seed .ne. 0_i4)) then
795 0 : if (allocated(x)) then
796 0 : deallocate(x)
797 0 : deallocate(i)
798 0 : deallocate(w)
799 : end if
800 :
801 0 : allocate(i(m))
802 0 : i = -1
803 0 : allocate(x(m, 0 : r - 1))
804 0 : allocate(w(m))
805 : end if
806 :
807 0 : Do j = 1, m !Loop over every stream
808 0 : If ((i(j) .lt. 0) .or. (seed(j) .ne. 0)) then ! Initialization necessary
809 0 : If (seed(j) .ne. 0) then ! v must be nonzero
810 0 : v(j) = seed(j)
811 : else
812 0 : v(j) = NOT(seed(j))
813 : end if
814 :
815 0 : do k = wlen, 1, -1 ! Avoid correlations for close seeds
816 : ! This recurrence has period of 2^32-1
817 0 : v(j) = IEOR(v(j), ISHFT(v(j), 13))
818 0 : v(j) = IEOR(v(j), ISHFT(v(j), -17))
819 0 : v(j) = IEOR(v(j), ISHFT(v(j), 5))
820 : end do
821 :
822 : ! Initialize circular array
823 0 : w(j) = v(j)
824 0 : do k = 0, r - 1
825 : ! w(j) = w(j) + weyl
826 0 : if (w(j) < 0_i4) then
827 0 : w(j) = w(j) + weyl
828 0 : else if ((huge(1_i4) - w(j)) > weyl) then
829 0 : w(j) = w(j) + weyl
830 : else
831 0 : tmp = -(huge(1_i4) - w(j) - weyl)
832 0 : w(j) = tmp - huge(1_i4) - 2_i4
833 : endif
834 0 : v(j) = IEOR(v(j), ISHFT(v(j), 13))
835 0 : v(j) = IEOR(v(j), ISHFT(v(j), -17))
836 0 : v(j) = IEOR(v(j), ISHFT(v(j), 5))
837 0 : x(j, k) = v(j) + w(j)
838 : end do
839 :
840 : ! Discard first 4*r results (Gimeno)
841 0 : i(j) = r - 1
842 0 : do k = 4 * r, 1, -1
843 0 : i(j) = IAND(i(j) + 1, r - 1)
844 0 : t(j) = x(j, i(j))
845 0 : v(j) = x(j, IAND(i(j) + (r - s), r - 1))
846 0 : t(j) = IEOR(t(j), ISHFT(t(j), a))
847 0 : t(j) = IEOR(t(j), ISHFT(t(j), -b))
848 0 : v(j) = IEOR(v(j), ISHFT(v(j), c))
849 0 : v(j) = IEOR(v(j), IEOR(t(j), ISHFT(v(j), -d)))
850 0 : x(j, i(j)) = v(j)
851 : end do
852 : end if ! end of initialization
853 : end do
854 :
855 : ! Apart from initialization (above), this is the generator
856 0 : v = 0_i4
857 0 : Do j = 1, m
858 0 : Do While (v(j) .eq. 0_i4)
859 0 : i(j) = IAND(i(j) + 1, r - 1)
860 0 : t(j) = x(j, i(j))
861 0 : v(j) = x(j, IAND(i(j) + (r - s), r - 1))
862 0 : t(j) = IEOR(t(j), ISHFT(t(j), a))
863 0 : t(j) = IEOR(t(j), ISHFT(t(j), -b))
864 0 : v(j) = IEOR(v(j), ISHFT(v(j), c))
865 0 : v(j) = IEOR(v(j), IEOR(t(j), ISHFT(v(j), -d)))
866 0 : x(j, i(j)) = v(j)
867 0 : w(j) = w(j) + weyl
868 0 : v(j) = v(j) + w(j)
869 0 : v(j) = ISHFT(v(j), -8)
870 : End Do
871 : End Do
872 :
873 0 : SingleRealRN = t24 * v
874 :
875 0 : if(present(save_state)) then
876 0 : save_state(:, 1 : r) = x(:, 0 : r - 1)
877 0 : save_state(:, r + 1) = i(:)
878 0 : save_state(:, r + 2) = w(:)
879 0 : if ((r + 3) <= n_save_state) save_state(:, r + 3 : n_save_state) = 0
880 : end if
881 :
882 0 : end subroutine xor4096f_1d
883 :
884 : ! -----------------------------------------------------------------------------
885 :
886 0 : subroutine xor4096l_0d(seed, DoubleIntegerRN, save_state)
887 :
888 : implicit none
889 :
890 : integer(i8), intent(in) :: seed
891 : integer(i8), intent(out) :: DoubleIntegerRN
892 : integer(i8), optional, dimension(n_save_state), intent(inout) :: save_state
893 :
894 : integer(i8) :: wlen, r, s, a, b, c, d
895 : integer(i8), save :: w
896 : integer(i8), save :: x(0 : 63) ! x(0) ... x(r-1)
897 : integer(i8) :: weyl = 7046029254386353131_i8
898 : integer(i8) :: t, v, tmp
899 : integer(i8), save :: i = -1 ! i<0 indicates first call
900 : integer(i8) :: k
901 :
902 : !$omp threadprivate(x,i,w)
903 :
904 0 : wlen = 64_i8
905 0 : r = 64_i8
906 0 : s = 53_i8
907 0 : a = 33_i8
908 0 : b = 26_i8
909 0 : c = 27_i8
910 0 : d = 29_i8
911 :
912 0 : if (present(save_state) .and. (seed .eq. 0_i8)) then
913 0 : x(0 : r - 1) = save_state(1 : r)
914 0 : i = save_state(r + 1)
915 0 : w = save_state(r + 2)
916 : end if
917 :
918 0 : If ((i .lt. 0) .or. (seed .ne. 0)) then ! Initialization necessary
919 0 : If (seed .ne. 0) then ! v must be nonzero
920 : v = seed
921 : else
922 0 : v = NOT(seed)
923 : end if
924 :
925 0 : do k = wlen, 1, -1 ! Avoid correlations for close seeds
926 : ! This recurrence has period of 2^64-1
927 0 : v = IEOR(v, ISHFT(v, 7))
928 0 : v = IEOR(v, ISHFT(v, -9))
929 : end do
930 :
931 : ! Initialize circular array
932 0 : w = v
933 0 : do k = 0, r - 1
934 : ! w = w + weyl
935 0 : if (w < 0_i8) then
936 0 : w = w + weyl
937 0 : else if ((huge(1_i8) - w) > weyl) then
938 0 : w = w + weyl
939 : else
940 0 : tmp = -(huge(1_i8) - w - weyl)
941 0 : w = tmp - huge(1_i8) - 2_i8
942 : endif
943 0 : v = IEOR(v, ISHFT(v, 7))
944 0 : v = IEOR(v, ISHFT(v, -9))
945 0 : x(k) = v + w
946 : end do
947 :
948 : ! Discard first 4*r results (Gimeno)
949 0 : i = r - 1
950 0 : do k = 4 * r, 1, -1
951 0 : i = IAND(i + 1, r - 1)
952 0 : t = x(i)
953 0 : v = x(IAND(i + (r - s), r - 1))
954 0 : t = IEOR(t, ISHFT(t, a))
955 0 : t = IEOR(t, ISHFT(t, -b))
956 0 : v = IEOR(v, ISHFT(v, c))
957 0 : v = IEOR(v, IEOR(t, ISHFT(v, -d)))
958 0 : x(i) = v
959 : end do
960 : end if ! end of initialization
961 :
962 : ! Apart from initialization (above), this is the generator
963 0 : i = IAND(i + 1, r - 1)
964 0 : t = x(i)
965 0 : v = x(IAND(i + (r - s), r - 1))
966 0 : t = IEOR(t, ISHFT(t, a))
967 0 : t = IEOR(t, ISHFT(t, -b))
968 0 : v = IEOR(v, ISHFT(v, c))
969 0 : v = IEOR(v, IEOR(t, ISHFT(v, -d)))
970 0 : x(i) = v
971 :
972 0 : w = w + weyl
973 :
974 0 : DoubleIntegerRN = v + w
975 :
976 0 : if(present(save_state)) then
977 0 : save_state(1 : r) = x(0 : r - 1)
978 0 : save_state(r + 1) = i
979 0 : save_state(r + 2) = w
980 0 : if ((r + 3) <= n_save_state) save_state(r + 3 : n_save_state) = 0
981 : end if
982 :
983 0 : end subroutine xor4096l_0d
984 :
985 : ! -----------------------------------------------------------------------------
986 :
987 0 : subroutine xor4096l_1d(seed, DoubleIntegerRN, save_state)
988 :
989 : implicit none
990 :
991 : integer(i8), dimension(:), intent(in) :: seed
992 : integer(i8), dimension(size(seed, 1)), intent(out) :: DoubleIntegerRN
993 : integer(i8), optional, dimension(size(seed, 1), n_save_state), intent(inout) :: save_state
994 :
995 : integer(i4) :: m
996 : integer(i8) :: wlen, r, s, a, b, c, d
997 : integer(i8) :: weyl = 7046029254386353131_i8
998 : integer(i8) :: k, j, tmp
999 0 : integer(i8), dimension(size(seed)) :: t, v
1000 : integer(i8), dimension(:, :), allocatable, save :: x ! x(0) ... x(r-1)
1001 : integer(i8), dimension(:), allocatable, save :: i, w ! i<0 indicates first call
1002 :
1003 : !$omp threadprivate(x,i,w)
1004 :
1005 0 : wlen = 64_i8
1006 0 : r = 64_i8
1007 0 : s = 53_i8
1008 0 : a = 33_i8
1009 0 : b = 26_i8
1010 0 : c = 27_i8
1011 0 : d = 29_i8
1012 :
1013 0 : m = size(seed, 1)
1014 :
1015 0 : if (any(seed .eq. 0_i8) .and. any(seed .ne. 0_i8)) then
1016 0 : stop 'xor4096: seeds have to be eigther all 0 or all larger than 0'
1017 : end if
1018 :
1019 0 : if (present(save_state) .and. all(seed .eq. 0_i4)) then
1020 0 : if (allocated(x)) then
1021 0 : if (size(x, 1) .ne. m) then
1022 0 : deallocate(x)
1023 0 : deallocate(i)
1024 0 : deallocate(w)
1025 0 : allocate(i(m))
1026 0 : allocate(x(m, 0 : r - 1))
1027 0 : allocate(w(m))
1028 : end if
1029 : end if
1030 0 : if (.not. allocated(x)) then
1031 0 : allocate(i(m))
1032 0 : allocate(x(m, 0 : r - 1))
1033 0 : allocate(w(m))
1034 : end if
1035 0 : x(:, 0 : r - 1) = save_state(:, 1 : r)
1036 0 : i(:) = save_state(:, r + 1)
1037 0 : w(:) = save_state(:, r + 2)
1038 : end if
1039 :
1040 0 : if(all(seed .ne. 0_i4)) then
1041 0 : if (allocated(x)) then
1042 0 : deallocate(x)
1043 0 : deallocate(i)
1044 0 : deallocate(w)
1045 : end if
1046 :
1047 0 : allocate(i(m))
1048 0 : i = -1
1049 0 : allocate(x(m, 0 : r - 1))
1050 0 : allocate(w(m))
1051 : end if
1052 :
1053 0 : Do j = 1, m
1054 0 : If ((i(j) .lt. 0) .or. (seed(j) .ne. 0)) then ! Initialization necessary
1055 0 : If (seed(j) .ne. 0) then ! v must be nonzero
1056 0 : v(j) = seed(j)
1057 : else
1058 0 : v(j) = NOT(seed(j))
1059 : end if
1060 :
1061 0 : do k = wlen, 1, -1 ! Avoid correlations for close seeds
1062 : ! This recurrence has period of 2^64-1
1063 0 : v(j) = IEOR(v(j), ISHFT(v(j), 7))
1064 0 : v(j) = IEOR(v(j), ISHFT(v(j), -9))
1065 : end do
1066 :
1067 : ! Initialize circular array
1068 0 : w(j) = v(j)
1069 0 : do k = 0, r - 1
1070 : ! w(j) = w(j) + weyl
1071 0 : if (w(j) < 0_i8) then
1072 0 : w(j) = w(j) + weyl
1073 0 : else if ((huge(1_i8) - w(j)) > weyl) then
1074 0 : w(j) = w(j) + weyl
1075 : else
1076 0 : tmp = -(huge(1_i8) - w(j) - weyl)
1077 0 : w(j) = tmp - huge(1_i8) - 2_i8
1078 : endif
1079 0 : v(j) = IEOR(v(j), ISHFT(v(j), 7))
1080 0 : v(j) = IEOR(v(j), ISHFT(v(j), -9))
1081 0 : x(j, k) = v(j) + w(j)
1082 : end do
1083 :
1084 : ! Discard first 4*r results (Gimeno)
1085 0 : i(j) = r - 1
1086 0 : do k = 4 * r, 1, -1
1087 0 : i(j) = IAND(i(j) + 1, r - 1)
1088 0 : t(j) = x(j, i(j))
1089 0 : v(j) = x(j, IAND(i(j) + (r - s), r - 1))
1090 0 : t(j) = IEOR(t(j), ISHFT(t(j), a))
1091 0 : t(j) = IEOR(t(j), ISHFT(t(j), -b))
1092 0 : v(j) = IEOR(v(j), ISHFT(v(j), c))
1093 0 : v(j) = IEOR(v(j), IEOR(t(j), ISHFT(v(j), -d)))
1094 0 : x(j, i(j)) = v(j)
1095 : end do
1096 : end if ! end of initialization
1097 : end do
1098 :
1099 : ! Apart from initialization (above), this is the generator
1100 0 : do j = 1, m
1101 0 : i(j) = IAND(i(j) + 1, r - 1)
1102 0 : t(j) = x(j, i(j))
1103 0 : v(j) = x(j, IAND(i(j) + (r - s), r - 1))
1104 0 : t(j) = IEOR(t(j), ISHFT(t(j), a))
1105 0 : t(j) = IEOR(t(j), ISHFT(t(j), -b))
1106 0 : v(j) = IEOR(v(j), ISHFT(v(j), c))
1107 0 : v(j) = IEOR(v(j), IEOR(t(j), ISHFT(v(j), -d)))
1108 0 : x(j, i(j)) = v(j)
1109 :
1110 0 : w(j) = w(j) + weyl
1111 : end do
1112 :
1113 0 : DoubleIntegerRN = v + w
1114 :
1115 0 : if(present(save_state)) then
1116 0 : save_state(:, 1 : r) = x(:, 0 : r - 1)
1117 0 : save_state(:, r + 1) = i(:)
1118 0 : save_state(:, r + 2) = w(:)
1119 0 : if ((r + 3) <= n_save_state) save_state(:, r + 3 : n_save_state) = 0
1120 : end if
1121 :
1122 0 : end subroutine xor4096l_1d
1123 :
1124 : ! -----------------------------------------------------------------------------
1125 :
1126 6673953 : subroutine xor4096d_0d(seed, DoubleRealRN, save_state)
1127 :
1128 : implicit none
1129 :
1130 : integer(i8), intent(in) :: seed
1131 : real(DP), intent(out) :: DoubleRealRN
1132 : integer(i8), optional, dimension(n_save_state), intent(inout) :: save_state
1133 :
1134 : integer(i8) :: wlen, r, s, a, b, c, d
1135 :
1136 : integer(i8), save :: w
1137 : integer(i8), save :: x(0 : 63) ! x(0) ... x(r-1)
1138 : integer(i8) :: weyl = 7046029254386353131_i8
1139 : integer(i8) :: t, v, tmp
1140 : integer(i8), save :: i = -1 ! i<0 indicates first call
1141 : integer(i8) :: k
1142 :
1143 : real(DP) :: t53 = 1.0_DP / 9007199254740992.0_DP ! = 0.5^53 = 1/2^53
1144 :
1145 : !$omp threadprivate(x,i,w)
1146 :
1147 : ! produces a 53bit Integer Random Number (0...9 007 199 254 740 992) and
1148 : ! scales it afterwards to (0.0,1.0)
1149 :
1150 6673953 : wlen = 64_i8
1151 6673953 : r = 64_i8
1152 6673953 : s = 53_i8
1153 6673953 : a = 33_i8
1154 6673953 : b = 26_i8
1155 6673953 : c = 27_i8
1156 6673953 : d = 29_i8
1157 :
1158 6673953 : if (present(save_state) .and. (seed .eq. 0_i8)) then
1159 291121935 : x(0 : r - 1) = save_state(1 : r)
1160 4478799 : i = save_state(r + 1)
1161 4478799 : w = save_state(r + 2)
1162 : end if
1163 :
1164 6673953 : If ((i .lt. 0) .or. (seed .ne. 0)) then ! Initialization necessary
1165 32 : If (seed .ne. 0) then ! v must be nonzero
1166 : v = seed
1167 : else
1168 0 : v = NOT(seed)
1169 : end if
1170 :
1171 2080 : do k = wlen, 1, -1 ! Avoid correlations for close seeds
1172 : ! This recurrence has period of 2^64-1
1173 2048 : v = IEOR(v, ISHFT(v, 7))
1174 2080 : v = IEOR(v, ISHFT(v, -9))
1175 : end do
1176 :
1177 : ! Initialize circular array
1178 32 : w = v
1179 2080 : do k = 0, r - 1
1180 : ! w = w + weyl
1181 2048 : if (w < 0_i8) then
1182 1020 : w = w + weyl
1183 1028 : else if ((huge(1_i8) - w) > weyl) then
1184 244 : w = w + weyl
1185 : else
1186 784 : tmp = -(huge(1_i8) - w - weyl)
1187 784 : w = tmp - huge(1_i8) - 2_i8
1188 : endif
1189 2048 : v = IEOR(v, ISHFT(v, 7))
1190 2048 : v = IEOR(v, ISHFT(v, -9))
1191 2080 : x(k) = v + w
1192 : end do
1193 :
1194 : ! Discard first 4*r results (Gimeno)
1195 32 : i = r - 1
1196 8224 : do k = 4 * r, 1, -1
1197 8192 : i = IAND(i + 1, r - 1)
1198 8192 : t = x(i)
1199 8192 : v = x(IAND(i + (r - s), r - 1))
1200 8192 : t = IEOR(t, ISHFT(t, a))
1201 8192 : t = IEOR(t, ISHFT(t, -b))
1202 8192 : v = IEOR(v, ISHFT(v, c))
1203 8192 : v = IEOR(v, IEOR(t, ISHFT(v, -d)))
1204 8224 : x(i) = v
1205 : end do
1206 : end if ! end of initialization
1207 :
1208 : ! Apart from initialization (above), this is the generator
1209 : v = 0_i8
1210 13347906 : Do While (v .eq. 0_i8)
1211 6673953 : i = IAND(i + 1, r - 1)
1212 6673953 : t = x(i)
1213 6673953 : v = x(IAND(i + (r - s), r - 1))
1214 6673953 : t = IEOR(t, ISHFT(t, a))
1215 6673953 : t = IEOR(t, ISHFT(t, -b))
1216 6673953 : v = IEOR(v, ISHFT(v, c))
1217 6673953 : v = IEOR(v, IEOR(t, ISHFT(v, -d)))
1218 6673953 : x(i) = v
1219 6673953 : w = w + weyl
1220 6673953 : v = v + w
1221 13347906 : v = ISHFT(v, -11)
1222 : End Do
1223 :
1224 6673953 : DoubleRealRN = t53 * v
1225 :
1226 6673953 : if(present(save_state)) then
1227 291123885 : save_state(1 : r) = x(0 : r - 1)
1228 4478829 : save_state(r + 1) = i
1229 4478829 : save_state(r + 2) = w
1230 300081543 : if ((r + 3) <= n_save_state) save_state(r + 3 : n_save_state) = 0
1231 : end if
1232 :
1233 0 : end subroutine xor4096d_0d
1234 :
1235 : ! -----------------------------------------------------------------------------
1236 :
1237 4 : subroutine xor4096d_1d(seed, DoubleRealRN, save_state)
1238 :
1239 : implicit none
1240 :
1241 : integer(i8), dimension(:), intent(in) :: seed
1242 : real(dp), dimension(size(seed, 1)), intent(out) :: DoubleRealRN
1243 : integer(i8), optional, dimension(size(seed, 1), n_save_state), intent(inout) :: save_state
1244 :
1245 : integer(i4) :: m
1246 : integer(i8) :: wlen, r, s, a, b, c, d
1247 : integer(i8) :: weyl = 7046029254386353131_i8
1248 : real(dp) :: t53 = 1.0_dp / 9007199254740992.0_dp ! = 0.5^53 = 1/2^53
1249 : integer(i8) :: k, j, tmp
1250 2 : integer(i8), dimension(size(seed)) :: t, v
1251 : integer(i8), dimension(:, :), allocatable, save :: x ! x(0) ... x(r-1)
1252 : integer(i8), dimension(:), allocatable, save :: w
1253 : integer(i8), dimension(:), allocatable, save :: i ! i<0 indicates first call
1254 :
1255 : !$omp threadprivate(x,i,w)
1256 :
1257 : ! produces a 53bit Integer Random Number (0...9 007 199 254 740 992) and
1258 : ! scales it afterwards to (0.0,1.0)
1259 :
1260 2 : wlen = 64_i8
1261 2 : r = 64_i8
1262 2 : s = 53_i8
1263 2 : a = 33_i8
1264 2 : b = 26_i8
1265 2 : c = 27_i8
1266 2 : d = 29_i8
1267 :
1268 2 : m = size(seed, 1)
1269 :
1270 6 : if (any(seed .eq. 0_i8) .and. any(seed .ne. 0_i8)) then
1271 0 : stop 'xor4096: seeds have to be eigther all 0 or all larger than 0'
1272 : end if
1273 :
1274 2 : if (present(save_state) .and. all(seed .eq. 0_i4)) then
1275 0 : if (allocated(x)) then
1276 0 : if (size(x, 1) .ne. m) then
1277 0 : deallocate(x)
1278 0 : deallocate(i)
1279 0 : deallocate(w)
1280 0 : allocate(i(m))
1281 0 : allocate(x(m, 0 : r - 1))
1282 0 : allocate(w(m))
1283 : end if
1284 : end if
1285 0 : if (.not. allocated(x)) then
1286 0 : allocate(i(m))
1287 0 : allocate(x(m, 0 : r - 1))
1288 0 : allocate(w(m))
1289 : end if
1290 0 : x(:, 0 : r - 1) = save_state(:, 1 : r)
1291 0 : i(:) = save_state(:, r + 1)
1292 0 : w(:) = save_state(:, r + 2)
1293 : end if
1294 :
1295 4 : if(all(seed .ne. 0_i4)) then
1296 2 : if (allocated(x)) then
1297 1 : deallocate(x)
1298 1 : deallocate(i)
1299 1 : deallocate(w)
1300 : end if
1301 :
1302 6 : allocate(i(m))
1303 6 : i = -1
1304 6 : allocate(x(m, 0 : r - 1))
1305 4 : allocate(w(m))
1306 : end if
1307 :
1308 4 : Do j = 1, m
1309 4 : If ((i(j) .lt. 0) .or. (seed(j) .ne. 0)) then ! Initialization necessary
1310 2 : If (seed(j) .ne. 0) then ! v must be nonzero
1311 2 : v(j) = seed(j)
1312 : else
1313 0 : v(j) = NOT(seed(j))
1314 : end if
1315 :
1316 130 : do k = wlen, 1, -1 ! Avoid correlations for close seeds
1317 : ! This recurrence has period of 2^64-1
1318 128 : v(j) = IEOR(v(j), ISHFT(v(j), 7))
1319 130 : v(j) = IEOR(v(j), ISHFT(v(j), -9))
1320 : end do
1321 :
1322 : ! Initialize circular array
1323 2 : w(j) = v(j)
1324 130 : do k = 0, r - 1
1325 : ! w(j) = w(j) + weyl
1326 128 : if (w(j) < 0_i8) then
1327 66 : w(j) = w(j) + weyl
1328 62 : else if ((huge(1_i8) - w(j)) > weyl) then
1329 14 : w(j) = w(j) + weyl
1330 : else
1331 48 : tmp = -(huge(1_i8) - w(j) - weyl)
1332 48 : w(j) = tmp - huge(1_i8) - 2_i8
1333 : endif
1334 128 : v(j) = IEOR(v(j), ISHFT(v(j), 7))
1335 128 : v(j) = IEOR(v(j), ISHFT(v(j), -9))
1336 130 : x(j, k) = v(j) + w(j)
1337 : end do
1338 :
1339 : ! Discard first 4*r results (Gimeno)
1340 2 : i(j) = r - 1
1341 514 : do k = 4 * r, 1, -1
1342 512 : i(j) = IAND(i(j) + 1, r - 1)
1343 512 : t(j) = x(j, i(j))
1344 512 : v(j) = x(j, IAND(i(j) + (r - s), r - 1))
1345 512 : t(j) = IEOR(t(j), ISHFT(t(j), a))
1346 512 : t(j) = IEOR(t(j), ISHFT(t(j), -b))
1347 512 : v(j) = IEOR(v(j), ISHFT(v(j), c))
1348 512 : v(j) = IEOR(v(j), IEOR(t(j), ISHFT(v(j), -d)))
1349 514 : x(j, i(j)) = v(j)
1350 : end do
1351 : end if ! end of initialization
1352 : end do
1353 :
1354 : ! Apart from initialization (above), this is the generator
1355 4 : v = 0_i8
1356 4 : Do j = 1, m
1357 6 : Do While (v(j) .eq. 0_i8)
1358 2 : i(j) = IAND(i(j) + 1, r - 1)
1359 2 : t(j) = x(j, i(j))
1360 2 : v(j) = x(j, IAND(i(j) + (r - s), r - 1))
1361 2 : t(j) = IEOR(t(j), ISHFT(t(j), a))
1362 2 : t(j) = IEOR(t(j), ISHFT(t(j), -b))
1363 2 : v(j) = IEOR(v(j), ISHFT(v(j), c))
1364 2 : v(j) = IEOR(v(j), IEOR(t(j), ISHFT(v(j), -d)))
1365 2 : x(j, i(j)) = v(j)
1366 2 : w(j) = w(j) + weyl
1367 2 : v(j) = v(j) + w(j)
1368 4 : v(j) = ISHFT(v(j), -11)
1369 : End Do
1370 : End Do
1371 :
1372 4 : DoubleRealRN = t53 * v
1373 :
1374 2 : if(present(save_state)) then
1375 258 : save_state(:, 1 : r) = x(:, 0 : r - 1)
1376 4 : save_state(:, r + 1) = i(:)
1377 4 : save_state(:, r + 2) = w(:)
1378 266 : if ((r + 3) <= n_save_state) save_state(:, r + 3 : n_save_state) = 0
1379 : end if
1380 :
1381 6673955 : end subroutine xor4096d_1d
1382 :
1383 : ! ------------------------------------------------------------------
1384 :
1385 0 : subroutine xor4096gf_0d(seed, SingleRealRN, save_state)
1386 :
1387 : implicit none
1388 :
1389 : integer(i4), intent(in) :: seed
1390 : real(SP), intent(out) :: SingleRealRN
1391 : integer(i4), optional, dimension(n_save_state), intent(inout) :: save_state
1392 :
1393 : integer(i4) :: wlen, r, s, a, b, c, d
1394 : integer(i4), save :: w
1395 : integer(i4), save :: x(0 : 127) ! x(0) ... x(r-1)
1396 : integer(i4) :: weyl = 1640531527_i4
1397 : integer(i4) :: t, v, tmp
1398 : integer(i4), save :: i = -1 ! i<0 indicates first call
1399 : integer(i4) :: k
1400 : real(SP) :: t24 = 1.0_SP / 16777216.0_SP ! = 0.5^24 = 1/2^24
1401 :
1402 0 : real(SP) :: rn1, rn2 ! uniform random numbers
1403 0 : real(SP) :: x1, x2, y1, ww ! for Box-Mueller transform
1404 : integer(i4), save :: Flag = 1 ! if Flag = 1 return y1 else return y2
1405 : real(SP), save :: y2
1406 :
1407 : !$omp threadprivate(x,i,w,y2,flag)
1408 :
1409 : ! produces a 24bit Integer Random Number (0...16777216) and
1410 : ! scales it afterwards to (0.0,1.0)
1411 : ! transform using polar-method
1412 :
1413 0 : wlen = 32
1414 0 : r = 128
1415 0 : s = 95
1416 0 : a = 17
1417 0 : b = 12
1418 0 : c = 13
1419 0 : d = 15
1420 :
1421 0 : if(present(save_state) .and. (seed .eq. 0)) then
1422 0 : x(0 : r - 1) = save_state(1 : r)
1423 0 : i = save_state(r + 1)
1424 0 : w = save_state(r + 2)
1425 0 : flag = save_state(r + 3)
1426 0 : y2 = transfer(save_state(r + 4), 1.0_sp)
1427 : end if
1428 :
1429 0 : If ((i .lt. 0) .or. (seed .ne. 0)) then ! Initialization necessary
1430 0 : If (seed .ne. 0) then ! v must be nonzero
1431 : v = seed
1432 : else
1433 0 : v = NOT(seed)
1434 : end if
1435 :
1436 0 : do k = wlen, 1, -1 ! Avoid correlations for close seeds
1437 : ! This recurrence has period of 2^32-1
1438 0 : v = IEOR(v, ISHFT(v, 13))
1439 0 : v = IEOR(v, ISHFT(v, -17))
1440 0 : v = IEOR(v, ISHFT(v, 5))
1441 : end do
1442 :
1443 : ! Initialize circular array
1444 0 : w = v
1445 0 : do k = 0, r - 1
1446 : ! w = w + weyl
1447 0 : if (w < 0_i4) then
1448 0 : w = w + weyl
1449 0 : else if ((huge(1_i4) - w) > weyl) then
1450 0 : w = w + weyl
1451 : else
1452 0 : tmp = -(huge(1_i4) - w - weyl)
1453 0 : w = tmp - huge(1_i4) - 2_i4
1454 : endif
1455 0 : v = IEOR(v, ISHFT(v, 13))
1456 0 : v = IEOR(v, ISHFT(v, -17))
1457 0 : v = IEOR(v, ISHFT(v, 5))
1458 0 : x(k) = v + w
1459 : end do
1460 :
1461 : ! Discard first 4*r results (Gimeno)
1462 0 : i = r - 1
1463 0 : do k = 4 * r, 1, -1
1464 0 : i = IAND(i + 1, r - 1)
1465 0 : t = x(i)
1466 0 : v = x(IAND(i + (r - s), r - 1))
1467 0 : t = IEOR(t, ISHFT(t, a))
1468 0 : t = IEOR(t, ISHFT(t, -b))
1469 0 : v = IEOR(v, ISHFT(v, c))
1470 0 : v = IEOR(v, IEOR(t, ISHFT(v, -d)))
1471 0 : x(i) = v
1472 : end do
1473 0 : Flag = 1
1474 : end if ! end of initialization
1475 :
1476 0 : If (Flag .eq. 1) then
1477 : ! Polar method of Box-Mueller-transform to generate Gaussian distributed random number
1478 : ww = 1.0_SP
1479 0 : do while (ww .ge. 1.0_SP)
1480 :
1481 : ! Apart from initialization (above), this is the generator
1482 : v = 0_i4
1483 0 : Do While (v .eq. 0_i4)
1484 0 : i = IAND(i + 1, r - 1)
1485 0 : t = x(i)
1486 0 : v = x(IAND(i + (r - s), r - 1))
1487 0 : t = IEOR(t, ISHFT(t, a))
1488 0 : t = IEOR(t, ISHFT(t, -b))
1489 0 : v = IEOR(v, ISHFT(v, c))
1490 0 : v = IEOR(v, IEOR(t, ISHFT(v, -d)))
1491 0 : x(i) = v
1492 0 : w = w + weyl
1493 0 : v = v + w
1494 0 : v = ISHFT(v, -8)
1495 : End Do
1496 :
1497 0 : rn1 = t24 * v
1498 :
1499 0 : v = 0_i4
1500 0 : Do While (v .eq. 0_i4)
1501 0 : i = IAND(i + 1, r - 1)
1502 0 : t = x(i)
1503 0 : v = x(IAND(i + (r - s), r - 1))
1504 0 : t = IEOR(t, ISHFT(t, a))
1505 0 : t = IEOR(t, ISHFT(t, -b))
1506 0 : v = IEOR(v, ISHFT(v, c))
1507 0 : v = IEOR(v, IEOR(t, ISHFT(v, -d)))
1508 0 : x(i) = v
1509 0 : w = w + weyl
1510 0 : v = v + w
1511 0 : v = ISHFT(v, -8)
1512 : End Do
1513 :
1514 0 : rn2 = t24 * v
1515 :
1516 0 : x1 = 2.0_SP * rn1 - 1.0_SP
1517 0 : x2 = 2.0_SP * rn2 - 1.0_SP
1518 :
1519 0 : ww = x1 * x1 + x2 * x2
1520 : end do ! end of polar method
1521 :
1522 0 : ww = Sqrt((-2.0_SP * Log(ww)) / ww)
1523 0 : y1 = x1 * ww
1524 0 : y2 = x2 * ww
1525 :
1526 : end if ! Only if Flag = 1
1527 :
1528 0 : If (Flag .eq. 1) then
1529 0 : Flag = 2
1530 0 : SingleRealRN = y1
1531 : else
1532 0 : Flag = 1
1533 0 : SingleRealRN = y2
1534 : end if
1535 :
1536 0 : if(present(save_state)) then
1537 0 : save_state(1 : r) = x(0 : r - 1)
1538 0 : save_state(r + 1) = i
1539 0 : save_state(r + 2) = w
1540 0 : save_state(r + 3) = flag
1541 0 : save_state(r + 4) = transfer(y2, 1_i4)
1542 : if ((r + 5) <= n_save_state) save_state(r + 5 : n_save_state) = 0
1543 : end if
1544 :
1545 2 : end subroutine xor4096gf_0d
1546 :
1547 : ! -----------------------------------------------------------------------------
1548 :
1549 0 : subroutine xor4096gf_1d(seed, SingleRealRN, save_state)
1550 :
1551 : implicit none
1552 :
1553 : integer(i4), dimension(:), intent(in) :: seed
1554 : real(SP), dimension(size(seed)), intent(out) :: SingleRealRN
1555 : integer(i4), optional, dimension(size(seed), n_save_state), intent(inout) :: save_state
1556 :
1557 : integer(i4) :: m
1558 : integer(i4) :: wlen, r, s, a, b, c, d
1559 : integer(i4) :: weyl = 1640531527_i4
1560 : integer(i4) :: k, j, tmp
1561 : real(SP) :: t24 = 1.0_SP / 16777216.0_SP ! = 0.5^24 = 1/2^24
1562 0 : integer(i4), dimension(size(seed)) :: t, v
1563 : integer(i4), dimension(:, :), allocatable, save :: x ! x(0) ... x(r-1)
1564 : integer(i4), dimension(:), allocatable, save :: i, w ! i<0 indicates first call
1565 :
1566 0 : real(SP), dimension(size(seed)) :: rn1, rn2 ! uniform random numbers
1567 0 : real(SP), dimension(size(seed)) :: x1, x2, y1, ww ! for Box-Mueller transform
1568 : real(SP), dimension(:), allocatable, save :: y2
1569 : integer(i4), dimension(:), allocatable, save :: Flag ! if Flag = 1 return y1 else return y2
1570 :
1571 : !$omp threadprivate(x,i,w,y2,flag)
1572 :
1573 : ! produces a 24bit Integer Random Number (0...16777216) and
1574 : ! scales it afterwards to (0.0,1.0)
1575 : ! transform using polar-method
1576 :
1577 0 : wlen = 32
1578 0 : r = 128
1579 0 : s = 95
1580 0 : a = 17
1581 0 : b = 12
1582 0 : c = 13
1583 0 : d = 15
1584 :
1585 0 : m = size(seed, 1)
1586 :
1587 0 : if (any(seed .eq. 0_i4) .and. any(seed .ne. 0_i4)) then
1588 0 : stop 'xor4096g: seeds have to be eigther all 0 or all larger than 0'
1589 : end if
1590 :
1591 0 : if (present(save_state) .and. all(seed .eq. 0_i4)) then
1592 0 : if (allocated(x)) then
1593 0 : if (size(x, 1) .ne. m) then
1594 0 : deallocate(x)
1595 0 : deallocate(i)
1596 0 : deallocate(w)
1597 0 : deallocate(flag)
1598 0 : deallocate(y2)
1599 0 : allocate(i(m))
1600 0 : allocate(x(m, 0 : r - 1))
1601 0 : allocate(w(m))
1602 0 : allocate(Flag(m))
1603 0 : allocate(y2(m))
1604 : end if
1605 : end if
1606 0 : if (.not. allocated(x)) then
1607 0 : allocate(i(m))
1608 0 : allocate(x(m, 0 : r - 1))
1609 0 : allocate(w(m))
1610 0 : allocate(Flag(m))
1611 0 : allocate(y2(m))
1612 : end if
1613 0 : x(:, 0 : r - 1) = save_state(:, 1 : r)
1614 0 : i(:) = save_state(:, r + 1)
1615 0 : w(:) = save_state(:, r + 2)
1616 0 : flag(:) = save_state(:, r + 3)
1617 0 : do j = 1, m
1618 0 : y2(j) = transfer(save_state(j, r + 4), 1.0_sp)
1619 : end do
1620 : end if
1621 :
1622 0 : if(all(seed .ne. 0_i4)) then
1623 0 : if (allocated(x)) then
1624 0 : deallocate(x)
1625 0 : deallocate(i)
1626 0 : deallocate(w)
1627 0 : deallocate(Flag)
1628 0 : deallocate(y2)
1629 : end if
1630 :
1631 0 : allocate(i(m))
1632 0 : i = -1
1633 0 : allocate(x(m, 0 : r - 1))
1634 0 : allocate(w(m))
1635 0 : allocate(Flag(m))
1636 0 : Flag = 1
1637 0 : allocate(y2(m))
1638 : end if
1639 :
1640 0 : Do j = 1, m !Loop over every stream
1641 0 : If ((i(j) .lt. 0) .or. (seed(j) .ne. 0)) then ! Initialization necessary
1642 0 : If (seed(j) .ne. 0) then ! v must be nonzero
1643 0 : v(j) = seed(j)
1644 : else
1645 0 : v(j) = NOT(seed(j))
1646 : end if
1647 :
1648 0 : do k = wlen, 1, -1 ! Avoid correlations for close seeds
1649 : ! This recurrence has period of 2^32-1
1650 0 : v(j) = IEOR(v(j), ISHFT(v(j), 13))
1651 0 : v(j) = IEOR(v(j), ISHFT(v(j), -17))
1652 0 : v(j) = IEOR(v(j), ISHFT(v(j), 5))
1653 : end do
1654 :
1655 : ! Initialize circular array
1656 0 : w(j) = v(j)
1657 0 : do k = 0, r - 1
1658 : ! w(j) = w(j) + weyl
1659 0 : if (w(j) < 0_i4) then
1660 0 : w(j) = w(j) + weyl
1661 0 : else if ((huge(1_i4) - w(j)) > weyl) then
1662 0 : w(j) = w(j) + weyl
1663 : else
1664 0 : tmp = -(huge(1_i4) - w(j) - weyl)
1665 0 : w(j) = tmp - huge(1_i4) - 2_i4
1666 : endif
1667 0 : v(j) = IEOR(v(j), ISHFT(v(j), 13))
1668 0 : v(j) = IEOR(v(j), ISHFT(v(j), -17))
1669 0 : v(j) = IEOR(v(j), ISHFT(v(j), 5))
1670 0 : x(j, k) = v(j) + w(j)
1671 : end do
1672 :
1673 : ! Discard first 4*r results (Gimeno)
1674 0 : i(j) = r - 1
1675 0 : do k = 4 * r, 1, -1
1676 0 : i(j) = IAND(i(j) + 1, r - 1)
1677 0 : t(j) = x(j, i(j))
1678 0 : v(j) = x(j, IAND(i(j) + (r - s), r - 1))
1679 0 : t(j) = IEOR(t(j), ISHFT(t(j), a))
1680 0 : t(j) = IEOR(t(j), ISHFT(t(j), -b))
1681 0 : v(j) = IEOR(v(j), ISHFT(v(j), c))
1682 0 : v(j) = IEOR(v(j), IEOR(t(j), ISHFT(v(j), -d)))
1683 0 : x(j, i(j)) = v(j)
1684 : end do
1685 0 : Flag(j) = 1
1686 : end if ! end of initialization
1687 : end do
1688 :
1689 0 : Do j = 1, m !Loop over every stream
1690 0 : If (Flag(j) .eq. 1) then
1691 : ! Polar method of Box-Mueller-transform to generate Gaussian distributed random number
1692 0 : ww(j) = 1.0_SP
1693 0 : do while (ww(j) .ge. 1.0_SP)
1694 :
1695 : ! Apart from initialization (above), this is the generator
1696 0 : v(j) = 0_i4
1697 0 : Do While (v(j) .eq. 0_i4)
1698 0 : i(j) = IAND(i(j) + 1, r - 1)
1699 0 : t(j) = x(j, i(j))
1700 0 : v(j) = x(j, IAND(i(j) + (r - s), r - 1))
1701 0 : t(j) = IEOR(t(j), ISHFT(t(j), a))
1702 0 : t(j) = IEOR(t(j), ISHFT(t(j), -b))
1703 0 : v(j) = IEOR(v(j), ISHFT(v(j), c))
1704 0 : v(j) = IEOR(v(j), IEOR(t(j), ISHFT(v(j), -d)))
1705 0 : x(j, i(j)) = v(j)
1706 0 : w(j) = w(j) + weyl
1707 0 : v(j) = v(j) + w(j)
1708 0 : v(j) = ISHFT(v(j), -8)
1709 : End Do
1710 :
1711 0 : rn1(j) = t24 * v(j)
1712 :
1713 0 : v(j) = 0_i4
1714 0 : Do While (v(j) .eq. 0_i4)
1715 0 : i(j) = IAND(i(j) + 1, r - 1)
1716 0 : t(j) = x(j, i(j))
1717 0 : v(j) = x(j, IAND(i(j) + (r - s), r - 1))
1718 0 : t(j) = IEOR(t(j), ISHFT(t(j), a))
1719 0 : t(j) = IEOR(t(j), ISHFT(t(j), -b))
1720 0 : v(j) = IEOR(v(j), ISHFT(v(j), c))
1721 0 : v(j) = IEOR(v(j), IEOR(t(j), ISHFT(v(j), -d)))
1722 0 : x(j, i(j)) = v(j)
1723 0 : w(j) = w(j) + weyl
1724 0 : v(j) = v(j) + w(j)
1725 0 : v(j) = ISHFT(v(j), -8)
1726 : End Do
1727 :
1728 0 : rn2(j) = t24 * v(j)
1729 :
1730 0 : x1(j) = 2.0_SP * rn1(j) - 1.0_SP
1731 0 : x2(j) = 2.0_SP * rn2(j) - 1.0_SP
1732 0 : ww(j) = x1(j) * x1(j) + x2(j) * x2(j)
1733 : end do ! end of polar method
1734 :
1735 0 : ww(j) = Sqrt((-2.0_SP * Log(ww(j))) / ww(j))
1736 0 : y1(j) = x1(j) * ww(j)
1737 0 : y2(j) = x2(j) * ww(j)
1738 : end if ! Only if Flag = 1
1739 : end do ! Loop over each stream
1740 :
1741 0 : Do j = 1, m
1742 0 : If (Flag(j) .eq. 1) then
1743 0 : Flag(j) = 2
1744 0 : SingleRealRN(j) = y1(j)
1745 : else
1746 0 : Flag(j) = 1
1747 0 : SingleRealRN(j) = y2(j)
1748 : end if
1749 : end Do
1750 :
1751 0 : if(present(save_state)) then
1752 0 : save_state(:, 1 : r) = x(:, 0 : r - 1)
1753 0 : save_state(:, r + 1) = i(:)
1754 0 : save_state(:, r + 2) = w(:)
1755 0 : save_state(:, r + 3) = flag(:)
1756 0 : do j = 1, m
1757 0 : save_state(j, r + 4) = transfer(y2(j), 1_i4)
1758 : end do
1759 : if ((r + 5) <= n_save_state) save_state(:, r + 5 : n_save_state) = 0
1760 : end if
1761 :
1762 0 : end subroutine xor4096gf_1d
1763 :
1764 : ! -----------------------------------------------------------------------------
1765 :
1766 553488 : subroutine xor4096gd_0d(seed, DoubleRealRN, save_state)
1767 :
1768 : implicit none
1769 :
1770 : integer(i8), intent(in) :: seed
1771 : real(DP), intent(out) :: DoubleRealRN
1772 : integer(i8), optional, dimension(n_save_state), intent(inout) :: save_state
1773 :
1774 : integer(i8) :: wlen, r, s, a, b, c, d
1775 :
1776 : integer(i8), save :: w
1777 : integer(i8), save :: x(0 : 63) ! x(0) ... x(r-1)
1778 : integer(i8) :: weyl = 7046029254386353131_i8
1779 : integer(i8) :: t, v, tmp
1780 : integer(i8), save :: i = -1_i8 ! i<0 indicates first call
1781 : integer(i8) :: k
1782 :
1783 : real(DP) :: t53 = 1.0_DP / 9007199254740992.0_DP ! = 0.5^53 = 1/2^53
1784 :
1785 553488 : real(DP) :: rn1, rn2 ! uniform random numbers
1786 553488 : real(DP) :: x1, x2, y1, ww ! for Box-Mueller transform
1787 : real(DP), save :: y2
1788 : integer(i8), save :: Flag = 1_i8 ! if Flag = 1 return y1 else return y2
1789 :
1790 : !$omp threadprivate(x,i,w,y2,flag)
1791 :
1792 : ! produces a 53bit Integer Random Number (0...9 007 199 254 740 992) and
1793 : ! scales it afterwards to (0.0,1.0)
1794 : ! transform using polar-method
1795 :
1796 553488 : wlen = 64_i8
1797 553488 : r = 64_i8
1798 553488 : s = 53_i8
1799 553488 : a = 33_i8
1800 553488 : b = 26_i8
1801 553488 : c = 27_i8
1802 553488 : d = 29_i8
1803 :
1804 553488 : if(present(save_state) .and. (seed .eq. 0)) then
1805 17691505 : x(0 : r - 1) = save_state(1 : r)
1806 272177 : i = save_state(r + 1)
1807 272177 : w = save_state(r + 2)
1808 272177 : flag = save_state(r + 3)
1809 272177 : y2 = transfer(save_state(r + 4), 1.0_dp)
1810 : end if
1811 :
1812 553488 : If ((i .lt. 0_i8) .or. (seed .ne. 0_i8)) then ! Initialization necessary
1813 16 : If (seed .ne. 0) then ! v must be nonzero
1814 : v = seed
1815 : else
1816 0 : v = NOT(seed)
1817 : end if
1818 :
1819 1040 : do k = wlen, 1, -1 ! Avoid correlations for close seeds
1820 : ! This recurrence has period of 2^64-1
1821 1024 : v = IEOR(v, ISHFT(v, 7))
1822 1040 : v = IEOR(v, ISHFT(v, -9))
1823 : end do
1824 :
1825 : ! Initialize circular array
1826 16 : w = v
1827 1040 : do k = 0, r - 1
1828 : ! w = w + weyl
1829 1024 : if (w < 0_i8) then
1830 513 : w = w + weyl
1831 511 : else if ((huge(1_i8) - w) > weyl) then
1832 122 : w = w + weyl
1833 : else
1834 389 : tmp = -(huge(1_i8) - w - weyl)
1835 389 : w = tmp - huge(1_i8) - 2_i8
1836 : endif
1837 1024 : v = IEOR(v, ISHFT(v, 7))
1838 1024 : v = IEOR(v, ISHFT(v, -9))
1839 1040 : x(k) = v + w
1840 : end do
1841 :
1842 : ! Discard first 4*r results (Gimeno)
1843 16 : i = r - 1
1844 4112 : do k = 4 * r, 1, -1
1845 4096 : i = IAND(i + 1, r - 1)
1846 4096 : t = x(i)
1847 4096 : v = x(IAND(i + (r - s), r - 1))
1848 4096 : t = IEOR(t, ISHFT(t, a))
1849 4096 : t = IEOR(t, ISHFT(t, -b))
1850 4096 : v = IEOR(v, ISHFT(v, c))
1851 4096 : v = IEOR(v, IEOR(t, ISHFT(v, -d)))
1852 4112 : x(i) = v
1853 : end do
1854 16 : Flag = 1_i8
1855 : end if ! end of initialization
1856 :
1857 553488 : If (Flag .eq. 1_i8) then
1858 : ! Polar method of Box-Mueller-transform to generate Gaussian distributed random number
1859 : ww = 1.0_DP
1860 628925 : do while (ww .ge. 1.0_DP)
1861 :
1862 : ! Apart from initialization (above), this is the generator
1863 : v = 0_i8
1864 704352 : Do While (v .eq. 0_i8)
1865 352176 : i = IAND(i + 1, r - 1)
1866 352176 : t = x(i)
1867 352176 : v = x(IAND(i + (r - s), r - 1))
1868 352176 : t = IEOR(t, ISHFT(t, a))
1869 352176 : t = IEOR(t, ISHFT(t, -b))
1870 352176 : v = IEOR(v, ISHFT(v, c))
1871 352176 : v = IEOR(v, IEOR(t, ISHFT(v, -d)))
1872 352176 : x(i) = v
1873 352176 : w = w + weyl
1874 352176 : v = v + w
1875 704352 : v = ISHFT(v, -11)
1876 : End Do
1877 :
1878 352176 : rn1 = t53 * v
1879 :
1880 352176 : v = 0_i8
1881 704352 : Do While (v .eq. 0_i8)
1882 352176 : i = IAND(i + 1, r - 1)
1883 352176 : t = x(i)
1884 352176 : v = x(IAND(i + (r - s), r - 1))
1885 352176 : t = IEOR(t, ISHFT(t, a))
1886 352176 : t = IEOR(t, ISHFT(t, -b))
1887 352176 : v = IEOR(v, ISHFT(v, c))
1888 352176 : v = IEOR(v, IEOR(t, ISHFT(v, -d)))
1889 352176 : x(i) = v
1890 352176 : w = w + weyl
1891 352176 : v = v + w
1892 704352 : v = ISHFT(v, -11)
1893 : End Do
1894 :
1895 352176 : rn2 = t53 * v
1896 :
1897 352176 : x1 = 2.0_DP * rn1 - 1.0_DP
1898 352176 : x2 = 2.0_DP * rn2 - 1.0_DP
1899 352176 : ww = x1 * x1 + x2 * x2
1900 : end do ! end of polar method
1901 :
1902 276749 : ww = Sqrt((-2.0_DP * Log(ww)) / ww)
1903 276749 : y1 = x1 * ww
1904 276749 : y2 = x2 * ww
1905 :
1906 : end if ! Only if Flag = 1
1907 :
1908 553488 : If (Flag .eq. 1) then
1909 276749 : Flag = 2
1910 276749 : DoubleRealRN = y1
1911 : else
1912 276739 : Flag = 1
1913 276739 : DoubleRealRN = y2
1914 : end if
1915 :
1916 553488 : if(present(save_state)) then
1917 17692415 : save_state(1 : r) = x(0 : r - 1)
1918 272191 : save_state(r + 1) = i
1919 272191 : save_state(r + 2) = w
1920 272191 : save_state(r + 3) = flag
1921 272191 : save_state(r + 4) = transfer(y2, 1_i8)
1922 17692415 : if ((r + 5) <= n_save_state) save_state(r + 5 : n_save_state) = 0
1923 : end if
1924 :
1925 0 : end subroutine xor4096gd_0d
1926 :
1927 : ! -----------------------------------------------------------------------------
1928 :
1929 4 : subroutine xor4096gd_1d(seed, DoubleRealRN, save_state)
1930 :
1931 : implicit none
1932 :
1933 : integer(i8), dimension(:), intent(in) :: seed
1934 : real(DP), dimension(size(seed)), intent(out) :: DoubleRealRN
1935 : integer(i8), optional, dimension(size(seed), n_save_state), intent(inout) :: save_state
1936 :
1937 : integer(i4) :: m
1938 : integer(i8) :: wlen, r, s, a, b, c, d
1939 : integer(i8) :: weyl = 7046029254386353131_i8 !Z'61C88647' = Hexadecimal notation
1940 : integer(i8) :: k, j, tmp
1941 : real(DP) :: t53 = 1.0_DP / 9007199254740992.0_DP ! = 0.5^24 = 1/2^24
1942 2 : integer(i8), dimension(size(seed)) :: t, v
1943 : integer(i8), dimension(:, :), allocatable, save :: x ! x(0) ... x(r-1)
1944 : integer(i8), dimension(:), allocatable, save :: i, w ! i<0 indicates first call
1945 :
1946 8 : real(DP), dimension(size(seed)) :: rn1, rn2 ! uniform random numbers
1947 10 : real(DP), dimension(size(seed)) :: x1, x2, y1, ww ! for Box-Mueller transform
1948 : real(DP), dimension(:), allocatable, save :: y2
1949 : integer(i8), dimension(:), allocatable, save :: Flag ! if Flag = 1 return y1 else return y2
1950 :
1951 : !$omp threadprivate(x,i,w,y2,flag)
1952 :
1953 : ! produces a 53bit Integer Random Number (0...9 007 199 254 740 992) and
1954 : ! scales it afterwards to (0.0,1.0)
1955 : ! transform using polar-method
1956 :
1957 2 : wlen = 64_i8
1958 2 : r = 64_i8
1959 2 : s = 53_i8
1960 2 : a = 33_i8
1961 2 : b = 26_i8
1962 2 : c = 27_i8
1963 2 : d = 29_i8
1964 :
1965 2 : m = size(seed, 1)
1966 :
1967 6 : if (any(seed .eq. 0_i8) .and. any(seed .ne. 0_i8)) then
1968 0 : stop 'xor4096g: seeds have to be eigther all 0 or all larger than 0'
1969 : end if
1970 :
1971 2 : if (present(save_state) .and. all(seed .eq. 0_i8)) then
1972 0 : if (allocated(x)) then
1973 0 : if (size(x, 1) .ne. m) then
1974 0 : deallocate(x)
1975 0 : deallocate(i)
1976 0 : deallocate(w)
1977 0 : deallocate(flag)
1978 0 : deallocate(y2)
1979 0 : allocate(i(m))
1980 0 : allocate(x(m, 0 : r - 1))
1981 0 : allocate(w(m))
1982 0 : allocate(Flag(m))
1983 0 : allocate(y2(m))
1984 : end if
1985 : end if
1986 0 : if (.not. allocated(x)) then
1987 0 : allocate(i(m))
1988 0 : allocate(x(m, 0 : r - 1))
1989 0 : allocate(w(m))
1990 0 : allocate(Flag(m))
1991 0 : allocate(y2(m))
1992 : end if
1993 0 : x(:, 0 : r - 1) = save_state(:, 1 : r)
1994 0 : i(:) = save_state(:, r + 1)
1995 0 : w(:) = save_state(:, r + 2)
1996 0 : flag(:) = save_state(:, r + 3)
1997 0 : do j = 1, m
1998 0 : y2(j) = transfer(save_state(j, r + 4), 1.0_dp)
1999 : end do
2000 : end if
2001 :
2002 4 : if(all(seed .ne. 0_i8)) then
2003 2 : if (allocated(x)) then
2004 1 : deallocate(x)
2005 1 : deallocate(i)
2006 1 : deallocate(w)
2007 1 : deallocate(Flag)
2008 1 : deallocate(y2)
2009 : end if
2010 :
2011 6 : allocate(i(m))
2012 4 : i = -1
2013 6 : allocate(x(m, 0 : r - 1))
2014 4 : allocate(w(m))
2015 4 : allocate(Flag(m))
2016 4 : Flag = 1
2017 4 : allocate(y2(m))
2018 : end if
2019 :
2020 4 : Do j = 1, m !Loop over every stream
2021 4 : If ((i(j) .lt. 0) .or. (seed(j) .ne. 0)) then ! Initialization necessary
2022 2 : If (seed(j) .ne. 0) then ! v must be nonzero
2023 2 : v(j) = seed(j)
2024 : else
2025 0 : v(j) = NOT(seed(j))
2026 : end if
2027 :
2028 130 : do k = wlen, 1, -1 ! Avoid correlations for close seeds
2029 : ! This recurrence has period of 2^64-1
2030 128 : v(j) = IEOR(v(j), ISHFT(v(j), 7))
2031 130 : v(j) = IEOR(v(j), ISHFT(v(j), -9))
2032 : end do
2033 :
2034 : ! Initialize circular array
2035 2 : w(j) = v(j)
2036 130 : do k = 0, r - 1
2037 : ! w(j) = w(j) + weyl
2038 128 : if (w(j) < 0_i8) then
2039 62 : w(j) = w(j) + weyl
2040 66 : else if ((huge(1_i8) - w(j)) > weyl) then
2041 16 : w(j) = w(j) + weyl
2042 : else
2043 50 : tmp = -(huge(1_i8) - w(j) - weyl)
2044 50 : w(j) = tmp - huge(1_i8) - 2_i8
2045 : endif
2046 128 : v(j) = IEOR(v(j), ISHFT(v(j), 7))
2047 128 : v(j) = IEOR(v(j), ISHFT(v(j), -9))
2048 130 : x(j, k) = v(j) + w(j)
2049 : end do
2050 :
2051 : ! Discard first 4*r results (Gimeno)
2052 2 : i(j) = r - 1
2053 514 : do k = 4 * r, 1, -1
2054 512 : i(j) = IAND(i(j) + 1, r - 1)
2055 512 : t(j) = x(j, i(j))
2056 512 : v(j) = x(j, IAND(i(j) + (r - s), r - 1))
2057 512 : t(j) = IEOR(t(j), ISHFT(t(j), a))
2058 512 : t(j) = IEOR(t(j), ISHFT(t(j), -b))
2059 512 : v(j) = IEOR(v(j), ISHFT(v(j), c))
2060 512 : v(j) = IEOR(v(j), IEOR(t(j), ISHFT(v(j), -d)))
2061 514 : x(j, i(j)) = v(j)
2062 : end do
2063 2 : Flag(j) = 1
2064 : end if ! end of initialization
2065 : end do
2066 :
2067 4 : Do j = 1, m !Loop over every stream
2068 4 : If (Flag(j) .eq. 1) then
2069 : ! Polar method of Box-Mueller-transform to generate Gaussian distributed random number
2070 2 : ww(j) = 1.0_DP
2071 4 : do while (ww(j) .ge. 1.0_DP)
2072 :
2073 : ! Apart from initialization (above), this is the generator
2074 2 : v(j) = 0_i8
2075 4 : Do While (v(j) .eq. 0_i8)
2076 2 : i(j) = IAND(i(j) + 1, r - 1)
2077 2 : t(j) = x(j, i(j))
2078 2 : v(j) = x(j, IAND(i(j) + (r - s), r - 1))
2079 2 : t(j) = IEOR(t(j), ISHFT(t(j), a))
2080 2 : t(j) = IEOR(t(j), ISHFT(t(j), -b))
2081 2 : v(j) = IEOR(v(j), ISHFT(v(j), c))
2082 2 : v(j) = IEOR(v(j), IEOR(t(j), ISHFT(v(j), -d)))
2083 2 : x(j, i(j)) = v(j)
2084 2 : w(j) = w(j) + weyl
2085 2 : v(j) = v(j) + w(j)
2086 4 : v(j) = ISHFT(v(j), -11)
2087 : End Do
2088 :
2089 2 : rn1(j) = t53 * v(j)
2090 :
2091 2 : v(j) = 0_i8
2092 4 : Do While (v(j) .eq. 0_i8)
2093 2 : i(j) = IAND(i(j) + 1, r - 1)
2094 2 : t(j) = x(j, i(j))
2095 2 : v(j) = x(j, IAND(i(j) + (r - s), r - 1))
2096 2 : t(j) = IEOR(t(j), ISHFT(t(j), a))
2097 2 : t(j) = IEOR(t(j), ISHFT(t(j), -b))
2098 2 : v(j) = IEOR(v(j), ISHFT(v(j), c))
2099 2 : v(j) = IEOR(v(j), IEOR(t(j), ISHFT(v(j), -d)))
2100 2 : x(j, i(j)) = v(j)
2101 2 : w(j) = w(j) + weyl
2102 2 : v(j) = v(j) + w(j)
2103 4 : v(j) = ISHFT(v(j), -11)
2104 : End Do
2105 :
2106 2 : rn2(j) = t53 * v(j)
2107 :
2108 2 : x1(j) = 2.0_DP * rn1(j) - 1.0_DP
2109 2 : x2(j) = 2.0_DP * rn2(j) - 1.0_DP
2110 2 : ww(j) = x1(j) * x1(j) + x2(j) * x2(j)
2111 : end do ! end of polar method
2112 :
2113 2 : ww(j) = Sqrt((-2.0_DP * Log(ww(j))) / ww(j))
2114 2 : y1(j) = x1(j) * ww(j)
2115 2 : y2(j) = x2(j) * ww(j)
2116 :
2117 : end if ! Only if Flag = 1
2118 : end do ! Loop over each stream
2119 :
2120 4 : Do j = 1, m
2121 4 : If (Flag(j) .eq. 1) then
2122 2 : Flag(j) = 2
2123 2 : DoubleRealRN(j) = y1(j)
2124 : else
2125 0 : Flag(j) = 1
2126 0 : DoubleRealRN(j) = y2(j)
2127 : end if
2128 : End do
2129 :
2130 2 : if(present(save_state)) then
2131 258 : save_state(:, 1 : r) = x(:, 0 : r - 1)
2132 4 : save_state(:, r + 1) = i(:)
2133 4 : save_state(:, r + 2) = w(:)
2134 4 : save_state(:, r + 3) = flag(:)
2135 4 : do j = 1, m
2136 4 : save_state(j, r + 4) = transfer(y2(j), 1_i8)
2137 : end do
2138 258 : if ((r + 5) <= n_save_state) save_state(:, r + 5 : n_save_state) = 0
2139 : end if
2140 :
2141 553490 : end subroutine xor4096gd_1d
2142 :
2143 : end module mo_xor4096
|
