|
| 1 | +#:include "common.fypp" |
| 2 | +module stdlib_stats_distribution_PRNG |
| 3 | + use stdlib_kinds, only: int8, int16, int32, int64 |
| 4 | + use stdlib_error, only: error_stop |
| 5 | + implicit none |
| 6 | + private |
| 7 | + integer, parameter :: MAX_INT_BIT_SIZE = bit_size(1_int64) |
| 8 | + integer(int64) :: st(4) ! internal states for xoshiro256ss function |
| 9 | + integer(int64) :: si = 614872703977525537_int64 ! default seed value |
| 10 | + logical :: seed_initialized = .false. |
| 11 | + |
| 12 | + public :: random_seed |
| 13 | + public :: dist_rand |
| 14 | + |
| 15 | + |
| 16 | + interface dist_rand |
| 17 | + !! Version experimental |
| 18 | + !! |
| 19 | + !! Generation of random integers with different kinds |
| 20 | + !! ([Specification](../page/specs/stdlib_stats_distribution_PRNG.html# |
| 21 | + !! description)) |
| 22 | + #:for k1, t1 in INT_KINDS_TYPES |
| 23 | + module procedure dist_rand_${t1[0]}$${k1}$ |
| 24 | + #:endfor |
| 25 | + end interface dist_rand |
| 26 | + |
| 27 | + interface random_seed |
| 28 | + !! Version experimental |
| 29 | + !! |
| 30 | + !! Set seed value for random number generator |
| 31 | + !! ([Specification](../page/specs/stdlib_stats_distribution_PRNG.html# |
| 32 | + !! description)) |
| 33 | + !! |
| 34 | + #:for k1, t1 in INT_KINDS_TYPES |
| 35 | + module procedure random_distribution_seed_${t1[0]}$${k1}$ |
| 36 | + #:endfor |
| 37 | + end interface random_seed |
| 38 | + |
| 39 | + |
| 40 | + contains |
| 41 | + |
| 42 | + #:for k1, t1 in INT_KINDS_TYPES |
| 43 | + function dist_rand_${t1[0]}$${k1}$(n) result(res) |
| 44 | + !! Random integer generation for various kinds |
| 45 | + !! result = [-2^k, 2^k - 1], k = 7, 15, 31, 63, depending on input kind |
| 46 | + !! Result will be operated by bitwise operators to generate desired integer |
| 47 | + !! and real pseudorandom numbers |
| 48 | + !! |
| 49 | + ${t1}$, intent(in) :: n |
| 50 | + ${t1}$ :: res |
| 51 | + integer :: k |
| 52 | + |
| 53 | + k = MAX_INT_BIT_SIZE - bit_size(n) |
| 54 | + if(k < 0) call error_stop("Error(dist_rand): Integer bit size is" & |
| 55 | + //" greater than 64bit") |
| 56 | + res = shiftr(xoshiro256ss( ), k) |
| 57 | + end function dist_rand_${t1[0]}$${k1}$ |
| 58 | + |
| 59 | + #:endfor |
| 60 | + |
| 61 | + function xoshiro256ss( ) result (res) |
| 62 | + ! Generate random 64-bit integers |
| 63 | + ! Written in 2018 by David Blackman and Sebastiano Vigna ( [email protected]) |
| 64 | + ! http://prng.di.unimi.it/xoshiro256starstar.c |
| 65 | + ! |
| 66 | + ! This is xoshiro256** 1.0, one of our all-purpose, rock-solid |
| 67 | + ! generators. It has excellent (sub-ns) speed, a state (256 bits) that is |
| 68 | + ! large enough for any parallel application, and it passes all tests we |
| 69 | + ! are aware of. |
| 70 | + ! |
| 71 | + ! The state must be seeded so that it is not everywhere zero. If you have |
| 72 | + ! a 64-bit seed, we suggest to seed a splitmix64 generator and use its |
| 73 | + ! output to fill st. |
| 74 | + ! |
| 75 | + ! Fortran 90 version translated from C by Jim-215-Fisher |
| 76 | + ! |
| 77 | + integer(int64) :: res, t |
| 78 | + |
| 79 | + if(.not. seed_initialized) call random_distribution_seed_iint64(si,t) |
| 80 | + res = rol64(st(2) * 5, 7) * 9 |
| 81 | + t = shiftl(st(2), 17) |
| 82 | + st(3) = ieor(st(3), st(1)) |
| 83 | + st(4) = ieor(st(4), st(2)) |
| 84 | + st(2) = ieor(st(2), st(3)) |
| 85 | + st(1) = ieor(st(1), st(4)) |
| 86 | + st(3) = ieor(st(3), t) |
| 87 | + st(4) = rol64(st(4), 45) |
| 88 | + end function xoshiro256ss |
| 89 | + |
| 90 | + pure function rol64(x, k) result(res) |
| 91 | + integer(int64), intent(in) :: x |
| 92 | + integer, intent(in) :: k |
| 93 | + integer(int64) :: t1, t2, res |
| 94 | + |
| 95 | + t1 = shiftr(x, (64 - k)) |
| 96 | + t2 = shiftl(x, k) |
| 97 | + res = ior(t1, t2) |
| 98 | + end function rol64 |
| 99 | + |
| 100 | + |
| 101 | + function splitmix64(s) result(res) |
| 102 | + ! Written in 2015 by Sebastiano Vigna ( [email protected]) |
| 103 | + ! This is a fixed-increment version of Java 8's SplittableRandom |
| 104 | + ! generator. |
| 105 | + ! See http://dx.doi.org/10.1145/2714064.2660195 and |
| 106 | + ! http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html |
| 107 | + ! |
| 108 | + ! It is a very fast generator passing BigCrush, and it can be useful if |
| 109 | + ! for some reason you absolutely want 64 bits of state. |
| 110 | + ! |
| 111 | + ! Fortran 90 translated from C by Jim-215-Fisher |
| 112 | + ! |
| 113 | + integer(int64) :: res |
| 114 | + integer(int64), intent(in), optional :: s |
| 115 | + integer(int64) :: int01 = -7046029254386353131_int64, & |
| 116 | + int02 = -4658895280553007687_int64, & |
| 117 | + int03 = -7723592293110705685_int64 |
| 118 | + ! Values are converted from C unsigned integer of 0x9e3779b97f4a7c15, |
| 119 | + ! 0xbf58476d1ce4e5b9, 0x94d049bb133111eb |
| 120 | + |
| 121 | + if(present(s)) si = s |
| 122 | + res = si |
| 123 | + si = res + int01 |
| 124 | + res = ieor(res, shiftr(res, 30)) * int02 |
| 125 | + res = ieor(res, shiftr(res, 27)) * int03 |
| 126 | + res = ieor(res, shiftr(res, 31)) |
| 127 | + end function splitmix64 |
| 128 | + |
| 129 | + #:for k1, t1 in INT_KINDS_TYPES |
| 130 | + subroutine random_distribution_seed_${t1[0]}$${k1}$(put, get) |
| 131 | + !! Set seed value for random number generator |
| 132 | + !! |
| 133 | + ${t1}$, intent(in) :: put |
| 134 | + ${t1}$, intent(out) :: get |
| 135 | + integer(int64) :: tmp |
| 136 | + integer :: i |
| 137 | + |
| 138 | + tmp = splitmix64(int(put, kind = int64)) |
| 139 | + do i = 1, 10 |
| 140 | + tmp = splitmix64( ) |
| 141 | + end do |
| 142 | + do i = 1, 4 |
| 143 | + tmp = splitmix64( ) |
| 144 | + st(i) = tmp |
| 145 | + end do |
| 146 | + get = int(tmp, kind = ${k1}$) |
| 147 | + seed_initialized = .true. |
| 148 | + end subroutine random_distribution_seed_${t1[0]}$${k1}$ |
| 149 | + |
| 150 | + #:endfor |
| 151 | +end module stdlib_stats_distribution_PRNG |
0 commit comments