From a7c84834f4ba4948b8102a05b603bf8d51484760 Mon Sep 17 00:00:00 2001 From: Hungry Blue Dev <46458389+hungrybluedev@users.noreply.github.com> Date: Tue, 2 Jun 2020 00:43:56 +0530 Subject: [PATCH] rand: reorganize (step 1) --- vlib/hash/wyhash/rand.v | 15 -- vlib/hash/wyhash/wyhash.v | 2 +- vlib/hash/wyhash/wyhash_test.v | 13 -- vlib/rand/mt19937.v | 320 +++++++++++++++++++++++++++++ vlib/rand/mt19937_test.v | 340 ++++++++++++++++++++++++++++++ vlib/rand/musl_rng.v | 236 +++++++++++++++++++++ vlib/rand/musl_rng_test.v | 330 +++++++++++++++++++++++++++++ vlib/rand/pcg32.v | 260 ++++++++++++++++++----- vlib/rand/pcg32_test.v | 340 +++++++++++++++++++++++++++--- vlib/rand/rand.v | 133 +++++++++++- vlib/rand/random_numbers_test.v | 151 +++++++++++++- vlib/rand/splitmix64.v | 232 ++++++++++++++++++--- vlib/rand/splitmix64_test.v | 338 ++++++++++++++++++++++++++++-- vlib/rand/system_rng.c.v | 291 ++++++++++++++++++++++++++ vlib/rand/system_rng.js.v | 15 ++ vlib/rand/system_rng_test.v | 354 ++++++++++++++++++++++++++++++++ vlib/rand/util.v | 42 ++++ vlib/rand/wyrand.v | 251 ++++++++++++++++++++++ vlib/rand/wyrand_test.v | 330 +++++++++++++++++++++++++++++ 19 files changed, 3833 insertions(+), 160 deletions(-) delete mode 100644 vlib/hash/wyhash/rand.v create mode 100644 vlib/rand/mt19937.v create mode 100644 vlib/rand/mt19937_test.v create mode 100644 vlib/rand/musl_rng.v create mode 100644 vlib/rand/musl_rng_test.v create mode 100644 vlib/rand/system_rng.c.v create mode 100644 vlib/rand/system_rng.js.v create mode 100644 vlib/rand/system_rng_test.v create mode 100644 vlib/rand/util.v create mode 100644 vlib/rand/wyrand.v create mode 100644 vlib/rand/wyrand_test.v diff --git a/vlib/hash/wyhash/rand.v b/vlib/hash/wyhash/rand.v deleted file mode 100644 index 0583a905dc..0000000000 --- a/vlib/hash/wyhash/rand.v +++ /dev/null @@ -1,15 +0,0 @@ -// Copyright (c) 2019-2020 Alexander Medvednikov. All rights reserved. -// Use of this source code is governed by an MIT license -// that can be found in the LICENSE file. -module wyhash - -pub fn rand_u64(seed &u64) u64 { - mut seed0 := seed - unsafe{ - mut seed1 := *seed0 - seed1 += wyp0 - *seed0 = seed1 - return wymum(seed1^wyp1, seed1) - } - return 0 -} diff --git a/vlib/hash/wyhash/wyhash.v b/vlib/hash/wyhash/wyhash.v index d2978a67c0..d4748b391c 100644 --- a/vlib/hash/wyhash/wyhash.v +++ b/vlib/hash/wyhash/wyhash.v @@ -92,7 +92,7 @@ fn wyrotr(v u64, k u32) u64 { } [inline] -fn wymum(a, b u64) u64 { +pub fn wymum(a, b u64) u64 { /* mut r := u128(a) r = r*b diff --git a/vlib/hash/wyhash/wyhash_test.v b/vlib/hash/wyhash/wyhash_test.v index 75c890c185..1d7bd3b675 100644 --- a/vlib/hash/wyhash/wyhash_test.v +++ b/vlib/hash/wyhash/wyhash_test.v @@ -27,16 +27,3 @@ fn test_wyhash() { assert got == test.expected } } - -fn test_rand_u64() { - seed := u64(111) - mut rand_nos := []u64{} - for _ in 0..40 { - rand_no := wyhash.rand_u64(&seed) - for r in rand_nos { - assert rand_no != r - } - rand_nos << rand_no - } - assert true -} diff --git a/vlib/rand/mt19937.v b/vlib/rand/mt19937.v new file mode 100644 index 0000000000..3fa7bf239b --- /dev/null +++ b/vlib/rand/mt19937.v @@ -0,0 +1,320 @@ +// Copyright (c) 2019-2020 Alexander Medvednikov. All rights reserved. +// Use of this source code is governed by an MIT license +// that can be found in the LICENSE file. +module rand + +import math.bits + +/* +C++ functions for MT19937, with initialization improved 2002/2/10. + Coded by Takuji Nishimura and Makoto Matsumoto. + This is a faster version by taking Shawn Cokus's optimization, + Matthe Bellew's simplification, Isaku Wada's real version. + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + + Any feedback is very welcome. + http://www.math.keio.ac.jp/matumoto/emt.html + email: matumoto@math.keio.ac.jp +*/ +const ( + nn = 312 + mm = 156 + matrix_a = 0xB5026F5AA96619E9 + um = 0xFFFFFFFF80000000 + lm = 0x7FFFFFFF + inv_f64_limit = 1.0 / 9007199254740992.0 +) + +// A generator that uses the Mersenne Twister algorithm with period 2^19937 +pub struct MT19937RNG { +mut: + state []u64 = calculate_state(time_seed_array(2), mut []u64{len: nn}) + mti int = nn + next_rnd u32 = 0 + has_next bool = false +} + +fn calculate_state(seed_data []u32, mut state []u64) []u64 { + lo := u64(seed_data[0]) + hi := u64(seed_data[1]) + state[0] = u64((hi << 32) | lo) + for j := 1; j < nn; j++ { + state[j] = u64(6364136223846793005) * (state[j - 1] ^ (state[j - 1] >> 62)) + u64(j) + } + return state +} + +// seed() - Set the seed, needs only two u32s in little endian format as [lower, higher] +pub fn (mut rng MT19937RNG) seed(seed_data []u32) { + if seed_data.len != 2 { + eprintln('mt19937 needs only two 32bit integers as seed: [lower, higher]') + exit(1) + } + rng.state = calculate_state(seed_data, mut rng.state) + rng.mti = nn + rng.next_rnd = 0 + rng.has_next = false +} + +// rng.u32() - return a pseudorandom 32bit int in [0, 2**32) +[inline] +pub fn (mut rng MT19937RNG) u32() u32 { + if rng.has_next { + rng.has_next = false + return rng.next_rnd + } + ans := rng.u64() + rng.next_rnd = u32(ans >> 32) + rng.has_next = true + return u32(ans & 0xffffffff) +} + +// rng.u64() - return a pseudorandom 64bit int in [0, 2**64) +[inline] +pub fn (mut rng MT19937RNG) u64() u64 { + mag01 := [u64(0), u64(matrix_a)] + mut x := u64(0) + mut i := int(0) + if rng.mti >= nn { + for i = 0; i < nn - mm; i++ { + x = (rng.state[i] & um) | (rng.state[i + 1] & lm) + rng.state[i] = rng.state[i + mm] ^ (x >> 1) ^ mag01[int(x & 1)] + } + for i < nn - 1 { + x = (rng.state[i] & um) | (rng.state[i + 1] & lm) + rng.state[i] = rng.state[i + (mm - nn)] ^ (x >> 1) ^ mag01[int(x & 1)] + i++ + } + x = (rng.state[nn - 1] & um) | (rng.state[0] & lm) + rng.state[nn - 1] = rng.state[mm - 1] ^ (x >> 1) ^ mag01[int(x & 1)] + rng.mti = 0 + } + x = rng.state[rng.mti] + rng.mti++ + x ^= (x >> 29) & 0x5555555555555555 + x ^= (x << 17) & 0x71D67FFFEDA60000 + x ^= (x << 37) & 0xFFF7EEE000000000 + x ^= (x >> 43) + return x +} + +// rng.int() - return a 32-bit signed (possibly negative) int +[inline] +pub fn (mut rng MT19937RNG) int() int { + return int(rng.u32()) +} + +// rng.i64() - return a 64-bit signed (possibly negative) i64 +[inline] +pub fn (mut rng MT19937RNG) i64() i64 { + return i64(rng.u64()) +} + +// rng.int31() - return a 31bit positive pseudorandom integer +[inline] +pub fn (mut rng MT19937RNG) int31() int { + return int(rng.u32() >> 1) +} + +// rng.int63() - return a 63bit positive pseudorandom integer +[inline] +pub fn (mut rng MT19937RNG) int63() i64 { + return i64(rng.u64() >> 1) +} + +// rng.u32n(max) - return a 32bit u32 in [0, max) +[inline] +pub fn (mut rng MT19937RNG) u32n(max u32) u32 { + if max == 0 { + eprintln('max must be positive integer') + exit(1) + } + // Check SysRNG in system_rng.c.v for explanation + bit_len := bits.len_32(max) + if bit_len == 32 { + for { + value := rng.u32() + if value < max { + return value + } + } + } else { + mask := (u32(1) << (bit_len + 1)) - 1 + for { + value := rng.u32() & mask + if value < max { + return value + } + } + } + return u32(0) +} + +// rng.u64n(max) - return a 64bit u64 in [0, max) +[inline] +pub fn (mut rng MT19937RNG) u64n(max u64) u64 { + if max == 0 { + eprintln('max must be positive integer') + exit(1) + } + bit_len := bits.len_64(max) + if bit_len == 64 { + for { + value := rng.u64() + if value < max { + return value + } + } + } else { + mask := (u64(1) << (bit_len + 1)) - 1 + for { + value := rng.u64() & mask + if value < max { + return value + } + } + } + return u64(0) +} + +// rng.u32n(min, max) returns a pseudorandom u32 value that is guaranteed to be in [min, max) +[inline] +pub fn (mut rng MT19937RNG) u32_in_range(min, max u32) u32 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.u32n(max - min) +} + +// rng.u64n(min, max) returns a pseudorandom u64 value that is guaranteed to be in [min, max) +[inline] +pub fn (mut rng MT19937RNG) u64_in_range(min, max u64) u64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.u64n(max - min) +} + +// rng.intn(max) - return a 32bit positive int in [0, max) +[inline] +pub fn (mut rng MT19937RNG) intn(max int) int { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return int(rng.u32n(max)) +} + +// rng.i64n(max) - return a 64bit positive i64 in [0, max) +[inline] +pub fn (mut rng MT19937RNG) i64n(max i64) i64 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return i64(rng.u64n(max)) +} + +// rng.int_in_range(min, max) - return a 32bit positive int in [0, max) +[inline] +pub fn (mut rng MT19937RNG) int_in_range(min, max int) int { + if max <= min { + eprintln('max must be greater than min.') + exit(1) + } + return min + rng.intn(max - min) +} + +// rng.i64_in_range(min, max) - return a 64bit positive i64 in [0, max) +[inline] +pub fn (mut rng MT19937RNG) i64_in_range(min, max i64) i64 { + if max <= min { + eprintln('max must be greater than min.') + exit(1) + } + return min + rng.i64n(max - min) +} + +// rng.f32() - return a 32bit real in [0, 1) +[inline] +pub fn (mut rng MT19937RNG) f32() f32 { + return f32(rng.f64()) +} + +// rng.f64() - return 64bit real in [0, 1) +[inline] +pub fn (mut rng MT19937RNG) f64() f64 { + return f64(rng.u64() >> 11) * inv_f64_limit +} + +// rng.f32n(max) - return 64bit real in [0, max) +[inline] +pub fn (mut rng MT19937RNG) f32n(max f32) f32 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return rng.f32() * max +} + +// rng.f64n(max) - return 64bit real in [0, max) +[inline] +pub fn (mut rng MT19937RNG) f64n(max f64) f64 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return rng.f64() * max +} + +// rng.f32_in_range(min, max) returns a pseudorandom f32 that lies in [min, max) +[inline] +pub fn (mut rng MT19937RNG) f32_in_range(min, max f32) f32 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.f32n(max - min) +} + +// rng.i64_in_range(min, max) returns a pseudorandom i64 that lies in [min, max) +[inline] +pub fn (mut rng MT19937RNG) f64_in_range(min, max f64) f64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.f64n(max - min) +} diff --git a/vlib/rand/mt19937_test.v b/vlib/rand/mt19937_test.v new file mode 100644 index 0000000000..a392073b9c --- /dev/null +++ b/vlib/rand/mt19937_test.v @@ -0,0 +1,340 @@ +import rand +import math + +const ( + range_limit = 40 + value_count = 1000 + seeds = [[u32(0xcafebabe), u32(0xdeadbeef)], [u32(0xc0de), u32(0xfeed)]] +) + +const ( + sample_size = 1000 + stats_epsilon = 0.05 + inv_sqrt_12 = 1.0 / math.sqrt(12) +) + +fn mt19937_basic_test() { + rng := rand.MT19937RNG{} + rng.seed([u32(0xdeadbeef)]) + target := [956529277, 3842322136, 3319553134, 1843186657, 2704993644, 595827513, 938518626, + 1676224337, 3221315650, 1819026461] + for i := 0; i < 10; i++ { + assert target[i] == rng.u32() + } +} + +fn gen_randoms(seed_data []u32, bound int) []u64 { + bound_u64 := u64(bound) + mut randoms := [u64(0)].repeat(20) + mut rnd := rand.MT19937RNG{} + rnd.seed(seed_data) + for i in 0 .. 20 { + randoms[i] = rnd.u64n(bound_u64) + } + return randoms +} + +fn test_mt19937_reproducibility() { + seed_data := rand.time_seed_array(2) + randoms1 := gen_randoms(seed_data, 1000) + randoms2 := gen_randoms(seed_data, 1000) + assert randoms1.len == randoms2.len + len := randoms1.len + for i in 0 .. len { + assert randoms1[i] == randoms2[i] + } +} + +// TODO: use the `in` syntax and remove this function +// after generics has been completely implemented +fn found(value u64, arr []u64) bool { + for item in arr { + if value == item { + return true + } + } + return false +} + +fn test_mt19937_variability() { + // If this test fails and if it is certainly not the implementation + // at fault, try changing the seed values. Repeated values are + // improbable but not impossible. + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + mut values := []u64{cap: value_count} + for i in 0 .. value_count { + value := rng.u64() + assert !found(value, values) + assert values.len == i + values << value + } + } +} + +fn check_uniformity_u64(rng rand.MT19937RNG, range u64) { + range_f64 := f64(range) + expected_mean := range_f64 / 2.0 + mut variance := 0.0 + for _ in 0 .. sample_size { + diff := f64(rng.u64n(range)) - expected_mean + variance += diff * diff + } + variance /= sample_size - 1 + sigma := math.sqrt(variance) + expected_sigma := range_f64 * inv_sqrt_12 + error := (sigma - expected_sigma) / expected_sigma + assert math.abs(error) < stats_epsilon +} + +fn test_mt19937_uniformity_u64() { + ranges := [14019545, 80240, 130] + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + for range in ranges { + check_uniformity_u64(rng, range) + } + } +} + +fn check_uniformity_f64(rng rand.MT19937RNG) { + expected_mean := 0.5 + mut variance := 0.0 + for _ in 0 .. sample_size { + diff := rng.f64() - expected_mean + variance += diff * diff + } + variance /= sample_size - 1 + sigma := math.sqrt(variance) + expected_sigma := inv_sqrt_12 + error := (sigma - expected_sigma) / expected_sigma + assert math.abs(error) < stats_epsilon +} + +fn test_mt19937_uniformity_f64() { + // The f64 version + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + check_uniformity_f64(rng) + } +} + +fn test_mt19937_u32n() { + max := 16384 + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u32n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_mt19937_u64n() { + max := u64(379091181005) + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u64n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_mt19937_u32_in_range() { + max := 484468466 + min := 316846 + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u32_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_mt19937_u64_in_range() { + max := u64(216468454685163) + min := u64(6848646868) + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u64_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_mt19937_int31() { + max_u31 := 0x7FFFFFFF + sign_mask := 0x80000000 + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int31() + assert value >= 0 + assert value <= max_u31 + // This statement ensures that the sign bit is zero + assert (value & sign_mask) == 0 + } + } +} + +fn test_mt19937_int63() { + max_u63 := i64(0x7FFFFFFFFFFFFFFF) + sign_mask := i64(0x8000000000000000) + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int63() + assert value >= 0 + assert value <= max_u63 + assert (value & sign_mask) == 0 + } + } +} + +fn test_mt19937_intn() { + max := 2525642 + for seed in seeds { + rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.intn(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_mt19937_i64n() { + max := i64(3246727724653636) + for seed in seeds { + rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.i64n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_mt19937_int_in_range() { + min := -4252 + max := 1034 + for seed in seeds { + rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_mt19937_i64_in_range() { + min := i64(-24095) + max := i64(324058) + for seed in seeds { + rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.i64_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_mt19937_f32() { + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= 0.0 + assert value < 1.0 + } + } +} + +fn test_mt19937_f64() { + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= 0.0 + assert value < 1.0 + } + } +} + +fn test_mt19937_f32n() { + max := f32(357.0) + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32n(max) + assert value >= 0.0 + assert value < max + } + } +} + +fn test_mt19937_f64n() { + max := 1.52e6 + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64n(max) + assert value >= 0.0 + assert value < max + } + } +} + +fn test_mt19937_f32_in_range() { + min := f32(-24.0) + max := f32(125.0) + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_mt19937_f64_in_range() { + min := -548.7 + max := 5015.2 + for seed in seeds { + mut rng := rand.MT19937RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64_in_range(min, max) + assert value >= min + assert value < max + } + } +} diff --git a/vlib/rand/musl_rng.v b/vlib/rand/musl_rng.v new file mode 100644 index 0000000000..0cb05ec88a --- /dev/null +++ b/vlib/rand/musl_rng.v @@ -0,0 +1,236 @@ +// Copyright (c) 2019-2020 Alexander Medvednikov. All rights reserved. +// Use of this source code is governed by an MIT license +// that can be found in the LICENSE file. +module rand + +import math.bits + +// Ported from https://git.musl-libc.org/cgit/musl/tree/src/prng/rand_r.c +pub struct MuslRNG { +mut: + state u32 = time_seed_32() +} + +pub fn (mut rng MuslRNG) seed(seed_data []u32) { + if seed_data.len != 1 { + eprintln('MuslRNG needs only one unsigned 32 bit integer as a seed.') + exit(1) + } + rng.state = seed_data[0] +} + +[inline] +fn temper(prev u32) u32 { + mut x := prev + x ^= x >> 11 + x ^= (x << 7) & 0x9D2C5680 + x ^= (x << 15) & 0xEFC60000 + x ^= (x >> 18) + return x +} + +// rng.u32() - return a pseudorandom 32 bit unsigned u32 +[inline] +pub fn (mut rng MuslRNG) u32() u32 { + rng.state = rng.state * 1103515245 + 12345 + // We are not dividing by 2 (or shifting right by 1) + // because we want all 32-bits of random data + return temper(rng.state) +} + +// rng.u64() - return a pseudorandom 64 bit unsigned u64 +[inline] +pub fn (mut rng MuslRNG) u64() u64 { + return u64(rng.u32()) | (u64(rng.u32()) << 32) +} + +// rn.u32n(max) - return a pseudorandom 32 bit unsigned u32 in [0, max) +[inline] +pub fn (mut rng MuslRNG) u32n(max u32) u32 { + if max == 0 { + eprintln('max must be positive integer') + exit(1) + } + // Check SysRNG in system_rng.c.v for explanation + bit_len := bits.len_32(max) + if bit_len == 32 { + for { + value := rng.u32() + if value < max { + return value + } + } + } else { + mask := (u32(1) << (bit_len + 1)) - 1 + for { + value := rng.u32() & mask + if value < max { + return value + } + } + } + return u32(0) +} + +// rn.u64n(max) - return a pseudorandom 64 bit unsigned u64 in [0, max) +[inline] +pub fn (mut rng MuslRNG) u64n(max u64) u64 { + if max == 0 { + eprintln('max must be positive integer') + exit(1) + } + bit_len := bits.len_64(max) + if bit_len == 64 { + for { + value := rng.u64() + if value < max { + return value + } + } + } else { + mask := (u64(1) << (bit_len + 1)) - 1 + for { + value := rng.u64() & mask + if value < max { + return value + } + } + } + return u64(0) +} + +// rn.u32_in_range(min, max) - return a pseudorandom 32 bit unsigned u32 in [min, max) +[inline] +pub fn (mut rng MuslRNG) u32_in_range(min, max u64) u64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.u32n(max - min) +} + +// rn.u64_in_range(min, max) - return a pseudorandom 64 bit unsigned u64 in [min, max) +[inline] +pub fn (mut rng MuslRNG) u64_in_range(min, max u64) u64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.u64n(max - min) +} + +// rng.int() - return a 32-bit signed (possibly negative) int +[inline] +pub fn (mut rng MuslRNG) int() int { + return int(rng.u32()) +} + +// rng.i64() - return a 64-bit signed (possibly negative) i64 +[inline] +pub fn (mut rng MuslRNG) i64() i64 { + return i64(rng.u64()) +} + +// rng.int31() - return a 31bit positive pseudorandom integer +[inline] +pub fn (mut rng MuslRNG) int31() int { + return int(rng.u32() >> 1) +} + +// rng.int63() - return a 63bit positive pseudorandom integer +[inline] +pub fn (mut rng MuslRNG) int63() i64 { + return i64(rng.u64() >> 1) +} + +// rng.intn(max) - return a 32bit positive int in [0, max) +[inline] +pub fn (mut rng MuslRNG) intn(max int) int { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return int(rng.u32n(max)) +} + +// rng.i64n(max) - return a 64bit positive i64 in [0, max) +[inline] +pub fn (mut rng MuslRNG) i64n(max i64) i64 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return i64(rng.u64n(max)) +} + +// rng.int_in_range(min, max) - return a 32bit positive int in [0, max) +[inline] +pub fn (mut rng MuslRNG) int_in_range(min, max int) int { + if max <= min { + eprintln('max must be greater than min.') + exit(1) + } + return min + rng.intn(max - min) +} + +// rng.i64_in_range(min, max) - return a 64bit positive i64 in [0, max) +[inline] +pub fn (mut rng MuslRNG) i64_in_range(min, max i64) i64 { + if max <= min { + eprintln('max must be greater than min.') + exit(1) + } + return min + rng.i64n(max - min) +} + +// rng.f32() returns a pseudorandom f32 value between 0.0 (inclusive) and 1.0 (exclusive) i.e [0, 1) +[inline] +pub fn (mut rng MuslRNG) f32() f32 { + return f32(rng.u32()) / max_u32_as_f32 +} + +// rng.f64() returns a pseudorandom f64 value between 0.0 (inclusive) and 1.0 (exclusive) i.e [0, 1) +[inline] +pub fn (mut rng MuslRNG) f64() f64 { + return f64(rng.u64()) / max_u64_as_f64 +} + +// rng.f32n() returns a pseudorandom f32 value in [0, max) +[inline] +pub fn (mut rng MuslRNG) f32n(max f32) f32 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return rng.f32() * max +} + +// rng.f64n() returns a pseudorandom f64 value in [0, max) +[inline] +pub fn (mut rng MuslRNG) f64n(max f64) f64 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return rng.f64() * max +} + +// rng.f32_in_range(min, max) returns a pseudorandom f32 that lies in [min, max) +[inline] +pub fn (mut rng MuslRNG) f32_in_range(min, max f32) f32 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.f32n(max - min) +} + +// rng.i64_in_range(min, max) returns a pseudorandom i64 that lies in [min, max) +[inline] +pub fn (mut rng MuslRNG) f64_in_range(min, max f64) f64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.f64n(max - min) +} diff --git a/vlib/rand/musl_rng_test.v b/vlib/rand/musl_rng_test.v new file mode 100644 index 0000000000..45c8d86a91 --- /dev/null +++ b/vlib/rand/musl_rng_test.v @@ -0,0 +1,330 @@ +import rand +import math + +const ( + range_limit = 40 + value_count = 1000 + seeds = [[u32(42)], [u32(256)]] +) + +const ( + sample_size = 1000 + stats_epsilon = 0.05 + inv_sqrt_12 = 1.0 / math.sqrt(12) +) + +fn gen_randoms(seed_data []u32, bound int) []u64 { + bound_u64 := u64(bound) + mut randoms := [u64(0)].repeat(20) + mut rnd := rand.MuslRNG{} + rnd.seed(seed_data) + for i in 0 .. 20 { + randoms[i] = rnd.u64n(bound_u64) + } + return randoms +} + +fn test_musl_reproducibility() { + seed_data := rand.time_seed_array(1) + randoms1 := gen_randoms(seed_data, 1000) + randoms2 := gen_randoms(seed_data, 1000) + assert randoms1.len == randoms2.len + len := randoms1.len + for i in 0 .. len { + assert randoms1[i] == randoms2[i] + } +} + +// TODO: use the `in` syntax and remove this function +// after generics has been completely implemented +fn found(value u64, arr []u64) bool { + for item in arr { + if value == item { + return true + } + } + return false +} + +fn test_musl_variability() { + // If this test fails and if it is certainly not the implementation + // at fault, try changing the seed values. Repeated values are + // improbable but not impossible. + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + mut values := []u64{cap: value_count} + for i in 0 .. value_count { + value := rng.u64() + assert !found(value, values) + assert values.len == i + values << value + } + } +} + +fn check_uniformity_u64(rng rand.MuslRNG, range u64) { + range_f64 := f64(range) + expected_mean := range_f64 / 2.0 + mut variance := 0.0 + for _ in 0 .. sample_size { + diff := f64(rng.u64n(range)) - expected_mean + variance += diff * diff + } + variance /= sample_size - 1 + sigma := math.sqrt(variance) + expected_sigma := range_f64 * inv_sqrt_12 + error := (sigma - expected_sigma) / expected_sigma + assert math.abs(error) < stats_epsilon +} + +fn test_musl_uniformity_u64() { + ranges := [14019545, 80240, 130] + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + for range in ranges { + check_uniformity_u64(rng, range) + } + } +} + +fn check_uniformity_f64(rng rand.MuslRNG) { + expected_mean := 0.5 + mut variance := 0.0 + for _ in 0 .. sample_size { + diff := rng.f64() - expected_mean + variance += diff * diff + } + variance /= sample_size - 1 + sigma := math.sqrt(variance) + expected_sigma := inv_sqrt_12 + error := (sigma - expected_sigma) / expected_sigma + assert math.abs(error) < stats_epsilon +} + +fn test_musl_uniformity_f64() { + // The f64 version + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + check_uniformity_f64(rng) + } +} + +fn test_musl_u32n() { + max := 16384 + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u32n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_musl_u64n() { + max := u64(379091181005) + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u64n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_musl_u32_in_range() { + max := 484468466 + min := 316846 + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u32_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_musl_u64_in_range() { + max := u64(216468454685163) + min := u64(6848646868) + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u64_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_musl_int31() { + max_u31 := 0x7FFFFFFF + sign_mask := 0x80000000 + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int31() + assert value >= 0 + assert value <= max_u31 + // This statement ensures that the sign bit is zero + assert (value & sign_mask) == 0 + } + } +} + +fn test_musl_int63() { + max_u63 := i64(0x7FFFFFFFFFFFFFFF) + sign_mask := i64(0x8000000000000000) + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int63() + assert value >= 0 + assert value <= max_u63 + assert (value & sign_mask) == 0 + } + } +} + +fn test_musl_intn() { + max := 2525642 + for seed in seeds { + rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.intn(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_musl_i64n() { + max := i64(3246727724653636) + for seed in seeds { + rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.i64n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_musl_int_in_range() { + min := -4252 + max := 1034 + for seed in seeds { + rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_musl_i64_in_range() { + min := i64(-24095) + max := i64(324058) + for seed in seeds { + rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.i64_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_musl_f32() { + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= 0.0 + assert value < 1.0 + } + } +} + +fn test_musl_f64() { + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= 0.0 + assert value < 1.0 + } + } +} + +fn test_musl_f32n() { + max := f32(357.0) + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= 0.0 + assert value < max + } + } +} + +fn test_musl_f64n() { + max := 1.52e6 + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= 0.0 + assert value < max + } + } +} + +fn test_musl_f32_in_range() { + min := f32(-24.0) + max := f32(125.0) + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= min + assert value < max + } + } +} + +fn test_musl_f64_in_range() { + min := -548.7 + max := 5015.2 + for seed in seeds { + mut rng := rand.MuslRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= min + assert value < max + } + } +} diff --git a/vlib/rand/pcg32.v b/vlib/rand/pcg32.v index 84c25c2c69..769ddb56fa 100644 --- a/vlib/rand/pcg32.v +++ b/vlib/rand/pcg32.v @@ -1,68 +1,234 @@ +// Copyright (c) 2019-2020 Alexander Medvednikov. All rights reserved. +// Use of this source code is governed by an MIT license +// that can be found in the LICENSE file. module rand -// Ported from http://www.pcg-random.org/download.html -// and https://github.com/imneme/pcg-c-basic/blob/master/pcg_basic.c -pub struct Pcg32 { + +// Ported from http://www.pcg-random.org/download.html, +// https://github.com/imneme/pcg-c-basic/blob/master/pcg_basic.c, and +// https://github.com/imneme/pcg-c-basic/blob/master/pcg_basic.h +pub struct PCG32RNG { mut: - state u64 - inc u64 + state u64 = u64(0x853c49e6748fea9b) ^ time_seed_64() + inc u64 = u64(0xda3e39cb94b95bdb) ^ time_seed_64() } -/** - * 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 = (initseq<> 32), u32(init_seq), u32(init_seq >> 32)]) return rng } -/** - * Pcg32.next - update the PRNG state and get back the next random number - * @return the generated pseudo random number -*/ - - -[inline] -pub fn (mut rng Pcg32) next() u32 { - oldstate := rng.state - rng.state = oldstate * (6364136223846793005) + rng.inc - xorshifted := u32(((oldstate>>u64(18)) ^ oldstate)>>u64(27)) - rot := u32(oldstate>>u64(59)) - return ((xorshifted>>rot) | (xorshifted<<((-rot) & u32(31)))) +pub fn (mut rng PCG32RNG) bounded_next(bound u32) u32 { + return rng.u32n(bound) } -/** - * 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 -*/ - +// rng.seed(seed_data) - seed the PCG32RNG with 4 u32 values. +// The first 2 represent the 64-bit initial state as [lower 32 bits, higher 32 bits] +// The last 2 represent the 64-bit stream/step of the PRNG. +pub fn (mut rng PCG32RNG) seed(seed_data []u32) { + if seed_data.len != 4 { + eprintln('PCG32RNG needs 4 u32s to be seeded. First two the initial state and the last two the stream/step. Both in little endian format: [lower, higher]') + exit(1) + } + init_state := u64(seed_data[0]) | (u64(seed_data[1]) << 32) + init_seq := u64(seed_data[2]) | (u64(seed_data[3]) << 32) + rng.state = u64(0) + rng.inc = (init_seq << u64(1)) | u64(1) + rng.u32() + rng.state += init_state + rng.u32() +} +// rng.u32() - return a pseudorandom 32 bit unsigned u32 [inline] -pub fn (mut rng Pcg32) bounded_next(bound u32) u32 { +pub fn (mut rng PCG32RNG) u32() u32 { + oldstate := rng.state + rng.state = oldstate * (6364136223846793005) + rng.inc + xorshifted := u32(((oldstate >> u64(18)) ^ oldstate) >> u64(27)) + rot := u32(oldstate >> u64(59)) + return ((xorshifted >> rot) | (xorshifted << ((-rot) & u32(31)))) +} + +// rng.u64() - return a pseudorandom 64 bit unsigned u64 +[inline] +pub fn (mut rng PCG32RNG) u64() u64 { + return u64(rng.u32()) | (u64(rng.u32()) << 32) +} + +// rn.u32n(max) - return a pseudorandom 32 bit unsigned u32 in [0, max) +[inline] +pub fn (mut rng PCG32RNG) u32n(max u32) u32 { + if max == 0 { + eprintln('max must be positive') + exit(1) + } // 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 := (-bound % bound) + // max, which we do by dropping output less than a threshold. + threshold := (-max % max) // Uniformity guarantees that loop below will terminate. In practice, it - // should usually terminate quickly; on average (assuming all bounds are + // should usually terminate quickly; on average (assuming all max's 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 + // one iteration. In practice, max's are typically small and only a // tiny amount of the range is eliminated. for { - r := rng.next() + r := rng.u32() if r >= threshold { - return (r % bound) + return (r % max) } } return u32(0) } + +// rn.u64n(max) - return a pseudorandom 64 bit unsigned u64 in [0, max) +[inline] +pub fn (mut rng PCG32RNG) u64n(max u64) u64 { + if max == 0 { + eprintln('max must be positive') + exit(1) + } + threshold := (-max % max) + for { + r := rng.u64() + if r >= threshold { + return (r % max) + } + } + return u64(0) +} + +// rn.u32_in_range(min, max) - return a pseudorandom 32 bit unsigned u32 in [min, max) +[inline] +pub fn (mut rng PCG32RNG) u32_in_range(min, max u64) u64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.u32n(max - min) +} + +// rn.u64_in_range(min, max) - return a pseudorandom 64 bit unsigned u64 in [min, max) +[inline] +pub fn (mut rng PCG32RNG) u64_in_range(min, max u64) u64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.u64n(max - min) +} + +// rng.int() - return a 32-bit signed (possibly negative) int +[inline] +pub fn (mut rng PCG32RNG) int() int { + return int(rng.u32()) +} + +// rng.i64() - return a 64-bit signed (possibly negative) i64 +[inline] +pub fn (mut rng PCG32RNG) i64() i64 { + return i64(rng.u64()) +} + +// rng.int31() - return a 31bit positive pseudorandom integer +[inline] +pub fn (mut rng PCG32RNG) int31() int { + return int(rng.u32() >> 1) +} + +// rng.int63() - return a 63bit positive pseudorandom integer +[inline] +pub fn (mut rng PCG32RNG) int63() i64 { + return i64(rng.u64() >> 1) +} + +// rng.intn(max) - return a 32bit positive int in [0, max) +[inline] +pub fn (mut rng PCG32RNG) intn(max int) int { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return int(rng.u32n(max)) +} + +// rng.i64n(max) - return a 64bit positive i64 in [0, max) +[inline] +pub fn (mut rng PCG32RNG) i64n(max i64) i64 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return i64(rng.u64n(max)) +} + +// rng.int_in_range(min, max) - return a 32bit positive int in [0, max) +[inline] +pub fn (mut rng PCG32RNG) int_in_range(min, max int) int { + if max <= min { + eprintln('max must be greater than min.') + exit(1) + } + return min + rng.intn(max - min) +} + +// rng.i64_in_range(min, max) - return a 64bit positive i64 in [0, max) +[inline] +pub fn (mut rng PCG32RNG) i64_in_range(min, max i64) i64 { + if max <= min { + eprintln('max must be greater than min.') + exit(1) + } + return min + rng.i64n(max - min) +} + +// rng.f32() returns a pseudorandom f32 value between 0.0 (inclusive) and 1.0 (exclusive) i.e [0, 1) +[inline] +pub fn (mut rng PCG32RNG) f32() f32 { + return f32(rng.u32()) / max_u32_as_f32 +} + +// rng.f64() returns a pseudorandom f64 value between 0.0 (inclusive) and 1.0 (exclusive) i.e [0, 1) +[inline] +pub fn (mut rng PCG32RNG) f64() f64 { + return f64(rng.u64()) / max_u64_as_f64 +} + +// rng.f32n() returns a pseudorandom f32 value in [0, max) +[inline] +pub fn (mut rng PCG32RNG) f32n(max f32) f32 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return rng.f32() * max +} + +// rng.f64n() returns a pseudorandom f64 value in [0, max) +[inline] +pub fn (mut rng PCG32RNG) f64n(max f64) f64 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return rng.f64() * max +} + +// rng.f32_in_range(min, max) returns a pseudorandom f32 that lies in [min, max) +[inline] +pub fn (mut rng PCG32RNG) f32_in_range(min, max f32) f32 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.f32n(max - min) +} + +// rng.i64_in_range(min, max) returns a pseudorandom i64 that lies in [min, max) +[inline] +pub fn (mut rng PCG32RNG) f64_in_range(min, max f64) f64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.f64n(max - min) +} diff --git a/vlib/rand/pcg32_test.v b/vlib/rand/pcg32_test.v index ad68ee771a..0476176069 100644 --- a/vlib/rand/pcg32_test.v +++ b/vlib/rand/pcg32_test.v @@ -1,36 +1,328 @@ - import rand -import time +import math -fn show_u32s(a []u32){ - mut res := []string{} - for x in a { - res << x.str() - } - print('[') - print(res.join(', ')) - println(']') -} -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)) +const ( + range_limit = 40 + value_count = 1000 + seeds = [[u32(42), 242, 267, 14195], [u32(256), 340, 1451, 1505]] +) + +const ( + sample_size = 1000 + stats_epsilon = 0.05 + inv_sqrt_12 = 1.0 / math.sqrt(12) +) + +fn gen_randoms(seed_data []u32, bound int) []u32 { + mut randoms := []u32{len: 20} + mut rng := rand.PCG32RNG{} + rng.seed(seed_data) + for i in 0 .. 20 { + randoms[i] = rng.u32n(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) + randoms1 := gen_randoms(rand.time_seed_array(4), 1000) + randoms2 := gen_randoms(rand.time_seed_array(4), 1000) assert randoms1.len == randoms2.len - show_u32s(randoms1) - show_u32s(randoms2) len := randoms1.len - for i in 0..len { + for i in 0 .. len { assert randoms1[i] == randoms2[i] } } + +// TODO: use the `in` syntax and remove this function +// after generics has been completely implemented +fn found(value u64, arr []u64) bool { + for item in arr { + if value == item { + return true + } + } + return false +} + +fn test_pcg32_variability() { + // If this test fails and if it is certainly not the implementation + // at fault, try changing the seed values. Repeated values are + // improbable but not impossible. + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + mut values := []u64{cap: value_count} + for i in 0 .. value_count { + value := rng.u64() + assert !found(value, values) + assert values.len == i + values << value + } + } +} + +fn check_uniformity_u64(rng rand.PCG32RNG, range u64) { + range_f64 := f64(range) + expected_mean := range_f64 / 2.0 + mut variance := 0.0 + for _ in 0 .. sample_size { + diff := f64(rng.u64n(range)) - expected_mean + variance += diff * diff + } + variance /= sample_size - 1 + sigma := math.sqrt(variance) + expected_sigma := range_f64 * inv_sqrt_12 + error := (sigma - expected_sigma) / expected_sigma + assert math.abs(error) < stats_epsilon +} + +fn test_pcg32_uniformity_u64() { + ranges := [14019545, 80240, 130] + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + for range in ranges { + check_uniformity_u64(rng, range) + } + } +} + +fn check_uniformity_f64(rng rand.PCG32RNG) { + expected_mean := 0.5 + mut variance := 0.0 + for _ in 0 .. sample_size { + diff := rng.f64() - expected_mean + variance += diff * diff + } + variance /= sample_size - 1 + sigma := math.sqrt(variance) + expected_sigma := inv_sqrt_12 + error := (sigma - expected_sigma) / expected_sigma + assert math.abs(error) < stats_epsilon +} + +fn test_pcg32_uniformity_f64() { + // The f64 version + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + check_uniformity_f64(rng) + } +} + +fn test_pcg32_u32n() { + max := 16384 + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u32n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_pcg32_u64n() { + max := u64(379091181005) + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u64n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_pcg32_u32_in_range() { + max := 484468466 + min := 316846 + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u32_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_pcg32_u64_in_range() { + max := u64(216468454685163) + min := u64(6848646868) + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u64_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_pcg32_int31() { + max_u31 := 0x7FFFFFFF + sign_mask := 0x80000000 + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int31() + assert value >= 0 + assert value <= max_u31 + // This statement ensures that the sign bit is zero + assert (value & sign_mask) == 0 + } + } +} + +fn test_pcg32_int63() { + max_u63 := i64(0x7FFFFFFFFFFFFFFF) + sign_mask := i64(0x8000000000000000) + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int63() + assert value >= 0 + assert value <= max_u63 + assert (value & sign_mask) == 0 + } + } +} + +fn test_pcg32_intn() { + max := 2525642 + for seed in seeds { + rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.intn(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_pcg32_i64n() { + max := i64(3246727724653636) + for seed in seeds { + rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.i64n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_pcg32_int_in_range() { + min := -4252 + max := 1034 + for seed in seeds { + rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_pcg32_i64_in_range() { + min := i64(-24095) + max := i64(324058) + for seed in seeds { + rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.i64_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_pcg32_f32() { + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= 0.0 + assert value < 1.0 + } + } +} + +fn test_pcg32_f64() { + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= 0.0 + assert value < 1.0 + } + } +} + +fn test_pcg32_f32n() { + max := f32(357.0) + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= 0.0 + assert value < max + } + } +} + +fn test_pcg32_f64n() { + max := 1.52e6 + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= 0.0 + assert value < max + } + } +} + +fn test_pcg32_f32_in_range() { + min := f32(-24.0) + max := f32(125.0) + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= min + assert value < max + } + } +} + +fn test_pcg32_f64_in_range() { + min := -548.7 + max := 5015.2 + for seed in seeds { + mut rng := rand.PCG32RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= min + assert value < max + } + } +} diff --git a/vlib/rand/rand.v b/vlib/rand/rand.v index 94dcdff0ff..cb39509267 100644 --- a/vlib/rand/rand.v +++ b/vlib/rand/rand.v @@ -3,8 +3,12 @@ // that can be found in the LICENSE file. module rand -fn C.rand() int - +// TODO: Remove these functions once done: +// 1. C.rand() +// 2. seed() +// 3. next() +// 4. rand_r() +// fn C.rand() int pub fn seed(s int) { C.srand(s) } @@ -23,32 +27,155 @@ pub fn rand_r(seed &int) int { return ns & 0x7fffffff } +const ( + default_rng = new_default({}) +) + +pub struct PRNGConfigStruct { + seed []u32 = time_seed_array(2) +} + +pub fn new_default(config PRNGConfigStruct) &WyRandRNG { + rng := &WyRandRNG{} + rng.seed(config.seed) + return rng +} + + +// u32() - returns a uniformly distributed pseudorandom 32-bit unsigned u32 +pub fn u32() u32 { + return default_rng.u32() +} + +// u64() - returns a uniformly distributed pseudorandom 64-bit unsigned u64 +pub fn u64() u64 { + return default_rng.u64() +} + +// u32n(max) - returns a uniformly distributed pseudorandom 32-bit unsigned u32 in [0, max) +pub fn u32n(max u32) u32 { + return default_rng.u32n(max) +} + +// u64n(max) - returns a uniformly distributed pseudorandom 64-bit unsigned u64 in [0, max) +pub fn u64n(max u64) u64 { + return default_rng.u64n(max) +} + +// u32n() - returns a uniformly distributed pseudorandom 32-bit unsigned u32 in [min, max) +pub fn u32_in_range(min, max u32) u32 { + return default_rng.u32_in_range(min, max) +} + +// u64_in_range(min, max) - returns a uniformly distributed pseudorandom 64-bit unsigned u64 in [min, max) +pub fn u64_in_range(min, max u64) u64 { + return default_rng.u64_in_range(min, max) +} + +// int() - returns a uniformly distributed pseudorandom 32-bit signed (possibly negative) int +pub fn int() int { + return default_rng.int() +} + +// intn(max) - returns a uniformly distributed pseudorandom 32-bit signed positive int in [0, max) +pub fn intn(max int) int { + return default_rng.intn(max) +} + +// int_in_range(min, max) - returns a uniformly distributed pseudorandom +// 32-bit signed int in [min, max) +pub fn int_in_range(min, max int) int { + return default_rng.int_in_range(min, max) +} + +// int31() - returns a uniformly distributed pseudorandom 31-bit signed positive int +pub fn int31() int { + return default_rng.int31() +} + +// i64() - returns a uniformly distributed pseudorandom 64-bit signed (possibly negative) i64 +pub fn i64() i64 { + return default_rng.i64() +} + +// i64n(max) - returns a uniformly distributed pseudorandom 64-bit signed positive i64 in [0, max) +pub fn i64n(max i64) i64 { + return default_rng.i64n(max) +} + +// i64_in_range(min, max) - returns a uniformly distributed pseudorandom +// 64-bit signed int in [min, max) +pub fn i64_in_range(min, max i64) i64 { + return default_rng.i64_in_range(min, max) +} + +// int63() - returns a uniformly distributed pseudorandom 63-bit signed positive int +pub fn int63() i64 { + return default_rng.int63() +} + +// f32() - returns a uniformly distributed 32-bit floating point in [0, 1) +pub fn f32() f32 { + return default_rng.f32() +} + +// f64() - returns a uniformly distributed 64-bit floating point in [0, 1) +pub fn f64() f64 { + return default_rng.f64() +} + +// f32n() - returns a uniformly distributed 32-bit floating point in [0, max) +pub fn f32n(max f32) f32 { + return default_rng.f32n(max) +} + +// f64n() - returns a uniformly distributed 64-bit floating point in [0, max) +pub fn f64n(max f64) f64 { + return default_rng.f64n(max) +} + +// f32_in_range(min, max) - returns a uniformly distributed 32-bit floating point in [min, max) +pub fn f32_in_range(min, max f32) f32 { + return default_rng.f32_in_range(min, max) +} + +// f64_in_range(min, max) - returns a uniformly distributed 64-bit floating point in [min, max) +pub fn f64_in_range(min, max f64) f64 { + return default_rng.f64_in_range(min, max) +} + // rand_f32 return a random f32 between 0 and max +[deprecated] pub fn rand_f32(max f32) f32 { return rand_uniform_f32() * max } // rand_f32 return a random f32 in range min and max +[deprecated] pub fn rand_f32_in_range(min, max f32) f32 { return min + rand_uniform_f32() * (max - min) } // rand_f64 return a random f64 between 0 (inclusive) and max (exclusive) +[deprecated] pub fn rand_f64(max f64) f64 { return rand_uniform_f64() * max } // rand_f64 return a random f64 in range min (inclusive) and max (exclusive) +[deprecated] pub fn rand_f64_in_range(min, max f64) f64 { return min + rand_uniform_f64() * (max - min) } // rand_uniform_f32 returns a uniformly distributed f32 in the range 0 (inclusive) and 1 (exclusive) +[deprecated] pub fn rand_uniform_f32() f32 { return f32(C.rand()) / f32(C.RAND_MAX) } // rand_uniform_f64 returns a uniformly distributed f64 in the range 0 (inclusive) and 1 (exclusive) +[deprecated] pub fn rand_uniform_f64() f64 { return f64(C.rand()) / f64(C.RAND_MAX) -} \ No newline at end of file +} diff --git a/vlib/rand/random_numbers_test.v b/vlib/rand/random_numbers_test.v index d7a758f3bb..605b8a0dea 100644 --- a/vlib/rand/random_numbers_test.v +++ b/vlib/rand/random_numbers_test.v @@ -1,4 +1,5 @@ import rand +import math const ( rnd_count = 40 @@ -53,7 +54,7 @@ fn assert_randoms_equal(r1, r2 []int) { } } -fn test_rand_f32() { +fn test_rand_f32_old() { for seed in seeds { rand.seed(seed) for _ in 0 .. rnd_count { @@ -64,7 +65,7 @@ fn test_rand_f32() { } } -fn test_rand_f32_in_range() { +fn test_rand_f32_in_range_old() { for seed in seeds { rand.seed(seed) for _ in 0 .. rnd_count { @@ -75,7 +76,7 @@ fn test_rand_f32_in_range() { } } -fn test_rand_f64() { +fn test_rand_f64_old() { for seed in seeds { rand.seed(seed) for _ in 0 .. rnd_count { @@ -86,7 +87,7 @@ fn test_rand_f64() { } } -fn test_rand_f64_in_range() { +fn test_rand_f64_in_range_old() { for seed in seeds { rand.seed(seed) for _ in 0 .. rnd_count { @@ -118,3 +119,145 @@ fn test_rand_uniform_f64() { } } } + +fn test_rand_u32n() { + max := u32(4287502) + for _ in 0 .. rnd_count { + assert rand.u32n(max) < max + } +} + +fn test_rand_u64n() { + max := u64(23442353534587502) + for _ in 0 .. rnd_count { + assert rand.u64n(max) < max + } +} + +fn test_rand_u32_in_range() { + min := u32(5256) + max := u32(4287502) + for _ in 0 .. rnd_count { + value := rand.u32_in_range(min, max) + assert value >= min + assert value < max + } +} + +fn test_rand_u64_in_range() { + min := u64(4265266246) + max := u64(23442353534587502) + for _ in 0 .. rnd_count { + value := rand.u64_in_range(min, max) + assert value >= min + assert value < max + } +} + +fn test_rand_intn() { + max := 720948723 + for _ in 0 .. rnd_count { + value := rand.intn(max) + assert value >= 0 + assert value < max + } +} + +fn test_rand_i64n() { + max := i64(209487239094) + for _ in 0 .. rnd_count { + value := rand.i64n(max) + assert value >= 0 + assert value < max + } +} + +fn test_rand_int_in_range() { + min := -34058 + max := -10542 + for _ in 0 .. rnd_count { + value := rand.int_in_range(min, max) + assert value >= min + assert value < max + } +} + +fn test_rand_i64_in_range() { + min := i64(-5026245) + max := i64(209487239094) + for _ in 0 .. rnd_count { + value := rand.i64_in_range(min, max) + assert value >= min + assert value < max + } +} + +fn test_rand_int31() { + for _ in 0 .. rnd_count { + value := rand.int31() + assert value >= 0 + assert value <= math.max_i32 + } +} + +fn test_rand_int63() { + for _ in 0 .. rnd_count { + value := rand.int63() + assert value >= 0 + assert value <= math.max_i64 + } +} + +fn test_rand_f32() { + for _ in 0 .. rnd_count { + value := rand.f32() + assert value >= 0.0 + assert value < 1.0 + } +} + +fn test_rand_f64() { + for _ in 0 .. rnd_count { + value := rand.f64() + assert value >= 0.0 + assert value < 1.0 + } +} + +fn test_rand_f32n() { + max := f32(34.52) + for _ in 0 .. rnd_count { + value := rand.f32n(max) + assert value >= 0.0 + assert value < max + } +} + +fn test_rand_f64n() { + max := 3495.2 + for _ in 0 .. rnd_count { + value := rand.f64n(max) + assert value >= 0.0 + assert value < max + } +} + +fn test_rand_f32_in_range() { + min := f32(-10.4) + max := f32(43.2) + for _ in 0 .. rnd_count { + value := rand.f32_in_range(min, max) + assert value >= min + assert value < max + } +} + +fn test_rand_f64_in_range() { + min := -10980.4 + max := -2.0 + for _ in 0 .. rnd_count { + value := rand.f64_in_range(min, max) + assert value >= min + assert value < max + } +} \ No newline at end of file diff --git a/vlib/rand/splitmix64.v b/vlib/rand/splitmix64.v index c1c26099a6..5ad00cb577 100644 --- a/vlib/rand/splitmix64.v +++ b/vlib/rand/splitmix64.v @@ -1,52 +1,222 @@ +// Copyright (c) 2019-2020 Alexander Medvednikov. All rights reserved. +// Use of this source code is governed by an MIT license +// that can be found in the LICENSE file. module rand + // Ported from http://xoshiro.di.unimi.it/splitmix64.c -struct Splitmix64 { +pub struct SplitMix64RNG { mut: - state u64 + state u64 = time_seed_64() + has_extra bool = false + extra u32 } -/** - * 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} +// rng.seed(seed_data) sets the seed of the accepting SplitMix64RNG to the given data +// in little-endian format (i.e. lower 32 bits are in [0] and higher 32 bits in [1]). +pub fn (mut rng SplitMix64RNG) seed(seed_data []u32) { + if seed_data.len != 2 { + eprintln('SplitMix64RNG needs 2 32-bit unsigned integers as the seed.') + exit(1) + } + rng.state = seed_data[0] | (u64(seed_data[1]) << 32) + rng.has_extra = false } -/** - * Splitmix64.next - update the PRNG state and get back the next random number - * @return the generated pseudo random number -*/ - - +// rng.u32() updates the PRNG state and returns the next pseudorandom u32 [inline] -pub fn (mut rng Splitmix64) next() u64 { +pub fn (mut rng SplitMix64RNG) u32() u32 { + if rng.has_extra { + rng.has_extra = false + return rng.extra + } + full_value := rng.u64() + lower := u32(full_value & lower_mask) + upper := u32(full_value >> 32) + rng.extra = upper + rng.has_extra = true + return lower +} + +// rng.u64() updates the PRNG state and returns the next pseudorandom u64 +[inline] +pub fn (mut rng SplitMix64RNG) u64() u64 { rng.state += (0x9e3779b97f4a7c15) mut z := rng.state - z = (z ^ ((z>>u64(30)))) * (0xbf58476d1ce4e5b9) - z = (z ^ ((z>>u64(27)))) * (0x94d049bb133111eb) - return z ^ (z>>(31)) + z = (z ^ ((z >> u64(30)))) * (0xbf58476d1ce4e5b9) + z = (z ^ ((z >> u64(27)))) * (0x94d049bb133111eb) + return z ^ (z >> (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 -*/ - - +// rng.u32n(bound) returns a pseudorandom u32 less than the bound [inline] -pub fn (mut rng Splitmix64) bounded_next(bound u64) u64 { +pub fn (mut rng SplitMix64RNG) u32n(bound u32) u32 { + // This function is kept similar to the u64 version + if bound == 0 { + eprintln('max must be non-zero') + exit(1) + } threshold := -bound % bound for { - r := rng.next() + r := rng.u32() + if r >= threshold { + return r % bound + } + } + return u32(0) +} + +// rng.u64n(bound) returns a pseudorandom u64 less than the bound +[inline] +pub fn (mut rng SplitMix64RNG) u64n(bound u64) u64 { + // See pcg32.v for explanation of comment. This algorithm + // existed before the refactoring. + if bound == 0 { + eprintln('max must be non-zero') + exit(1) + } + threshold := -bound % bound + for { + r := rng.u64() if r >= threshold { return r % bound } } return u64(0) } + +// rng.u32n(min, max) returns a pseudorandom u32 value that is guaranteed to be in [min, max) +[inline] +pub fn (mut rng SplitMix64RNG) u32_in_range(min, max u32) u32 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.u32n(max - min) +} + +// rng.u64n(min, max) returns a pseudorandom u64 value that is guaranteed to be in [min, max) +[inline] +pub fn (mut rng SplitMix64RNG) u64_in_range(min, max u64) u64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.u64n(max - min) +} + +// rng.int() returns a pseudorandom 32-bit int (which may be negative) +[inline] +pub fn (mut rng SplitMix64RNG) int() int { + return int(rng.u32()) +} + +// rng.i64() returns a pseudorandom 64-bit i64 (which may be negative) +[inline] +pub fn (mut rng SplitMix64RNG) i64() i64 { + return i64(rng.u64()) +} + +// rng.int31() returns a pseudorandom 31-bit int which is non-negative +[inline] +pub fn (mut rng SplitMix64RNG) int31() int { + return int(rng.u32() & u31_mask) // Set the 32nd bit to 0. +} + +// rng.int63() returns a pseudorandom 63-bit int which is non-negative +[inline] +pub fn (mut rng SplitMix64RNG) int63() i64 { + return i64(rng.u64() & u63_mask) // Set the 64th bit to 0. +} + +// rng.intn(max) returns a pseudorandom int that lies in [0, max) +[inline] +pub fn (mut rng SplitMix64RNG) intn(max int) int { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return int(rng.u32n(max)) +} + +// rng.i64n(max) returns a pseudorandom int that lies in [0, max) +[inline] +pub fn (mut rng SplitMix64RNG) i64n(max i64) i64 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return i64(rng.u64n(max)) +} + +// rng.int_in_range(min, max) returns a pseudorandom int that lies in [min, max) +[inline] +pub fn (mut rng SplitMix64RNG) int_in_range(min, max int) int { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + // This supports negative ranges like [-10, -5) because the difference is positive + return min + rng.intn(max - min) +} + +// rng.i64_in_range(min, max) returns a pseudorandom i64 that lies in [min, max) +[inline] +pub fn (mut rng SplitMix64RNG) i64_in_range(min, max i64) i64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.i64n(max - min) +} + +// rng.f32() returns a pseudorandom f32 value between 0.0 (inclusive) and 1.0 (exclusive) i.e [0, 1) +[inline] +pub fn (mut rng SplitMix64RNG) f32() f32 { + return f32(rng.u32()) / max_u32_as_f32 +} + +// rng.f64() returns a pseudorandom f64 value between 0.0 (inclusive) and 1.0 (exclusive) i.e [0, 1) +[inline] +pub fn (mut rng SplitMix64RNG) f64() f64 { + return f64(rng.u64()) / max_u64_as_f64 +} + +// rng.f32n() returns a pseudorandom f32 value in [0, max) +[inline] +pub fn (mut rng SplitMix64RNG) f32n(max f32) f32 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return rng.f32() * max +} + +// rng.f64n() returns a pseudorandom f64 value in [0, max) +[inline] +pub fn (mut rng SplitMix64RNG) f64n(max f64) f64 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return rng.f64() * max +} + +// rng.f32_in_range(min, max) returns a pseudorandom f32 that lies in [min, max) +[inline] +pub fn (mut rng SplitMix64RNG) f32_in_range(min, max f32) f32 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.f32n(max - min) +} + +// rng.i64_in_range(min, max) returns a pseudorandom i64 that lies in [min, max) +[inline] +pub fn (mut rng SplitMix64RNG) f64_in_range(min, max f64) f64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.f64n(max - min) +} diff --git a/vlib/rand/splitmix64_test.v b/vlib/rand/splitmix64_test.v index 6f3adebd58..0099a1e403 100644 --- a/vlib/rand/splitmix64_test.v +++ b/vlib/rand/splitmix64_test.v @@ -1,36 +1,330 @@ - import rand -import time +import math -fn show_u64s(a []u64){ - mut res := []string{} - for x in a { - res << x.str() - } - print('[') - print(res.join(', ')) - println(']') -} +const ( + range_limit = 40 + value_count = 1000 + seeds = [[u32(42), 0], [u32(256), 0]] +) -fn gen_randoms(seed i64, bound int) []u64 { +const ( + sample_size = 1000 + stats_epsilon = 0.05 + inv_sqrt_12 = 1.0 / math.sqrt(12) +) + +fn gen_randoms(seed_data []u32, bound int) []u64 { + bound_u64 := u64(bound) 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)) + mut rnd := rand.SplitMix64RNG{} + rnd.seed(seed_data) + for i in 0 .. 20 { + randoms[i] = rnd.u64n(bound_u64) } return randoms } fn test_splitmix64_reproducibility() { - t := time.ticks() - println('t: $t') - randoms1 := gen_randoms(t, 1000) - randoms2 := gen_randoms(t, 1000) + seed_data := rand.time_seed_array(2) + randoms1 := gen_randoms(seed_data, 1000) + randoms2 := gen_randoms(seed_data, 1000) assert randoms1.len == randoms2.len - show_u64s( randoms1 ) - show_u64s( randoms2 ) len := randoms1.len - for i in 0..len { + for i in 0 .. len { assert randoms1[i] == randoms2[i] } } + +// TODO: use the `in` syntax and remove this function +// after generics has been completely implemented +fn found(value u64, arr []u64) bool { + for item in arr { + if value == item { + return true + } + } + return false +} + +fn test_splitmix64_variability() { + // If this test fails and if it is certainly not the implementation + // at fault, try changing the seed values. Repeated values are + // improbable but not impossible. + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + mut values := []u64{cap: value_count} + for i in 0 .. value_count { + value := rng.u64() + assert !found(value, values) + assert values.len == i + values << value + } + } +} + +fn check_uniformity_u64(rng rand.SplitMix64RNG, range u64) { + range_f64 := f64(range) + expected_mean := range_f64 / 2.0 + mut variance := 0.0 + for _ in 0 .. sample_size { + diff := f64(rng.u64n(range)) - expected_mean + variance += diff * diff + } + variance /= sample_size - 1 + sigma := math.sqrt(variance) + expected_sigma := range_f64 * inv_sqrt_12 + error := (sigma - expected_sigma) / expected_sigma + assert math.abs(error) < stats_epsilon +} + +fn test_splitmix64_uniformity_u64() { + ranges := [14019545, 80240, 130] + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + for range in ranges { + check_uniformity_u64(rng, range) + } + } +} + +fn check_uniformity_f64(rng rand.SplitMix64RNG) { + expected_mean := 0.5 + mut variance := 0.0 + for _ in 0 .. sample_size { + diff := rng.f64() - expected_mean + variance += diff * diff + } + variance /= sample_size - 1 + sigma := math.sqrt(variance) + expected_sigma := inv_sqrt_12 + error := (sigma - expected_sigma) / expected_sigma + assert math.abs(error) < stats_epsilon +} + +fn test_splitmix64_uniformity_f64() { + // The f64 version + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + check_uniformity_f64(rng) + } +} + +fn test_splitmix64_u32n() { + max := 16384 + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u32n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_splitmix64_u64n() { + max := u64(379091181005) + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u64n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_splitmix64_u32_in_range() { + max := 484468466 + min := 316846 + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u32_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_splitmix64_u64_in_range() { + max := u64(216468454685163) + min := u64(6848646868) + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u64_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_splitmix64_int31() { + max_u31 := 0x7FFFFFFF + sign_mask := 0x80000000 + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int31() + assert value >= 0 + assert value <= max_u31 + // This statement ensures that the sign bit is zero + assert (value & sign_mask) == 0 + } + } +} + +fn test_splitmix64_int63() { + max_u63 := i64(0x7FFFFFFFFFFFFFFF) + sign_mask := i64(0x8000000000000000) + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int63() + assert value >= 0 + assert value <= max_u63 + assert (value & sign_mask) == 0 + } + } +} + +fn test_splimix64_intn() { + max := 2525642 + for seed in seeds { + rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.intn(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_splimix64_i64n() { + max := i64(3246727724653636) + for seed in seeds { + rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.i64n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_splimix64_int_in_range() { + min := -4252 + max := 230549862 + for seed in seeds { + rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_splimix64_i64_in_range() { + min := i64(-24095) + max := i64(324058) + for seed in seeds { + rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.i64_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_splitmix64_f32() { + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= 0.0 + assert value < 1.0 + } + } +} + +fn test_splitmix64_f64() { + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= 0.0 + assert value < 1.0 + } + } +} + +fn test_splitmix64_f32n() { + max := f32(357.0) + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= 0.0 + assert value < max + } + } +} + +fn test_splitmix64_f64n() { + max := 1.52e6 + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= 0.0 + assert value < max + } + } +} + +fn test_splitmix64_f32_in_range() { + min := f32(-24.0) + max := f32(125.0) + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= min + assert value < max + } + } +} + +fn test_splitmix64_f64_in_range() { + min := -548.7 + max := 5015.2 + for seed in seeds { + mut rng := rand.SplitMix64RNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= min + assert value < max + } + } +} diff --git a/vlib/rand/system_rng.c.v b/vlib/rand/system_rng.c.v new file mode 100644 index 0000000000..f957458982 --- /dev/null +++ b/vlib/rand/system_rng.c.v @@ -0,0 +1,291 @@ +// Copyright (c) 2019-2020 Alexander Medvednikov. All rights reserved. +// Use of this source code is governed by an MIT license +// that can be found in the LICENSE file. +module rand + +import math.bits + +// Implementation note: +// ==================== +// C.rand() is okay to use within its defined range of C.RAND_MAX. +// (See: https://web.archive.org/web/20180801210127/http://eternallyconfuzzled.com/arts/jsw_art_rand.aspx) +// The problem is, this value varies with the libc implementation. On windows, +// for example, RAND_MAX is usually a measly 32767, whereas on (newer) linux it's generaly +// 2147483647. The repetition period also varies wildly. In order to provide more entropy +// without altering the underlying algorithm too much, this implementation simply +// requests for more random bits until the necessary width for the integers is achieved. +const ( + rand_limit = u64(C.RAND_MAX) + rand_bitsize = bits.len_64(rand_limit) + u32_iter_count = calculate_iterations_for(32) + u64_iter_count = calculate_iterations_for(64) +) + +fn calculate_iterations_for(bits int) int { + base := bits / rand_bitsize + extra := if bits % rand_bitsize == 0 { 0 } else { 1 } + return base + extra +} + +// Size constants to avoid importing the entire math module +const ( + max_u32 = 0xFFFFFFFF + max_u64 = 0xFFFFFFFFFFFFFFFF + max_u32_as_f32 = f32(max_u32) + max_u64_as_f64 = f64(max_u64) +) + +// Masks for fast modular division +const ( + u31_mask = u32(0x7FFFFFFF) + u63_mask = u64(0x7FFFFFFFFFFFFFFF) +) + +// C.rand returns a pseudorandom integer from 0 (inclusive) to C.RAND_MAX (exclusive) +fn C.rand() int + +// C.srand seeds the internal PRNG with the given int seed. +// fn C.srand(seed int) +// SysRNG is the PRNG provided by default in the libc implementiation that V uses. +pub struct SysRNG { +mut: + seed u32 = time_seed_32() +} + +// r.seed() sets the seed of the accepting SysRNG to the given data. +pub fn (mut r SysRNG) seed(seed_data []u32) { + if seed_data.len != 1 { + eprintln('SysRNG needs one 32-bit unsigned integer as the seed.') + exit(1) + } + r.seed = seed_data[0] + C.srand(int(r.seed)) +} + +// r.default_rand() exposes the default behavior of the system's RNG +// (equivalent to calling C.rand()). Recommended for testing/comparison +// b/w V and other languages using libc and not for regular use. +// This is also a one-off feature of SysRNG, similar to the global seed +// situation. Other generators will not have this. +[inline] +pub fn (r SysRNG) default_rand() int { + return C.rand() +} + +// r.u32() returns a pseudorandom u32 value less than 2^32 +[inline] +pub fn (r SysRNG) u32() u32 { + mut result := u32(C.rand()) + for i in 1 .. u32_iter_count { + result = result ^ (u32(C.rand()) << (rand_bitsize * i)) + } + return result +} + +// r.u64() returns a pseudorandom u64 value less than 2^64 +[inline] +pub fn (r SysRNG) u64() u64 { + mut result := u64(C.rand()) + for i in 1 .. u64_iter_count { + result = result ^ (u64(C.rand()) << (rand_bitsize * i)) + } + return result +} + +// r.u32n(max) returns a pseudorandom u32 value that is guaranteed to be less than max +[inline] +pub fn (r SysRNG) u32n(max u32) u32 { + if max == 0 { + eprintln('max must be positive integer') + exit(1) + } + // Owing to the pigeon-hole principle, we can't simply do + // val := rng.u32() % max. + // It'll wreck the properties of the distribution unless + // max evenly divides 2^32. So we divide evenly to + // the closest power of two. Then we loop until we find + // an int in the required range + bit_len := bits.len_32(max) + if bit_len == 32 { + for { + value := r.u32() + if value < max { + return value + } + } + } else { + mask := (u32(1) << (bit_len + 1)) - 1 + for { + value := r.u32() & mask + if value < max { + return value + } + } + } + return u32(0) +} + +// r.u64n(max) returns a pseudorandom u64 value that is guaranteed to be less than max +[inline] +pub fn (r SysRNG) u64n(max u64) u64 { + if max == 0 { + eprintln('max must be positive integer') + exit(1) + } + // Similar procedure for u64s + bit_len := bits.len_64(max) + if bit_len == 64 { + for { + value := r.u64() + if value < max { + return value + } + } + } else { + mask := (u64(1) << (bit_len + 1)) - 1 + for { + value := r.u64() & mask + if value < max { + return value + } + } + } + return u64(0) +} + +// r.u32n(min, max) returns a pseudorandom u32 value that is guaranteed to be in [min, max) +[inline] +pub fn (r SysRNG) u32_in_range(min, max u32) u32 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + r.u32n(max - min) +} + +// r.u64n(min, max) returns a pseudorandom u64 value that is guaranteed to be in [min, max) +[inline] +pub fn (r SysRNG) u64_in_range(min, max u64) u64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + r.u64n(max - min) +} + +// r.int() returns a pseudorandom 32-bit int (which may be negative) +[inline] +pub fn (r SysRNG) int() int { + return int(r.u32()) +} + +// r.i64() returns a pseudorandom 64-bit i64 (which may be negative) +[inline] +pub fn (r SysRNG) i64() i64 { + return i64(r.u64()) +} + +// r.int31() returns a pseudorandom 31-bit int which is non-negative +[inline] +pub fn (r SysRNG) int31() int { + return int(r.u32() & u31_mask) // Set the 32nd bit to 0. +} + +// r.int63() returns a pseudorandom 63-bit int which is non-negative +[inline] +pub fn (r SysRNG) int63() i64 { + return i64(r.u64() & u63_mask) // Set the 64th bit to 0. +} + +// r.intn(max) returns a pseudorandom int that lies in [0, max) +[inline] +pub fn (r SysRNG) intn(max int) int { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return int(r.u32n(max)) +} + +// r.i64n(max) returns a pseudorandom i64 that lies in [0, max) +[inline] +pub fn (r SysRNG) i64n(max i64) i64 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return i64(r.u64n(max)) +} + +// r.int_in_range(min, max) returns a pseudorandom int that lies in [min, max) +[inline] +pub fn (r SysRNG) int_in_range(min, max int) int { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + // This supports negative ranges like [-10, -5) because the difference is positive + return min + r.intn(max - min) +} + +// r.i64_in_range(min, max) returns a pseudorandom i64 that lies in [min, max) +[inline] +pub fn (r SysRNG) i64_in_range(min, max i64) i64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + r.i64n(max - min) +} + +// r.f32() returns a pseudorandom f32 value between 0.0 (inclusive) and 1.0 (exclusive) i.e [0, 1) +[inline] +pub fn (r SysRNG) f32() f32 { + return f32(r.u32()) / max_u32_as_f32 +} + +// r.f64() returns a pseudorandom f64 value between 0.0 (inclusive) and 1.0 (exclusive) i.e [0, 1) +[inline] +pub fn (r SysRNG) f64() f64 { + return f64(r.u64()) / max_u64_as_f64 +} + +// r.f32n() returns a pseudorandom f32 value in [0, max) +[inline] +pub fn (r SysRNG) f32n(max f32) f32 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return r.f32() * max +} + +// r.f64n() returns a pseudorandom f64 value in [0, max) +[inline] +pub fn (r SysRNG) f64n(max f64) f64 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return r.f64() * max +} + +// r.f32_in_range(min, max) returns a pseudorandom f32 that lies in [min, max) +[inline] +pub fn (r SysRNG) f32_in_range(min, max f32) f32 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + r.f32n(max - min) +} + +// r.i64_in_range(min, max) returns a pseudorandom i64 that lies in [min, max) +[inline] +pub fn (r SysRNG) f64_in_range(min, max f64) f64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + r.f64n(max - min) +} diff --git a/vlib/rand/system_rng.js.v b/vlib/rand/system_rng.js.v new file mode 100644 index 0000000000..223ab3a0ac --- /dev/null +++ b/vlib/rand/system_rng.js.v @@ -0,0 +1,15 @@ +// Copyright (c) 2019-2020 Alexander Medvednikov. All rights reserved. +// Use of this source code is governed by an MIT license +// that can be found in the LICENSE file. +module rand + +// Until there's a portable, JS has a seeded way to produce random numbers +// and not just Math.random(), use any of the existing implementations +// as the System's RNG +type SysRNG WyRandRNG + +// In the JS version, we simply return the same int as is normally generated. +[inline] +pub fn (r SysRNG) default_rand() int { + return r.int() +} \ No newline at end of file diff --git a/vlib/rand/system_rng_test.v b/vlib/rand/system_rng_test.v new file mode 100644 index 0000000000..bc4077fa5d --- /dev/null +++ b/vlib/rand/system_rng_test.v @@ -0,0 +1,354 @@ +import rand +import math + +const ( + range_limit = 40 + value_count = 1000 + seeds = [u32(42), 256] +) + +const ( + sample_size = 1000 + stats_epsilon = 0.05 + inv_sqrt_12 = 1.0 / math.sqrt(12) +) + +fn get_n_randoms(n int, r rand.SysRNG) []int { + mut ints := []int{cap: n} + for _ in 0 .. n { + ints << r.int() + } + return ints +} + +fn test_sys_rng_reproducibility() { + // Note that C.srand() sets the seed globally. + // So the order of seeding matters. It is recommended + // to obtain all necessary data first, then set the + // seed for another batch of data. + for seed in seeds { + seed_data := [seed] + r1 := rand.SysRNG{} + r2 := rand.SysRNG{} + r1.seed(seed_data) + ints1 := get_n_randoms(value_count, r1) + r2.seed(seed_data) + ints2 := get_n_randoms(value_count, r2) + assert ints1 == ints2 + } +} + +// TODO: use the `in` syntax and remove this function +// after generics has been completely implemented +fn found(value u64, arr []u64) bool { + for item in arr { + if value == item { + return true + } + } + return false +} + +fn test_sys_rng_variability() { + // If this test fails and if it is certainly not the implementation + // at fault, try changing the seed values. Repeated values are + // improbable but not impossible. + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + mut values := []u64{cap: value_count} + for i in 0 .. value_count { + value := rng.u64() + assert !found(value, values) + assert values.len == i + values << value + } + } +} + +fn check_uniformity_u64(rng rand.SysRNG, range u64) { + range_f64 := f64(range) + expected_mean := range_f64 / 2.0 + mut variance := 0.0 + for _ in 0 .. sample_size { + diff := f64(rng.u64n(range)) - expected_mean + variance += diff * diff + } + variance /= sample_size - 1 + sigma := math.sqrt(variance) + expected_sigma := range_f64 * inv_sqrt_12 + error := (sigma - expected_sigma) / expected_sigma + assert math.abs(error) < stats_epsilon +} + +fn test_sys_rng_uniformity_u64() { + // This assumes that C.rand() produces uniform results to begin with. + // If the failure persists, report an issue on GitHub + ranges := [14019545, 80240, 130] + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for range in ranges { + check_uniformity_u64(rng, range) + } + } +} + +fn check_uniformity_f64(rng rand.SysRNG) { + expected_mean := 0.5 + mut variance := 0.0 + for _ in 0 .. sample_size { + diff := rng.f64() - expected_mean + variance += diff * diff + } + variance /= sample_size - 1 + sigma := math.sqrt(variance) + expected_sigma := inv_sqrt_12 + error := (sigma - expected_sigma) / expected_sigma + assert math.abs(error) < stats_epsilon +} + +fn test_sys_rng_uniformity_f64() { + // The f64 version + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + check_uniformity_f64(rng) + } +} + +fn test_sys_rng_u32n() { + max := 16384 + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.u32n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_sys_rng_u64n() { + max := u64(379091181005) + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.u64n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_sys_rng_u32_in_range() { + max := 484468466 + min := 316846 + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.u32_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_sys_rng_u64_in_range() { + max := u64(216468454685163) + min := u64(6848646868) + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.u64_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_sys_rng_intn() { + max := 2525642 + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.intn(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_sys_rng_i64n() { + max := i64(3246727724653636) + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.i64n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_sys_rng_int_in_range() { + min := -4252 + max := 23054962 + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.int_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_sys_rng_i64_in_range() { + min := i64(-24095) + max := i64(324058) + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.i64_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_sys_rng_int31() { + max_u31 := 0x7FFFFFFF + sign_mask := 0x80000000 + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.int31() + assert value >= 0 + assert value <= max_u31 + // This statement ensures that the sign bit is zero + assert (value & sign_mask) == 0 + } + } +} + +fn test_sys_rng_int63() { + max_u63 := i64(0x7FFFFFFFFFFFFFFF) + sign_mask := i64(0x8000000000000000) + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.int63() + assert value >= 0 + assert value <= max_u63 + assert (value & sign_mask) == 0 + } + } +} + +fn test_sys_rng_f32() { + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= 0.0 + assert value < 1.0 + } + } +} + +fn test_sys_rng_f64() { + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= 0.0 + assert value < 1.0 + } + } +} + +fn test_sys_rng_f32n() { + max := f32(357.0) + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= 0.0 + assert value < max + } + } +} + +fn test_sys_rng_f64n() { + max := 1.52e6 + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= 0.0 + assert value < max + } + } +} + +fn test_sys_rng_f32_in_range() { + min := f32(-24.0) + max := f32(125.0) + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= min + assert value < max + } + } +} + +fn test_sys_rng_f64_in_range() { + min := -548.7 + max := 5015.2 + for seed in seeds { + seed_data := [seed] + rng := rand.SysRNG{} + rng.seed(seed_data) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= min + assert value < max + } + } +} diff --git a/vlib/rand/util.v b/vlib/rand/util.v new file mode 100644 index 0000000000..cc021181cd --- /dev/null +++ b/vlib/rand/util.v @@ -0,0 +1,42 @@ +// Copyright (c) 2019-2020 Alexander Medvednikov. All rights reserved. +// Use of this source code is governed by an MIT license +// that can be found in the LICENSE file. +module rand + +import time + +// Commonly used constants across RNGs +const ( + lower_mask = u64(0x00000000ffffffff) +) + +// Constants taken from Numerical Recipes +[inline] +fn nr_next(prev u32) u32 { + return prev * 1664525 + 1013904223 +} + +// utility function that return the required number of u32s generated from system time +[inline] +pub fn time_seed_array(count int) []u32 { + mut seed := u32(time.now().unix_time()) + mut seed_data := []u32{cap: count} + for _ in 0 .. count { + seed = nr_next(seed) + seed_data << nr_next(seed) + } + return seed_data +} + +[inline] +fn time_seed_32() u32 { + return time_seed_array(1)[0] +} + +[inline] +fn time_seed_64() u64 { + seed_data := time_seed_array(2) + lower := u64(seed_data[0]) + upper := u64(seed_data[1]) + return lower | (upper << 32) +} diff --git a/vlib/rand/wyrand.v b/vlib/rand/wyrand.v new file mode 100644 index 0000000000..3ec2263211 --- /dev/null +++ b/vlib/rand/wyrand.v @@ -0,0 +1,251 @@ +// Copyright (c) 2019-2020 Alexander Medvednikov. All rights reserved. +// Use of this source code is governed by an MIT license +// that can be found in the LICENSE file. +module rand + +import math.bits +import hash.wyhash + +// Redefinition of some constants that we will need for pseudorandom number generation +const ( + wyp0 = u64(0xa0761d6478bd642f) + wyp1 = u64(0xe7037ed1a0b428db) +) + +// RNG based on the WyHash hashing algorithm +pub struct WyRandRNG { +mut: + state u64 = time_seed_64() + has_extra bool = false + extra u32 +} + +// seed() - Set the seed, needs only two u32s in little endian format as [lower, higher] +pub fn (mut rng WyRandRNG) seed(seed_data []u32) { + if seed_data.len != 2 { + eprintln('WyRandRNG needs 2 32-bit unsigned integers as the seed.') + exit(1) + } + rng.state = seed_data[0] | (u64(seed_data[1]) << 32) + rng.has_extra = false +} + + +// rng.u32() updates the PRNG state and returns the next pseudorandom u32 +[inline] +pub fn (mut rng WyRandRNG) u32() u32 { + if rng.has_extra { + rng.has_extra = false + return rng.extra + } + full_value := rng.u64() + lower := u32(full_value & lower_mask) + upper := u32(full_value >> 32) + rng.extra = upper + rng.has_extra = true + return lower +} + +// rng.u64() updates the PRNG state and returns the next pseudorandom u64 +[inline] +pub fn (mut rng WyRandRNG) u64() u64 { + unsafe { + mut seed1 := rng.state + seed1 += wyp0 + rng.state = seed1 + return wyhash.wymum(seed1 ^ wyp1, seed1) + } + return 0 +} + +// rng.u32n(max) returns a pseudorandom u32 less than the max +[inline] +pub fn (mut rng WyRandRNG) u32n(max u32) u32 { + if max == 0 { + eprintln('max must be positive integer') + exit(1) + } + // Check SysRNG in system_rng.c.v for explanation + bit_len := bits.len_32(max) + if bit_len == 32 { + for { + value := rng.u32() + if value < max { + return value + } + } + } else { + mask := (u32(1) << (bit_len + 1)) - 1 + for { + value := rng.u32() & mask + if value < max { + return value + } + } + } + return u32(0) +} + +// rng.u64n(max) returns a pseudorandom u64 less than the max +[inline] +pub fn (mut rng WyRandRNG) u64n(max u64) u64 { + if max == 0 { + eprintln('max must be positive integer') + exit(1) + } + bit_len := bits.len_64(max) + if bit_len == 64 { + for { + value := rng.u64() + if value < max { + return value + } + } + } else { + mask := (u64(1) << (bit_len + 1)) - 1 + for { + value := rng.u64() & mask + if value < max { + return value + } + } + } + return u64(0) +} + +// rng.u32n(min, max) returns a pseudorandom u32 value that is guaranteed to be in [min, max) +[inline] +pub fn (mut rng WyRandRNG) u32_in_range(min, max u32) u32 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.u32n(max - min) +} + +// rng.u64n(min, max) returns a pseudorandom u64 value that is guaranteed to be in [min, max) +[inline] +pub fn (mut rng WyRandRNG) u64_in_range(min, max u64) u64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.u64n(max - min) +} + +// rng.int() returns a pseudorandom 32-bit int (which may be negative) +[inline] +pub fn (mut rng WyRandRNG) int() int { + return int(rng.u32()) +} + +// rng.i64() returns a pseudorandom 64-bit i64 (which may be negative) +[inline] +pub fn (mut rng WyRandRNG) i64() i64 { + return i64(rng.u64()) +} + +// rng.int31() returns a pseudorandom 31-bit int which is non-negative +[inline] +pub fn (mut rng WyRandRNG) int31() int { + return int(rng.u32() & u31_mask) // Set the 32nd bit to 0. +} + +// rng.int63() returns a pseudorandom 63-bit int which is non-negative +[inline] +pub fn (mut rng WyRandRNG) int63() i64 { + return i64(rng.u64() & u63_mask) // Set the 64th bit to 0. +} + +// rng.intn(max) returns a pseudorandom int that lies in [0, max) +[inline] +pub fn (mut rng WyRandRNG) intn(max int) int { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return int(rng.u32n(max)) +} + +// rng.i64n(max) returns a pseudorandom int that lies in [0, max) +[inline] +pub fn (mut rng WyRandRNG) i64n(max i64) i64 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return i64(rng.u64n(max)) +} + +// rng.int_in_range(min, max) returns a pseudorandom int that lies in [min, max) +[inline] +pub fn (mut rng WyRandRNG) int_in_range(min, max int) int { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + // This supports negative ranges like [-10, -5) because the difference is positive + return min + rng.intn(max - min) +} + +// rng.i64_in_range(min, max) returns a pseudorandom i64 that lies in [min, max) +[inline] +pub fn (mut rng WyRandRNG) i64_in_range(min, max i64) i64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.i64n(max - min) +} + +// rng.f32() returns a pseudorandom f32 value between 0.0 (inclusive) and 1.0 (exclusive) i.e [0, 1) +[inline] +pub fn (mut rng WyRandRNG) f32() f32 { + return f32(rng.u32()) / max_u32_as_f32 +} + +// rng.f64() returns a pseudorandom f64 value between 0.0 (inclusive) and 1.0 (exclusive) i.e [0, 1) +[inline] +pub fn (mut rng WyRandRNG) f64() f64 { + return f64(rng.u64()) / max_u64_as_f64 +} + +// rng.f32n() returns a pseudorandom f32 value in [0, max) +[inline] +pub fn (mut rng WyRandRNG) f32n(max f32) f32 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return rng.f32() * max +} + +// rng.f64n() returns a pseudorandom f64 value in [0, max) +[inline] +pub fn (mut rng WyRandRNG) f64n(max f64) f64 { + if max <= 0 { + eprintln('max has to be positive.') + exit(1) + } + return rng.f64() * max +} + +// rng.f32_in_range(min, max) returns a pseudorandom f32 that lies in [min, max) +[inline] +pub fn (mut rng WyRandRNG) f32_in_range(min, max f32) f32 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.f32n(max - min) +} + +// rng.i64_in_range(min, max) returns a pseudorandom i64 that lies in [min, max) +[inline] +pub fn (mut rng WyRandRNG) f64_in_range(min, max f64) f64 { + if max <= min { + eprintln('max must be greater than min') + exit(1) + } + return min + rng.f64n(max - min) +} diff --git a/vlib/rand/wyrand_test.v b/vlib/rand/wyrand_test.v new file mode 100644 index 0000000000..20350023c0 --- /dev/null +++ b/vlib/rand/wyrand_test.v @@ -0,0 +1,330 @@ +import rand +import math + +const ( + range_limit = 40 + value_count = 1000 + seeds = [[u32(42), 0], [u32(256), 0]] +) + +const ( + sample_size = 1000 + stats_epsilon = 0.05 + inv_sqrt_12 = 1.0 / math.sqrt(12) +) + +fn gen_randoms(seed_data []u32, bound int) []u64 { + bound_u64 := u64(bound) + mut randoms := [u64(0)].repeat(20) + mut rnd := rand.WyRandRNG{} + rnd.seed(seed_data) + for i in 0 .. 20 { + randoms[i] = rnd.u64n(bound_u64) + } + return randoms +} + +fn test_wyrand_reproducibility() { + seed_data := rand.time_seed_array(2) + randoms1 := gen_randoms(seed_data, 1000) + randoms2 := gen_randoms(seed_data, 1000) + assert randoms1.len == randoms2.len + len := randoms1.len + for i in 0 .. len { + assert randoms1[i] == randoms2[i] + } +} + +// TODO: use the `in` syntax and remove this function +// after generics has been completely implemented +fn found(value u64, arr []u64) bool { + for item in arr { + if value == item { + return true + } + } + return false +} + +fn test_wyrand_variability() { + // If this test fails and if it is certainly not the implementation + // at fault, try changing the seed values. Repeated values are + // improbable but not impossible. + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + mut values := []u64{cap: value_count} + for i in 0 .. value_count { + value := rng.u64() + assert !found(value, values) + assert values.len == i + values << value + } + } +} + +fn check_uniformity_u64(rng rand.WyRandRNG, range u64) { + range_f64 := f64(range) + expected_mean := range_f64 / 2.0 + mut variance := 0.0 + for _ in 0 .. sample_size { + diff := f64(rng.u64n(range)) - expected_mean + variance += diff * diff + } + variance /= sample_size - 1 + sigma := math.sqrt(variance) + expected_sigma := range_f64 * inv_sqrt_12 + error := (sigma - expected_sigma) / expected_sigma + assert math.abs(error) < stats_epsilon +} + +fn test_wyrand_uniformity_u64() { + ranges := [14019545, 80240, 130] + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + for range in ranges { + check_uniformity_u64(rng, range) + } + } +} + +fn check_uniformity_f64(rng rand.WyRandRNG) { + expected_mean := 0.5 + mut variance := 0.0 + for _ in 0 .. sample_size { + diff := rng.f64() - expected_mean + variance += diff * diff + } + variance /= sample_size - 1 + sigma := math.sqrt(variance) + expected_sigma := inv_sqrt_12 + error := (sigma - expected_sigma) / expected_sigma + assert math.abs(error) < stats_epsilon +} + +fn test_wyrand_uniformity_f64() { + // The f64 version + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + check_uniformity_f64(rng) + } +} + +fn test_wyrand_u32n() { + max := 16384 + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u32n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_wyrand_u64n() { + max := u64(379091181005) + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u64n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_wyrand_u32_in_range() { + max := 484468466 + min := 316846 + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u32_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_wyrand_u64_in_range() { + max := u64(216468454685163) + min := u64(6848646868) + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.u64_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_wyrand_int31() { + max_u31 := 0x7FFFFFFF + sign_mask := 0x80000000 + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int31() + assert value >= 0 + assert value <= max_u31 + // This statement ensures that the sign bit is zero + assert (value & sign_mask) == 0 + } + } +} + +fn test_wyrand_int63() { + max_u63 := i64(0x7FFFFFFFFFFFFFFF) + sign_mask := i64(0x8000000000000000) + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int63() + assert value >= 0 + assert value <= max_u63 + assert (value & sign_mask) == 0 + } + } +} + +fn test_wyrand_intn() { + max := 2525642 + for seed in seeds { + rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.intn(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_wyrand_i64n() { + max := i64(3246727724653636) + for seed in seeds { + rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.i64n(max) + assert value >= 0 + assert value < max + } + } +} + +fn test_wyrand_int_in_range() { + min := -4252 + max := 1034 + for seed in seeds { + rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.int_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_wyrand_i64_in_range() { + min := i64(-24095) + max := i64(324058) + for seed in seeds { + rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.i64_in_range(min, max) + assert value >= min + assert value < max + } + } +} + +fn test_wyrand_f32() { + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= 0.0 + assert value < 1.0 + } + } +} + +fn test_wyrand_f64() { + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= 0.0 + assert value < 1.0 + } + } +} + +fn test_wyrand_f32n() { + max := f32(357.0) + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= 0.0 + assert value < max + } + } +} + +fn test_wyrand_f64n() { + max := 1.52e6 + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= 0.0 + assert value < max + } + } +} + +fn test_wyrand_f32_in_range() { + min := f32(-24.0) + max := f32(125.0) + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f32() + assert value >= min + assert value < max + } + } +} + +fn test_wyrand_f64_in_range() { + min := -548.7 + max := 5015.2 + for seed in seeds { + mut rng := rand.WyRandRNG{} + rng.seed(seed) + for _ in 0 .. range_limit { + value := rng.f64() + assert value >= min + assert value < max + } + } +}