JuliaSIMD / VectorizedRNG.jl

Vectorized uniform and normal random samplers.
MIT License
33 stars 7 forks source link

Generate integer numbers #3

Open PallHaraldsson opened 4 years ago

PallHaraldsson commented 4 years ago

Hi,

It seems these (and all?) RNGs generate integer numbers behind the scenes, but I couldn't get them to generate such:

julia> rand(local_rng(), 1:3) ERROR: MethodError: no method matching rng_native_52(::VectorizedRNG.Xoshift{4})

Is it just that floating point is currently implemented (and exported) but integer isn't yet? Or did I overlook how, or even not meant for integer? I found:

"xoroshiro128+ 1.0, our best and fastest small-state generator for floating-point numbers. [..] If you are concerned, use xoroshiro128** or xoshiro256+."

chriselrod commented 4 years ago

I'll need to implement scalar rngs. This library samples many random numbers at a time. To let people sample one at a time, I'll have it store these in a buffer, and then load from the buffer one at a time.

I'll also implement missing methods from the interface, so that it can sample from 1:3 through that.

Under the hood, it samples many random unsigned integers, which it turns into floating point numbers. To get from random integers to 1:3 would also take work, something like:

randrange(r::AbstractUnitRange) = rand(UInt) % length(r) + first(r)

Base might have a more clever algorithm, and there is a PR about using a division-less algorithm.

It's using Xoshiro-256++.

chriselrod commented 4 years ago

Adding random integer support should be really easy. If you want to try making a PR yourself:

Note how random_uniform and random_normal work.

They take a state object, call next_state to advance the state and get random unsigned integers, and then call a function to turn those integers into random uniforms or normals, respectively. To get random integers, just have random_integers return u directly.

Then to get random integers, define a function like randn! that calls random_sample! with your newly defined function as the first argument. You could even just define random_integers as an anonymous function inside the rand!(rng::AbstractVRNG, ::AbstractArray{Int64}) method if you want.

If I recall correctly, u is a NTuple{N,NTuple{W,Core.VecElement{UInt64}}}. If you want to do something more complicated, you can turn the individual NTuple{W,Core.VecElement{UInt64}} into SVecs using SIMDPirates, which VectorizedRNG already depends on. The SVec type overloads a lot of base functions.

chriselrod commented 4 years ago

This what you wanted?

julia> x = zeros(UInt64, 200);

julia> @benchmark rand!(local_rng(), $x)
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     34.476 ns (0.00% GC)
  median time:      34.513 ns (0.00% GC)
  mean time:        34.557 ns (0.00% GC)
  maximum time:     67.631 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     993

julia> x'
1×200 LinearAlgebra.Adjoint{UInt64,Array{UInt64,1}}:
 0x9dca86b0f67b8940  0xd2b7260f7d4c405a  0xf6e9e52c2393b636  0xb588fe6931404fe7  0x15ecca994eedba28  …  0x36dbf4cadaadb2d1  0x35ab7b37a6689a18  0x824ccbff468599c4  0x7c9d89641210d3ff