rand: add pcg32 and splitmix64 implementations
parent
ab528bb043
commit
1b91b31698
|
@ -0,0 +1,58 @@
|
||||||
|
module rand
|
||||||
|
|
||||||
|
// Ported from http://www.pcg-random.org/download.html
|
||||||
|
// and https://github.com/imneme/pcg-c-basic/blob/master/pcg_basic.c
|
||||||
|
|
||||||
|
struct Pcg32 {
|
||||||
|
mut:
|
||||||
|
state u64
|
||||||
|
inc u64
|
||||||
|
}
|
||||||
|
/**
|
||||||
|
* new_pcg32 - a Pcg32 PRNG generator
|
||||||
|
* @param initstate - the initial state of the PRNG.
|
||||||
|
* @param initseq - the stream/step of the PRNG.
|
||||||
|
* @return a new Pcg32 PRNG instance
|
||||||
|
*/
|
||||||
|
pub fn new_pcg32(initstate u64, initseq u64) Pcg32 {
|
||||||
|
mut rng := Pcg32{}
|
||||||
|
rng.state = u64(0)
|
||||||
|
rng.inc = u64( u64(initseq << u64(1)) | u64(1) )
|
||||||
|
rng.next()
|
||||||
|
rng.state += initstate
|
||||||
|
rng.next()
|
||||||
|
return rng
|
||||||
|
}
|
||||||
|
/**
|
||||||
|
* Pcg32.next - update the PRNG state and get back the next random number
|
||||||
|
* @return the generated pseudo random number
|
||||||
|
*/
|
||||||
|
[inline] pub fn (rng mut Pcg32) next() u32 {
|
||||||
|
oldstate := rng.state
|
||||||
|
rng.state = oldstate * u64(6364136223846793005) + rng.inc
|
||||||
|
xorshifted := u32( u64( u64(oldstate >> u64(18)) ^ oldstate) >> u64(27) )
|
||||||
|
rot := u32( oldstate >> u64(59) )
|
||||||
|
return u32( (xorshifted >> rot) | (xorshifted << ((-rot) & u32(31))) )
|
||||||
|
}
|
||||||
|
/**
|
||||||
|
* Pcg32.bounded_next - update the PRNG state. Get the next number < bound
|
||||||
|
* @param bound - the returned random number will be < bound
|
||||||
|
* @return the generated pseudo random number
|
||||||
|
*/
|
||||||
|
[inline] pub fn (rng mut Pcg32) bounded_next(bound u32) u32 {
|
||||||
|
// To avoid bias, we need to make the range of the RNG a multiple of
|
||||||
|
// bound, which we do by dropping output less than a threshold.
|
||||||
|
threshold := u32( -bound % bound )
|
||||||
|
// Uniformity guarantees that loop below will terminate. In practice, it
|
||||||
|
// should usually terminate quickly; on average (assuming all bounds are
|
||||||
|
// equally likely), 82.25% of the time, we can expect it to require just
|
||||||
|
// one iteration. In practice, bounds are typically small and only a
|
||||||
|
// tiny amount of the range is eliminated.
|
||||||
|
for {
|
||||||
|
r := rng.next()
|
||||||
|
if r >= threshold {
|
||||||
|
return u32( r % bound )
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return u32(0)
|
||||||
|
}
|
|
@ -0,0 +1,27 @@
|
||||||
|
|
||||||
|
import rand
|
||||||
|
import time
|
||||||
|
|
||||||
|
fn gen_randoms(initstate i64, initseq i64, bound int) []u32 {
|
||||||
|
mut randoms := [u32(0)].repeat(20)
|
||||||
|
mut rnd := rand.new_pcg32( u64(initstate), u64(initseq) )
|
||||||
|
for i in 0..20 {
|
||||||
|
randoms[i] = rnd.bounded_next(u32(bound))
|
||||||
|
}
|
||||||
|
return randoms
|
||||||
|
}
|
||||||
|
|
||||||
|
fn test_pcg32_reproducibility() {
|
||||||
|
t := time.ticks()
|
||||||
|
tseq := t % 23237671
|
||||||
|
println('t: $t | tseq: $tseq')
|
||||||
|
randoms1 := gen_randoms(t, tseq, 1000)
|
||||||
|
randoms2 := gen_randoms(t, tseq, 1000)
|
||||||
|
assert randoms1.len == randoms2.len
|
||||||
|
println( randoms1 )
|
||||||
|
println( randoms2 )
|
||||||
|
len := randoms1.len
|
||||||
|
for i in 0..len {
|
||||||
|
assert randoms1[i] == randoms2[i]
|
||||||
|
}
|
||||||
|
}
|
|
@ -0,0 +1,42 @@
|
||||||
|
module rand
|
||||||
|
|
||||||
|
// Ported from http://xoshiro.di.unimi.it/splitmix64.c
|
||||||
|
|
||||||
|
struct Splitmix64 {
|
||||||
|
mut:
|
||||||
|
state u64
|
||||||
|
}
|
||||||
|
/**
|
||||||
|
* new_splitmix64 - a Splitmix64 PRNG generator
|
||||||
|
* @param seed the initial seed of the PRNG.
|
||||||
|
* @return a new Splitmix64 PRNG instance
|
||||||
|
*/
|
||||||
|
pub fn new_splitmix64(seed u64) Splitmix64 {
|
||||||
|
return Splitmix64{ seed }
|
||||||
|
}
|
||||||
|
/**
|
||||||
|
* Splitmix64.next - update the PRNG state and get back the next random number
|
||||||
|
* @return the generated pseudo random number
|
||||||
|
*/
|
||||||
|
[inline] pub fn (rng mut Splitmix64) next() u64 {
|
||||||
|
rng.state += u64(0x9e3779b97f4a7c15)
|
||||||
|
mut z := rng.state
|
||||||
|
z = (z ^ u64((z >> u64(30)))) * u64(0xbf58476d1ce4e5b9)
|
||||||
|
z = (z ^ u64((z >> u64(27)))) * u64(0x94d049bb133111eb)
|
||||||
|
return z ^ u64(z >> u64(31))
|
||||||
|
}
|
||||||
|
/**
|
||||||
|
* Splitmix64.bounded_next - Get the next random number < bound
|
||||||
|
* @param bound - the returned random number will be < bound
|
||||||
|
* @return the generated pseudo random number
|
||||||
|
*/
|
||||||
|
[inline] pub fn (rng mut Splitmix64) bounded_next(bound u64) u64 {
|
||||||
|
threshold := u64( -bound % bound )
|
||||||
|
for {
|
||||||
|
r := rng.next()
|
||||||
|
if r >= threshold {
|
||||||
|
return u64( r % bound )
|
||||||
|
}
|
||||||
|
}
|
||||||
|
return u64(0)
|
||||||
|
}
|
|
@ -0,0 +1,26 @@
|
||||||
|
|
||||||
|
import rand
|
||||||
|
import time
|
||||||
|
|
||||||
|
fn gen_randoms(seed i64, bound int) []u64 {
|
||||||
|
mut randoms := [u64(0)].repeat(20)
|
||||||
|
mut rnd := rand.new_splitmix64( u64(seed) )
|
||||||
|
for i in 0..20 {
|
||||||
|
randoms[i] = rnd.bounded_next(u64(bound))
|
||||||
|
}
|
||||||
|
return randoms
|
||||||
|
}
|
||||||
|
|
||||||
|
fn test_splitmix64_reproducibility() {
|
||||||
|
t := time.ticks()
|
||||||
|
println('t: $t')
|
||||||
|
randoms1 := gen_randoms(t, 1000)
|
||||||
|
randoms2 := gen_randoms(t, 1000)
|
||||||
|
assert randoms1.len == randoms2.len
|
||||||
|
println( randoms1 )
|
||||||
|
println( randoms2 )
|
||||||
|
len := randoms1.len
|
||||||
|
for i in 0..len {
|
||||||
|
assert randoms1[i] == randoms2[i]
|
||||||
|
}
|
||||||
|
}
|
Loading…
Reference in New Issue