|
| 1 | +# This file is a part of Julia. License is MIT: https://julialang.org/license |
| 2 | + |
| 3 | +## Xoshiro RNG |
| 4 | +# Lots of implementation is shared with TaskLocalRNG |
| 5 | + |
| 6 | +""" |
| 7 | + Xoshiro |
| 8 | +
|
| 9 | +Xoshiro256++ is a fast pseudorandom number generator originally developed by Sebastian Vigna. |
| 10 | +Reference implementation is available at http://prng.di.unimi.it |
| 11 | +
|
| 12 | +Apart from the high speed, Xoshiro has a small memory footprint, making it suitable for |
| 13 | +applications where many different random states need to be held for long time. |
| 14 | +
|
| 15 | +Julia's Xoshiro implementation has a bulk-generation mode; this seeds new virtual PRNGs |
| 16 | +from the parent, and uses SIMD to generate in parallel (i.e. the bulk stream consists of |
| 17 | +multiple interleaved xoshiro instances). |
| 18 | +The virtual PRNGs are discarded once the bulk request has been serviced (and should cause |
| 19 | +no heap allocations). |
| 20 | +""" |
| 21 | +mutable struct Xoshiro <: AbstractRNG |
| 22 | + s0::UInt64 |
| 23 | + s1::UInt64 |
| 24 | + s2::UInt64 |
| 25 | + s3::UInt64 |
| 26 | +end |
| 27 | + |
| 28 | +Xoshiro(::Nothing) = Xoshiro() |
| 29 | + |
| 30 | +function Xoshiro() |
| 31 | + parent = RandomDevice() |
| 32 | + # Constants have nothing up their sleeve, see task.c |
| 33 | + # 0x02011ce34bce797f == hash(UInt(1))|0x01 |
| 34 | + # 0x5a94851fb48a6e05 == hash(UInt(2))|0x01 |
| 35 | + # 0x3688cf5d48899fa7 == hash(UInt(3))|0x01 |
| 36 | + # 0x867b4bb4c42e5661 == hash(UInt(4))|0x01 |
| 37 | + |
| 38 | + Xoshiro(0x02011ce34bce797f * rand(parent, UInt64), |
| 39 | + 0x5a94851fb48a6e05 * rand(parent, UInt64), |
| 40 | + 0x3688cf5d48899fa7 * rand(parent, UInt64), |
| 41 | + 0x867b4bb4c42e5661 * rand(parent, UInt64)) |
| 42 | +end |
| 43 | + |
| 44 | +copy(rng::Xoshiro) = Xoshiro(rng.s0, rng.s1, rng.s2, rng.s3) |
| 45 | + |
| 46 | +function copy!(dst::Xoshiro, src::Xoshiro) |
| 47 | + dst.s0, dst.s1, dst.s2, dst.s3 = src.s0, src.s1, src.s2, src.s3 |
| 48 | + dst |
| 49 | +end |
| 50 | + |
| 51 | +function ==(a::Xoshiro, b::Xoshiro) |
| 52 | + a.s0 == b.s0 && a.s1 == b.s1 && a.s2 == b.s2 && a.s3 == b.s3 |
| 53 | +end |
| 54 | + |
| 55 | +rng_native_52(::Xoshiro) = UInt64 |
| 56 | + |
| 57 | +function seed!(rng::Xoshiro, s0::UInt64, s1::UInt64, s2::UInt64, s3::UInt64) |
| 58 | + # see task.c |
| 59 | + s = Base.hash_uint64(s0) |
| 60 | + rng.s0 = s |
| 61 | + s += Base.hash_uint64(s1) |
| 62 | + rng.s1 = s |
| 63 | + s += Base.hash_uint64(s2) |
| 64 | + rng.s2 = s |
| 65 | + s += Base.hash_uint64(s3) |
| 66 | + rng.s3 = s |
| 67 | + rng |
| 68 | +end |
| 69 | + |
| 70 | +@inline function rand(rng::Xoshiro, ::SamplerType{UInt64}) |
| 71 | + s0, s1, s2, s3 = rng.s0, rng.s1, rng.s2, rng.s3 |
| 72 | + tmp = s0 + s3 |
| 73 | + res = tmp << 23 | tmp >> 41 |
| 74 | + t = s1 << 17 |
| 75 | + s2 = xor(s2, s0) |
| 76 | + s3 = xor(s3, s1) |
| 77 | + s1 = xor(s1, s2) |
| 78 | + s0 = xor(s0, s3) |
| 79 | + s2 = xor(s2, t) |
| 80 | + s3 = s3 << 45 | s3 >> 19 |
| 81 | + rng.s0, rng.s1, rng.s2, rng.s3 = s0, s1, s2, s3 |
| 82 | + res |
| 83 | +end |
| 84 | + |
| 85 | + |
| 86 | +## Task local RNG |
| 87 | + |
| 88 | +""" |
| 89 | + TaskLocalRNG |
| 90 | +
|
| 91 | +The `TaskLocalRNG` has state that is local to its task, not its thread. |
| 92 | +It is seeded upon task creation, from the state of its parent task. |
| 93 | +Therefore, task creation is an event that changes the parent's RNG state. |
| 94 | +
|
| 95 | +As an upside, the `TaskLocalRNG` is pretty fast, and permits reproducible |
| 96 | +multithreaded simulations (barring race conditions), independent of scheduler |
| 97 | +decisions. As long as the number of threads is not used to make decisions on |
| 98 | +task creation, simulation results are also independent of the number of available |
| 99 | +threads / CPUs. The random stream should not depend on hardware specifics, up to |
| 100 | +endianness and possibly word size. |
| 101 | +
|
| 102 | +Using or seeding the RNG of any other task than the one returned by `current_task()` |
| 103 | +is undefined behavior: it will work most of the time, and may sometimes fail silently. |
| 104 | +""" |
| 105 | +struct TaskLocalRNG <: AbstractRNG end |
| 106 | +TaskLocalRNG(::Nothing) = TaskLocalRNG() |
| 107 | +rng_native_52(::TaskLocalRNG) = UInt64 |
| 108 | + |
| 109 | +function seed!(rng::TaskLocalRNG, s0::UInt64, s1::UInt64, s2::UInt64, s3::UInt64) |
| 110 | + # TODO: Consider a less ad-hoc construction |
| 111 | + # We can afford burning a handful of cycles here, and we don't want any |
| 112 | + # surprises with respect to bad seeds / bad interactions. |
| 113 | + t = current_task() |
| 114 | + s = hash(s0) |
| 115 | + t.rngState0 = s |
| 116 | + s += hash(s1) |
| 117 | + t.rngState1 = s |
| 118 | + s += hash(s2) |
| 119 | + t.rngState2 = s |
| 120 | + s += hash(s3) |
| 121 | + t.rngState3 = s |
| 122 | + rand(rng, UInt64) |
| 123 | + rand(rng, UInt64) |
| 124 | + rand(rng, UInt64) |
| 125 | + rand(rng, UInt64) |
| 126 | + rng |
| 127 | +end |
| 128 | + |
| 129 | +@inline function rand(::TaskLocalRNG, ::SamplerType{UInt64}) |
| 130 | + task = current_task() |
| 131 | + s0, s1, s2, s3 = task.rngState0, task.rngState1, task.rngState2, task.rngState3 |
| 132 | + tmp = s0 + s3 |
| 133 | + res = tmp << 23 | tmp >> 41 |
| 134 | + t = s1 << 17 |
| 135 | + s2 = xor(s2, s0) |
| 136 | + s3 = xor(s3, s1) |
| 137 | + s1 = xor(s1, s2) |
| 138 | + s0 = xor(s0, s3) |
| 139 | + s2 = xor(s2, t) |
| 140 | + s3 = s3 << 45 | s3 >> 19 |
| 141 | + task.rngState0, task.rngState1, task.rngState2, task.rngState3 = s0, s1, s2, s3 |
| 142 | + res |
| 143 | +end |
| 144 | + |
| 145 | +# Shared implementation between Xoshiro and TaskLocalRNG -- seeding |
| 146 | +function seed!(rng::Union{TaskLocalRNG, Xoshiro}, seed::UInt128) |
| 147 | + seed0 = seed % UInt64 |
| 148 | + seed1 = (seed>>>64) % UInt64 |
| 149 | + seed!(rng, seed0, seed1, zero(UInt64), zero(UInt64)) |
| 150 | +end |
| 151 | +seed!(rng::Union{TaskLocalRNG, Xoshiro}, seed::Integer) = seed!(rng, UInt128(seed)) |
| 152 | + |
| 153 | +seed!(rng::Union{TaskLocalRNG, Xoshiro}) = |
| 154 | + seed!(rng, rand(RandomDevice(), UInt64), rand(RandomDevice(), UInt64), |
| 155 | + rand(RandomDevice(), UInt64), rand(RandomDevice(), UInt64)) |
| 156 | + |
| 157 | +function seed!(rng::Union{TaskLocalRNG, Xoshiro}, seed::AbstractVector{UInt64}) |
| 158 | + if length(seed) > 4 |
| 159 | + throw(ArgumentError("seed should have no more than 256 bits")) |
| 160 | + end |
| 161 | + seed0 = length(seed)>0 ? seed[1] : UInt64(0) |
| 162 | + seed1 = length(seed)>1 ? seed[2] : UInt64(0) |
| 163 | + seed2 = length(seed)>2 ? seed[3] : UInt64(0) |
| 164 | + seed3 = length(seed)>3 ? seed[4] : UInt64(0) |
| 165 | + seed!(rng, seed0, seed1, seed2, seed3) |
| 166 | +end |
| 167 | + |
| 168 | +function seed!(rng::Union{TaskLocalRNG, Xoshiro}, seed::AbstractVector{UInt32}) |
| 169 | + if iseven(length(seed)) |
| 170 | + seed!(rng, reinterpret(UInt64, seed)) |
| 171 | + else |
| 172 | + seed!(rng, UInt64[reinterpret(UInt64, @view(seed[begin:end-1])); seed[end] % UInt64]) |
| 173 | + end |
| 174 | +end |
| 175 | + |
| 176 | +@inline function rand(rng::Union{TaskLocalRNG, Xoshiro}, ::SamplerType{UInt128}) |
| 177 | + first = rand(rng, UInt64) |
| 178 | + second = rand(rng,UInt64) |
| 179 | + second + UInt128(first)<<64 |
| 180 | +end |
| 181 | + |
| 182 | +@inline rand(rng::Union{TaskLocalRNG, Xoshiro}, ::SamplerType{Int128}) = rand(rng, UInt128) % Int128 |
| 183 | + |
| 184 | +@inline rand(rng::Union{TaskLocalRNG, Xoshiro}, ::SamplerType{T}) where {T<:Union{Bool, UInt8, Int8, UInt16, Int16, UInt32, Int32, Int64}} = rand(rng, UInt64) % T |
| 185 | + |
| 186 | +function copy(rng::TaskLocalRNG) |
| 187 | + t = current_task() |
| 188 | + Xoshiro(t.rngState0, t.rngState1, t.rngState2, t.rngState3) |
| 189 | +end |
| 190 | + |
| 191 | +function copy!(dst::TaskLocalRNG, src::Xoshiro) |
| 192 | + t = current_task() |
| 193 | + t.rngState0, t.rngState1, t.rngState2, t.rngState3 = src.s0, src.s1, src.s2, src.s3 |
| 194 | + dst |
| 195 | +end |
| 196 | + |
| 197 | +function copy!(dst::Xoshiro, src::TaskLocalRNG) |
| 198 | + t = current_task() |
| 199 | + dst.s0, dst.s1, dst.s2, dst.s3 = t.rngState0, t.rngState1, t.rngState2, t.rngState3 |
| 200 | + dst |
| 201 | +end |
| 202 | + |
| 203 | +function ==(a::Xoshiro, b::TaskLocalRNG) |
| 204 | + t = current_task() |
| 205 | + a.s0 == t.rngState0 && a.s1 == t.rngState1 && a.s2 == t.rngState2 && a.s3 == t.rngState3 |
| 206 | +end |
| 207 | + |
| 208 | +==(a::TaskLocalRNG, b::Xoshiro) = b == a |
0 commit comments