From 1b91b316980d80b2fab9d4e2b8f21bf2e8c1bf63 Mon Sep 17 00:00:00 2001 From: Delyan Angelov Date: Mon, 23 Sep 2019 00:50:22 +0300 Subject: [PATCH] rand: add pcg32 and splitmix64 implementations --- vlib/rand/pcg32.v | 58 +++++++++++++++++++++++++++++++++++++ vlib/rand/pcg32_test.v | 27 +++++++++++++++++ vlib/rand/splitmix64.v | 42 +++++++++++++++++++++++++++ vlib/rand/splitmix64_test.v | 26 +++++++++++++++++ 4 files changed, 153 insertions(+) create mode 100644 vlib/rand/pcg32.v create mode 100644 vlib/rand/pcg32_test.v create mode 100644 vlib/rand/splitmix64.v create mode 100644 vlib/rand/splitmix64_test.v diff --git a/vlib/rand/pcg32.v b/vlib/rand/pcg32.v new file mode 100644 index 0000000000..f2a2954de2 --- /dev/null +++ b/vlib/rand/pcg32.v @@ -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) +} diff --git a/vlib/rand/pcg32_test.v b/vlib/rand/pcg32_test.v new file mode 100644 index 0000000000..2ac64e081d --- /dev/null +++ b/vlib/rand/pcg32_test.v @@ -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] + } +} diff --git a/vlib/rand/splitmix64.v b/vlib/rand/splitmix64.v new file mode 100644 index 0000000000..9ead3dcd4e --- /dev/null +++ b/vlib/rand/splitmix64.v @@ -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) +} diff --git a/vlib/rand/splitmix64_test.v b/vlib/rand/splitmix64_test.v new file mode 100644 index 0000000000..dabae5bd9b --- /dev/null +++ b/vlib/rand/splitmix64_test.v @@ -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] + } +}