11 integer,
parameter :: dp = kind(0.0d0)
14 integer,
parameter :: i4 = selected_int_kind(9)
17 integer,
parameter :: i8 = selected_int_kind(18)
22 integer(i8),
private :: s(2) = [123456789_i8, 987654321_i8]
24 procedure, non_overridable :: set_seed
25 procedure, non_overridable :: jump
26 procedure, non_overridable :: int_4
27 procedure, non_overridable :: int_8
28 procedure, non_overridable :: unif_01
29 procedure, non_overridable :: unif_01_vec
30 procedure, non_overridable :: normal
31 procedure, non_overridable :: two_normals
32 procedure, non_overridable :: poisson
33 procedure, non_overridable :: circle
34 procedure, non_overridable :: sphere
35 procedure, non_overridable :: next
40 type(
rng_t),
allocatable :: rngs(:)
42 procedure, non_overridable :: init_parallel
51#if defined(__NVCOMPILER) || (defined(USE_INTRINSIC_SHIFT) && USE_INTRINSIC_SHIFT==0)
54 pure function shiftl(val, shift)
result(res_value)
55 integer(i8),
intent(in) :: val
56 integer,
intent(in) :: shift
57 integer(i8) :: res_value
58 integer(i8) :: bit_mask1, bit_mask2
62 bit_mask2=-bit_mask1-1
66 res_value = ior(lshift(iand(val, bit_mask1), shift),bit_mask2)
68 res_value = lshift(val, shift)
73 pure function shiftr(val, shift)
result(res_value)
74 integer(i8),
intent(in) :: val
75 integer,
intent(in) :: shift
76 integer(i8) :: res_value
77 integer(i8) :: bit_mask1, bit_mask2
81 bit_mask2=-bit_mask1-1
84 res_value = ior(rshift(iand(val, bit_mask1), shift),bit_mask2)
86 res_value = rshift(val, shift)
94 subroutine init_parallel(self, n_proc, rng)
95 class(
prng_t),
intent(inout) :: self
96 type(
rng_t),
intent(inout) :: rng
97 integer,
intent(in) :: n_proc
100 allocate(self%rngs(n_proc))
104 self%rngs(n) = self%rngs(n-1)
105 call self%rngs(n)%jump()
107 end subroutine init_parallel
110 subroutine set_seed(self, the_seed)
111 class(
rng_t),
intent(inout) :: self
112 integer(i8),
intent(in) :: the_seed(2)
118 end subroutine set_seed
123 subroutine jump(self)
124 class(
rng_t),
intent(inout) :: self
126 integer(i8) :: global_time(2), dummy
129 integer(i8),
parameter :: jmp_c(2) = &
130 (/-4707382666127344949_i8, -2852180941702784734_i8/)
135 if (iand(jmp_c(i),
shiftl(1_i8, b)) /= 0)
then
136 global_time = ieor(global_time, self%s)
146 integer(i4) function int_4(self)
147 class(
rng_t),
intent(inout) :: self
148 int_4 = int(self%next(), i4)
152 integer(i8) function int_8(self)
153 class(
rng_t),
intent(inout) :: self
158 real(dp) function unif_01(self)
159 class(
rng_t),
intent(inout) :: self
164 x = ior(
shiftl(1023_i8, 52), shiftr(x, 12))
165 unif_01 = transfer(x, tmp) - 1.0_dp
169 subroutine unif_01_vec(self, rr)
170 class(
rng_t),
intent(inout) :: self
171 real(dp),
intent(out) :: rr(:)
175 rr(i) = self%unif_01()
177 end subroutine unif_01_vec
181 function two_normals(self)
result(rands)
182 class(
rng_t),
intent(inout) :: self
183 real(dp) :: rands(2), sum_sq
186 rands(1) = 2 * self%unif_01() - 1
187 rands(2) = 2 * self%unif_01() - 1
188 sum_sq = sum(rands**2)
189 if (sum_sq < 1.0_dp .and. sum_sq > 0.0_dp)
exit
191 rands = rands * sqrt(-2 * log(sum_sq) / sum_sq)
192 end function two_normals
195 real(dp) function normal(self)
196 class(
rng_t),
intent(inout) :: self
199 rands = self%two_normals()
205 function poisson(self, lambda)
result(rr)
206 class(
rng_t),
intent(inout) :: self
207 real(dp),
intent(in) :: lambda
217 p = p * self%unif_01()
222 function circle(self, radius)
result(xy)
223 class(
rng_t),
intent(inout) :: self
224 real(dp),
intent(in) :: radius
225 real(dp) :: rands(2), xy(2)
230 rands(1) = 2 * self%unif_01() - 1
231 rands(2) = 2 * self%unif_01() - 1
232 sum_sq = sum(rands**2)
233 if (sum_sq <= 1)
exit
236 xy(1) = (rands(1)**2 - rands(2)**2) / sum_sq
237 xy(2) = 2 * rands(1) * rands(2) / sum_sq
242 function sphere(self, radius)
result(xyz)
243 class(
rng_t),
intent(inout) :: self
244 real(dp),
intent(in) :: radius
245 real(dp) :: rands(2), xyz(3)
246 real(dp) :: sum_sq, tmp_sqrt
250 rands(1) = 2 * self%unif_01() - 1
251 rands(2) = 2 * self%unif_01() - 1
252 sum_sq = sum(rands**2)
253 if (sum_sq <= 1)
exit
256 tmp_sqrt = sqrt(1 - sum_sq)
257 xyz(1:2) = 2 * rands(1:2) * tmp_sqrt
258 xyz(3) = 1 - 2 * sum_sq
263 function next(self)
result(res)
264 class(
rng_t),
intent(inout) :: self
266 integer(i8) :: global_time(2)
269 res = global_time(1) + global_time(2)
270 global_time(2) = ieor(global_time(1), global_time(2))
271 self%s(1) = ieor(ieor(rotl(global_time(1), 55), global_time(2)),
shiftl(global_time(2), 14))
272 self%s(2) = rotl(global_time(2), 36)
276 pure function rotl(x, k)
result(res)
277 integer(i8),
intent(in) :: x
278 integer,
intent(in) :: k
281 res = ior(
shiftl(x, k), shiftr(x, 64 - k))
Module for pseudo random number generation. The internal pseudo random generator is the xoroshiro128p...
pure integer(i8) function shiftl(val, shift)
added for nvidia compilers
Parallel random number generator type.
Random number generator type, which contains the state.