zig

fork of https://codeberg.org/ziglang/zig
Log | Files | Refs | README | LICENSE

commit 32f30399e5cd42af2d670321979aa424803b1819 (tree)
parent 54bbc73f8502fe073d385361ddb34a43d12eec39
Author: Andrew Kelley <andrew@ziglang.org>
Date:   Fri,  9 Feb 2024 13:42:04 -0800

Merge pull request #18867 from e4m2/random

std.rand: Move to std.Random
Diffstat:
MCMakeLists.txt | 2+-
Mlib/build_runner.zig | 4++--
Mlib/compiler_rt/paritydi2_test.zig | 2+-
Mlib/compiler_rt/paritysi2_test.zig | 2+-
Mlib/compiler_rt/parityti2_test.zig | 2+-
Mlib/compiler_rt/popcountdi2_test.zig | 2+-
Mlib/compiler_rt/popcountsi2_test.zig | 2+-
Mlib/compiler_rt/popcountti2_test.zig | 2+-
Mlib/compiler_rt/udivmodei4.zig | 2+-
Alib/std/Random.zig | 438+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/std/Random/Ascon.zig | 58++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/std/Random/ChaCha.zig | 97+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/std/Random/Isaac64.zig | 234+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/std/Random/Pcg.zig | 121+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/std/Random/RomuTrio.zig | 131+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/std/Random/Sfc64.zig | 131+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/std/Random/SplitMix64.zig | 21+++++++++++++++++++++
Alib/std/Random/Xoroshiro128.zig | 146+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/std/Random/Xoshiro256.zig | 145+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/std/Random/benchmark.zig | 217+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/std/Random/test.zig | 473+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/std/Random/ziggurat.zig | 170+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Mlib/std/Thread/RwLock.zig | 2+-
Mlib/std/crypto/benchmark.zig | 2+-
Mlib/std/crypto/kyber_d00.zig | 2+-
Mlib/std/crypto/tlcsprng.zig | 4++--
Mlib/std/hash/benchmark.zig | 2+-
Mlib/std/hash_map.zig | 4++--
Mlib/std/heap/arena_allocator.zig | 2+-
Mlib/std/io/test.zig | 2+-
Mlib/std/math/big/rational.zig | 2+-
Mlib/std/mem.zig | 2+-
Mlib/std/priority_dequeue.zig | 14+++++++-------
Dlib/std/rand.zig | 460-------------------------------------------------------------------------------
Dlib/std/rand/Ascon.zig | 59-----------------------------------------------------------
Dlib/std/rand/ChaCha.zig | 98-------------------------------------------------------------------------------
Dlib/std/rand/Isaac64.zig | 235-------------------------------------------------------------------------------
Dlib/std/rand/Pcg.zig | 122-------------------------------------------------------------------------------
Dlib/std/rand/RomuTrio.zig | 132-------------------------------------------------------------------------------
Dlib/std/rand/Sfc64.zig | 132-------------------------------------------------------------------------------
Dlib/std/rand/Xoroshiro128.zig | 147-------------------------------------------------------------------------------
Dlib/std/rand/Xoshiro256.zig | 146-------------------------------------------------------------------------------
Dlib/std/rand/benchmark.zig | 217-------------------------------------------------------------------------------
Dlib/std/rand/test.zig | 473-------------------------------------------------------------------------------
Dlib/std/rand/ziggurat.zig | 170-------------------------------------------------------------------------------
Mlib/std/sort.zig | 2+-
Mlib/std/std.zig | 5+++--
Mlib/std/treap.zig | 8++++----
Msrc/reduce.zig | 4++--
Msrc/resinator/compile.zig | 2+-
50 files changed, 2422 insertions(+), 2430 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt @@ -299,7 +299,7 @@ set(ZIG_STAGE2_SOURCES "${CMAKE_SOURCE_DIR}/lib/std/Progress.zig" "${CMAKE_SOURCE_DIR}/lib/std/pdb.zig" "${CMAKE_SOURCE_DIR}/lib/std/process.zig" - "${CMAKE_SOURCE_DIR}/lib/std/rand.zig" + "${CMAKE_SOURCE_DIR}/lib/std/Random.zig" "${CMAKE_SOURCE_DIR}/lib/std/sort.zig" "${CMAKE_SOURCE_DIR}/lib/compiler_rt.zig" "${CMAKE_SOURCE_DIR}/lib/compiler_rt/absv.zig" diff --git a/lib/build_runner.zig b/lib/build_runner.zig @@ -420,7 +420,7 @@ fn runStepNames( const starting_steps = try arena.dupe(*Step, step_stack.keys()); - var rng = std.rand.DefaultPrng.init(seed); + var rng = std.Random.DefaultPrng.init(seed); const rand = rng.random(); rand.shuffle(*Step, starting_steps); @@ -836,7 +836,7 @@ fn constructGraphAndCheckForDependencyLoop( b: *std.Build, s: *Step, step_stack: *std.AutoArrayHashMapUnmanaged(*Step, void), - rand: std.rand.Random, + rand: std.Random, ) !void { switch (s.state) { .precheck_started => { diff --git a/lib/compiler_rt/paritydi2_test.zig b/lib/compiler_rt/paritydi2_test.zig @@ -26,7 +26,7 @@ test "paritydi2" { try test__paritydi2(@bitCast(@as(u64, 0xffffffff_fffffffe))); try test__paritydi2(@bitCast(@as(u64, 0xffffffff_ffffffff))); - const RndGen = std.rand.DefaultPrng; + const RndGen = std.Random.DefaultPrng; var rnd = RndGen.init(42); var i: u32 = 0; while (i < 10_000) : (i += 1) { diff --git a/lib/compiler_rt/paritysi2_test.zig b/lib/compiler_rt/paritysi2_test.zig @@ -26,7 +26,7 @@ test "paritysi2" { try test__paritysi2(@bitCast(@as(u32, 0xfffffffe))); try test__paritysi2(@bitCast(@as(u32, 0xffffffff))); - const RndGen = std.rand.DefaultPrng; + const RndGen = std.Random.DefaultPrng; var rnd = RndGen.init(42); var i: u32 = 0; while (i < 10_000) : (i += 1) { diff --git a/lib/compiler_rt/parityti2_test.zig b/lib/compiler_rt/parityti2_test.zig @@ -26,7 +26,7 @@ test "parityti2" { try test__parityti2(@bitCast(@as(u128, 0xffffffff_ffffffff_ffffffff_fffffffe))); try test__parityti2(@bitCast(@as(u128, 0xffffffff_ffffffff_ffffffff_ffffffff))); - const RndGen = std.rand.DefaultPrng; + const RndGen = std.Random.DefaultPrng; var rnd = RndGen.init(42); var i: u32 = 0; while (i < 10_000) : (i += 1) { diff --git a/lib/compiler_rt/popcountdi2_test.zig b/lib/compiler_rt/popcountdi2_test.zig @@ -25,7 +25,7 @@ test "popcountdi2" { try test__popcountdi2(@as(i64, @bitCast(@as(u64, 0xffffffff_fffffffe)))); try test__popcountdi2(@as(i64, @bitCast(@as(u64, 0xffffffff_ffffffff)))); - const RndGen = std.rand.DefaultPrng; + const RndGen = std.Random.DefaultPrng; var rnd = RndGen.init(42); var i: u32 = 0; while (i < 10_000) : (i += 1) { diff --git a/lib/compiler_rt/popcountsi2_test.zig b/lib/compiler_rt/popcountsi2_test.zig @@ -25,7 +25,7 @@ test "popcountsi2" { try test__popcountsi2(@as(i32, @bitCast(@as(u32, 0xfffffffe)))); try test__popcountsi2(@as(i32, @bitCast(@as(u32, 0xffffffff)))); - const RndGen = std.rand.DefaultPrng; + const RndGen = std.Random.DefaultPrng; var rnd = RndGen.init(42); var i: u32 = 0; while (i < 10_000) : (i += 1) { diff --git a/lib/compiler_rt/popcountti2_test.zig b/lib/compiler_rt/popcountti2_test.zig @@ -25,7 +25,7 @@ test "popcountti2" { try test__popcountti2(@as(i128, @bitCast(@as(u128, 0xffffffff_ffffffff_ffffffff_fffffffe)))); try test__popcountti2(@as(i128, @bitCast(@as(u128, 0xffffffff_ffffffff_ffffffff_ffffffff)))); - const RndGen = std.rand.DefaultPrng; + const RndGen = std.Random.DefaultPrng; var rnd = RndGen.init(42); var i: u32 = 0; while (i < 10_000) : (i += 1) { diff --git a/lib/compiler_rt/udivmodei4.zig b/lib/compiler_rt/udivmodei4.zig @@ -132,7 +132,7 @@ test "__udivei4/__umodei4" { if (builtin.zig_backend == .stage2_c) return error.SkipZigTest; if (builtin.zig_backend == .stage2_x86_64) return error.SkipZigTest; - const RndGen = std.rand.DefaultPrng; + const RndGen = std.Random.DefaultPrng; var rnd = RndGen.init(42); var i: usize = 10000; while (i > 0) : (i -= 1) { diff --git a/lib/std/Random.zig b/lib/std/Random.zig @@ -0,0 +1,438 @@ +//! The engines provided here should be initialized from an external source. +//! For a thread-local cryptographically secure pseudo random number generator, +//! use `std.crypto.random`. +//! Be sure to use a CSPRNG when required, otherwise using a normal PRNG will +//! be faster and use substantially less stack space. + +const std = @import("std.zig"); +const math = std.math; +const mem = std.mem; +const assert = std.debug.assert; +const maxInt = std.math.maxInt; +pub const Random = @This(); // Remove pub when `std.rand` namespace is removed. + +/// Fast unbiased random numbers. +pub const DefaultPrng = Xoshiro256; + +/// Cryptographically secure random numbers. +pub const DefaultCsprng = ChaCha; + +pub const Ascon = @import("Random/Ascon.zig"); +pub const ChaCha = @import("Random/ChaCha.zig"); + +pub const Isaac64 = @import("Random/Isaac64.zig"); +pub const Pcg = @import("Random/Pcg.zig"); +pub const Xoroshiro128 = @import("Random/Xoroshiro128.zig"); +pub const Xoshiro256 = @import("Random/Xoshiro256.zig"); +pub const Sfc64 = @import("Random/Sfc64.zig"); +pub const RomuTrio = @import("Random/RomuTrio.zig"); +pub const SplitMix64 = @import("Random/SplitMix64.zig"); +pub const ziggurat = @import("Random/ziggurat.zig"); + +ptr: *anyopaque, +fillFn: *const fn (ptr: *anyopaque, buf: []u8) void, + +pub fn init(pointer: anytype, comptime fillFn: fn (ptr: @TypeOf(pointer), buf: []u8) void) Random { + const Ptr = @TypeOf(pointer); + assert(@typeInfo(Ptr) == .Pointer); // Must be a pointer + assert(@typeInfo(Ptr).Pointer.size == .One); // Must be a single-item pointer + assert(@typeInfo(@typeInfo(Ptr).Pointer.child) == .Struct); // Must point to a struct + const gen = struct { + fn fill(ptr: *anyopaque, buf: []u8) void { + const self: Ptr = @ptrCast(@alignCast(ptr)); + fillFn(self, buf); + } + }; + + return .{ + .ptr = pointer, + .fillFn = gen.fill, + }; +} + +/// Read random bytes into the specified buffer until full. +pub fn bytes(r: Random, buf: []u8) void { + r.fillFn(r.ptr, buf); +} + +pub fn boolean(r: Random) bool { + return r.int(u1) != 0; +} + +/// Returns a random value from an enum, evenly distributed. +/// +/// Note that this will not yield consistent results across all targets +/// due to dependence on the representation of `usize` as an index. +/// See `enumValueWithIndex` for further commentary. +pub inline fn enumValue(r: Random, comptime EnumType: type) EnumType { + return r.enumValueWithIndex(EnumType, usize); +} + +/// Returns a random value from an enum, evenly distributed. +/// +/// An index into an array of all named values is generated using the +/// specified `Index` type to determine the return value. +/// This allows for results to be independent of `usize` representation. +/// +/// Prefer `enumValue` if this isn't important. +/// +/// See `uintLessThan`, which this function uses in most cases, +/// for commentary on the runtime of this function. +pub fn enumValueWithIndex(r: Random, comptime EnumType: type, comptime Index: type) EnumType { + comptime assert(@typeInfo(EnumType) == .Enum); + + // We won't use int -> enum casting because enum elements can have + // arbitrary values. Instead we'll randomly pick one of the type's values. + const values = comptime std.enums.values(EnumType); + comptime assert(values.len > 0); // can't return anything + comptime assert(maxInt(Index) >= values.len - 1); // can't access all values + comptime if (values.len == 1) return values[0]; + + const index = if (comptime values.len - 1 == maxInt(Index)) + r.int(Index) + else + r.uintLessThan(Index, values.len); + + const MinInt = MinArrayIndex(Index); + return values[@as(MinInt, @intCast(index))]; +} + +/// Returns a random int `i` such that `minInt(T) <= i <= maxInt(T)`. +/// `i` is evenly distributed. +pub fn int(r: Random, comptime T: type) T { + const bits = @typeInfo(T).Int.bits; + const UnsignedT = std.meta.Int(.unsigned, bits); + const ceil_bytes = comptime std.math.divCeil(u16, bits, 8) catch unreachable; + const ByteAlignedT = std.meta.Int(.unsigned, ceil_bytes * 8); + + var rand_bytes: [ceil_bytes]u8 = undefined; + r.bytes(&rand_bytes); + + // use LE instead of native endian for better portability maybe? + // TODO: endian portability is pointless if the underlying prng isn't endian portable. + // TODO: document the endian portability of this library. + const byte_aligned_result = mem.readInt(ByteAlignedT, &rand_bytes, .little); + const unsigned_result: UnsignedT = @truncate(byte_aligned_result); + return @bitCast(unsigned_result); +} + +/// Constant-time implementation off `uintLessThan`. +/// The results of this function may be biased. +pub fn uintLessThanBiased(r: Random, comptime T: type, less_than: T) T { + comptime assert(@typeInfo(T).Int.signedness == .unsigned); + assert(0 < less_than); + return limitRangeBiased(T, r.int(T), less_than); +} + +/// Returns an evenly distributed random unsigned integer `0 <= i < less_than`. +/// This function assumes that the underlying `fillFn` produces evenly distributed values. +/// Within this assumption, the runtime of this function is exponentially distributed. +/// If `fillFn` were backed by a true random generator, +/// the runtime of this function would technically be unbounded. +/// However, if `fillFn` is backed by any evenly distributed pseudo random number generator, +/// this function is guaranteed to return. +/// If you need deterministic runtime bounds, use `uintLessThanBiased`. +pub fn uintLessThan(r: Random, comptime T: type, less_than: T) T { + comptime assert(@typeInfo(T).Int.signedness == .unsigned); + const bits = @typeInfo(T).Int.bits; + assert(0 < less_than); + + // adapted from: + // http://www.pcg-random.org/posts/bounded-rands.html + // "Lemire's (with an extra tweak from me)" + var x = r.int(T); + var m = math.mulWide(T, x, less_than); + var l: T = @truncate(m); + if (l < less_than) { + var t = -%less_than; + + if (t >= less_than) { + t -= less_than; + if (t >= less_than) { + t %= less_than; + } + } + while (l < t) { + x = r.int(T); + m = math.mulWide(T, x, less_than); + l = @truncate(m); + } + } + return @intCast(m >> bits); +} + +/// Constant-time implementation off `uintAtMost`. +/// The results of this function may be biased. +pub fn uintAtMostBiased(r: Random, comptime T: type, at_most: T) T { + assert(@typeInfo(T).Int.signedness == .unsigned); + if (at_most == maxInt(T)) { + // have the full range + return r.int(T); + } + return r.uintLessThanBiased(T, at_most + 1); +} + +/// Returns an evenly distributed random unsigned integer `0 <= i <= at_most`. +/// See `uintLessThan`, which this function uses in most cases, +/// for commentary on the runtime of this function. +pub fn uintAtMost(r: Random, comptime T: type, at_most: T) T { + assert(@typeInfo(T).Int.signedness == .unsigned); + if (at_most == maxInt(T)) { + // have the full range + return r.int(T); + } + return r.uintLessThan(T, at_most + 1); +} + +/// Constant-time implementation off `intRangeLessThan`. +/// The results of this function may be biased. +pub fn intRangeLessThanBiased(r: Random, comptime T: type, at_least: T, less_than: T) T { + assert(at_least < less_than); + const info = @typeInfo(T).Int; + if (info.signedness == .signed) { + // Two's complement makes this math pretty easy. + const UnsignedT = std.meta.Int(.unsigned, info.bits); + const lo: UnsignedT = @bitCast(at_least); + const hi: UnsignedT = @bitCast(less_than); + const result = lo +% r.uintLessThanBiased(UnsignedT, hi -% lo); + return @bitCast(result); + } else { + // The signed implementation would work fine, but we can use stricter arithmetic operators here. + return at_least + r.uintLessThanBiased(T, less_than - at_least); + } +} + +/// Returns an evenly distributed random integer `at_least <= i < less_than`. +/// See `uintLessThan`, which this function uses in most cases, +/// for commentary on the runtime of this function. +pub fn intRangeLessThan(r: Random, comptime T: type, at_least: T, less_than: T) T { + assert(at_least < less_than); + const info = @typeInfo(T).Int; + if (info.signedness == .signed) { + // Two's complement makes this math pretty easy. + const UnsignedT = std.meta.Int(.unsigned, info.bits); + const lo: UnsignedT = @bitCast(at_least); + const hi: UnsignedT = @bitCast(less_than); + const result = lo +% r.uintLessThan(UnsignedT, hi -% lo); + return @bitCast(result); + } else { + // The signed implementation would work fine, but we can use stricter arithmetic operators here. + return at_least + r.uintLessThan(T, less_than - at_least); + } +} + +/// Constant-time implementation off `intRangeAtMostBiased`. +/// The results of this function may be biased. +pub fn intRangeAtMostBiased(r: Random, comptime T: type, at_least: T, at_most: T) T { + assert(at_least <= at_most); + const info = @typeInfo(T).Int; + if (info.signedness == .signed) { + // Two's complement makes this math pretty easy. + const UnsignedT = std.meta.Int(.unsigned, info.bits); + const lo: UnsignedT = @bitCast(at_least); + const hi: UnsignedT = @bitCast(at_most); + const result = lo +% r.uintAtMostBiased(UnsignedT, hi -% lo); + return @bitCast(result); + } else { + // The signed implementation would work fine, but we can use stricter arithmetic operators here. + return at_least + r.uintAtMostBiased(T, at_most - at_least); + } +} + +/// Returns an evenly distributed random integer `at_least <= i <= at_most`. +/// See `uintLessThan`, which this function uses in most cases, +/// for commentary on the runtime of this function. +pub fn intRangeAtMost(r: Random, comptime T: type, at_least: T, at_most: T) T { + assert(at_least <= at_most); + const info = @typeInfo(T).Int; + if (info.signedness == .signed) { + // Two's complement makes this math pretty easy. + const UnsignedT = std.meta.Int(.unsigned, info.bits); + const lo: UnsignedT = @bitCast(at_least); + const hi: UnsignedT = @bitCast(at_most); + const result = lo +% r.uintAtMost(UnsignedT, hi -% lo); + return @bitCast(result); + } else { + // The signed implementation would work fine, but we can use stricter arithmetic operators here. + return at_least + r.uintAtMost(T, at_most - at_least); + } +} + +/// Return a floating point value evenly distributed in the range [0, 1). +pub fn float(r: Random, comptime T: type) T { + // Generate a uniformly random value for the mantissa. + // Then generate an exponentially biased random value for the exponent. + // This covers every possible value in the range. + switch (T) { + f32 => { + // Use 23 random bits for the mantissa, and the rest for the exponent. + // If all 41 bits are zero, generate additional random bits, until a + // set bit is found, or 126 bits have been generated. + const rand = r.int(u64); + var rand_lz = @clz(rand); + if (rand_lz >= 41) { + // TODO: when #5177 or #489 is implemented, + // tell the compiler it is unlikely (1/2^41) to reach this point. + // (Same for the if branch and the f64 calculations below.) + rand_lz = 41 + @clz(r.int(u64)); + if (rand_lz == 41 + 64) { + // It is astronomically unlikely to reach this point. + rand_lz += @clz(r.int(u32) | 0x7FF); + } + } + const mantissa: u23 = @truncate(rand); + const exponent = @as(u32, 126 - rand_lz) << 23; + return @bitCast(exponent | mantissa); + }, + f64 => { + // Use 52 random bits for the mantissa, and the rest for the exponent. + // If all 12 bits are zero, generate additional random bits, until a + // set bit is found, or 1022 bits have been generated. + const rand = r.int(u64); + var rand_lz: u64 = @clz(rand); + if (rand_lz >= 12) { + rand_lz = 12; + while (true) { + // It is astronomically unlikely for this loop to execute more than once. + const addl_rand_lz = @clz(r.int(u64)); + rand_lz += addl_rand_lz; + if (addl_rand_lz != 64) { + break; + } + if (rand_lz >= 1022) { + rand_lz = 1022; + break; + } + } + } + const mantissa = rand & 0xFFFFFFFFFFFFF; + const exponent = (1022 - rand_lz) << 52; + return @bitCast(exponent | mantissa); + }, + else => @compileError("unknown floating point type"), + } +} + +/// Return a floating point value normally distributed with mean = 0, stddev = 1. +/// +/// To use different parameters, use: floatNorm(...) * desiredStddev + desiredMean. +pub fn floatNorm(r: Random, comptime T: type) T { + const value = ziggurat.next_f64(r, ziggurat.NormDist); + switch (T) { + f32 => return @floatCast(value), + f64 => return value, + else => @compileError("unknown floating point type"), + } +} + +/// Return an exponentially distributed float with a rate parameter of 1. +/// +/// To use a different rate parameter, use: floatExp(...) / desiredRate. +pub fn floatExp(r: Random, comptime T: type) T { + const value = ziggurat.next_f64(r, ziggurat.ExpDist); + switch (T) { + f32 => return @floatCast(value), + f64 => return value, + else => @compileError("unknown floating point type"), + } +} + +/// Shuffle a slice into a random order. +/// +/// Note that this will not yield consistent results across all targets +/// due to dependence on the representation of `usize` as an index. +/// See `shuffleWithIndex` for further commentary. +pub inline fn shuffle(r: Random, comptime T: type, buf: []T) void { + r.shuffleWithIndex(T, buf, usize); +} + +/// Shuffle a slice into a random order, using an index of a +/// specified type to maintain distribution across targets. +/// Asserts the index type can represent `buf.len`. +/// +/// Indexes into the slice are generated using the specified `Index` +/// type, which determines distribution properties. This allows for +/// results to be independent of `usize` representation. +/// +/// Prefer `shuffle` if this isn't important. +/// +/// See `intRangeLessThan`, which this function uses, +/// for commentary on the runtime of this function. +pub fn shuffleWithIndex(r: Random, comptime T: type, buf: []T, comptime Index: type) void { + const MinInt = MinArrayIndex(Index); + if (buf.len < 2) { + return; + } + + // `i <= j < max <= maxInt(MinInt)` + const max: MinInt = @intCast(buf.len); + var i: MinInt = 0; + while (i < max - 1) : (i += 1) { + const j: MinInt = @intCast(r.intRangeLessThan(Index, i, max)); + mem.swap(T, &buf[i], &buf[j]); + } +} + +/// Randomly selects an index into `proportions`, where the likelihood of each +/// index is weighted by that proportion. +/// It is more likely for the index of the last proportion to be returned +/// than the index of the first proportion in the slice, and vice versa. +/// +/// This is useful for selecting an item from a slice where weights are not equal. +/// `T` must be a numeric type capable of holding the sum of `proportions`. +pub fn weightedIndex(r: Random, comptime T: type, proportions: []const T) usize { + // This implementation works by summing the proportions and picking a + // random point in [0, sum). We then loop over the proportions, + // accumulating until our accumulator is greater than the random point. + + const sum = s: { + var sum: T = 0; + for (proportions) |v| sum += v; + break :s sum; + }; + + const point = switch (@typeInfo(T)) { + .Int => |int_info| switch (int_info.signedness) { + .signed => r.intRangeLessThan(T, 0, sum), + .unsigned => r.uintLessThan(T, sum), + }, + // take care that imprecision doesn't lead to a value slightly greater than sum + .Float => @min(r.float(T) * sum, sum - std.math.floatEps(T)), + else => @compileError("weightedIndex does not support proportions of type " ++ + @typeName(T)), + }; + + assert(point < sum); + + var accumulator: T = 0; + for (proportions, 0..) |p, index| { + accumulator += p; + if (point < accumulator) return index; + } else unreachable; +} + +/// Convert a random integer 0 <= random_int <= maxValue(T), +/// into an integer 0 <= result < less_than. +/// This function introduces a minor bias. +pub fn limitRangeBiased(comptime T: type, random_int: T, less_than: T) T { + comptime assert(@typeInfo(T).Int.signedness == .unsigned); + const bits = @typeInfo(T).Int.bits; + + // adapted from: + // http://www.pcg-random.org/posts/bounded-rands.html + // "Integer Multiplication (Biased)" + const m = math.mulWide(T, random_int, less_than); + return @intCast(m >> bits); +} + +/// Returns the smallest of `Index` and `usize`. +fn MinArrayIndex(comptime Index: type) type { + const index_info = @typeInfo(Index).Int; + assert(index_info.signedness == .unsigned); + return if (index_info.bits >= @typeInfo(usize).Int.bits) usize else Index; +} + +test { + std.testing.refAllDecls(@This()); + _ = @import("Random/test.zig"); +} diff --git a/lib/std/Random/Ascon.zig b/lib/std/Random/Ascon.zig @@ -0,0 +1,58 @@ +//! CSPRNG based on the Reverie construction, a permutation-based PRNG +//! with forward security, instantiated with the Ascon(128,12,8) permutation. +//! +//! Compared to ChaCha, this PRNG has a much smaller state, and can be +//! a better choice for constrained environments. +//! +//! References: +//! - A Robust and Sponge-Like PRNG with Improved Efficiency https://eprint.iacr.org/2016/886.pdf +//! - Ascon https://ascon.iaik.tugraz.at/files/asconv12-nist.pdf + +const std = @import("std"); +const mem = std.mem; +const Self = @This(); + +const Ascon = std.crypto.core.Ascon(.little); + +state: Ascon, + +const rate = 16; +pub const secret_seed_length = 32; + +/// The seed must be uniform, secret and `secret_seed_length` bytes long. +pub fn init(secret_seed: [secret_seed_length]u8) Self { + var self = Self{ .state = Ascon.initXof() }; + self.addEntropy(&secret_seed); + return self; +} + +/// Inserts entropy to refresh the internal state. +pub fn addEntropy(self: *Self, bytes: []const u8) void { + comptime std.debug.assert(secret_seed_length % rate == 0); + var i: usize = 0; + while (i + rate < bytes.len) : (i += rate) { + self.state.addBytes(bytes[i..][0..rate]); + self.state.permuteR(8); + } + if (i != bytes.len) self.state.addBytes(bytes[i..]); + self.state.permute(); +} + +/// Returns a `std.Random` structure backed by the current RNG. +pub fn random(self: *Self) std.Random { + return std.Random.init(self, fill); +} + +/// Fills the buffer with random bytes. +pub fn fill(self: *Self, buf: []u8) void { + var i: usize = 0; + while (true) { + const left = buf.len - i; + const n = @min(left, rate); + self.state.extractBytes(buf[i..][0..n]); + if (left == 0) break; + self.state.permuteR(8); + i += n; + } + self.state.permuteRatchet(6, rate); +} diff --git a/lib/std/Random/ChaCha.zig b/lib/std/Random/ChaCha.zig @@ -0,0 +1,97 @@ +//! CSPRNG based on the ChaCha8 stream cipher, with forward security. +//! +//! References: +//! - Fast-key-erasure random-number generators https://blog.cr.yp.to/20170723-random.html + +const std = @import("std"); +const mem = std.mem; +const Self = @This(); + +const Cipher = std.crypto.stream.chacha.ChaCha8IETF; + +const State = [8 * Cipher.block_length]u8; + +state: State, +offset: usize, + +const nonce = [_]u8{0} ** Cipher.nonce_length; + +pub const secret_seed_length = Cipher.key_length; + +/// The seed must be uniform, secret and `secret_seed_length` bytes long. +pub fn init(secret_seed: [secret_seed_length]u8) Self { + var self = Self{ .state = undefined, .offset = 0 }; + Cipher.stream(&self.state, 0, secret_seed, nonce); + return self; +} + +/// Inserts entropy to refresh the internal state. +pub fn addEntropy(self: *Self, bytes: []const u8) void { + var i: usize = 0; + while (i + Cipher.key_length <= bytes.len) : (i += Cipher.key_length) { + Cipher.xor( + self.state[0..Cipher.key_length], + self.state[0..Cipher.key_length], + 0, + bytes[i..][0..Cipher.key_length].*, + nonce, + ); + } + if (i < bytes.len) { + var k = [_]u8{0} ** Cipher.key_length; + const src = bytes[i..]; + @memcpy(k[0..src.len], src); + Cipher.xor( + self.state[0..Cipher.key_length], + self.state[0..Cipher.key_length], + 0, + k, + nonce, + ); + } + self.refill(); +} + +/// Returns a `std.Random` structure backed by the current RNG. +pub fn random(self: *Self) std.Random { + return std.Random.init(self, fill); +} + +// Refills the buffer with random bytes, overwriting the previous key. +fn refill(self: *Self) void { + Cipher.stream(&self.state, 0, self.state[0..Cipher.key_length].*, nonce); + self.offset = 0; +} + +/// Fills the buffer with random bytes. +pub fn fill(self: *Self, buf_: []u8) void { + const bytes = self.state[Cipher.key_length..]; + var buf = buf_; + + const avail = bytes.len - self.offset; + if (avail > 0) { + // Bytes from the current block + const n = @min(avail, buf.len); + @memcpy(buf[0..n], bytes[self.offset..][0..n]); + @memset(bytes[self.offset..][0..n], 0); + buf = buf[n..]; + self.offset += n; + } + if (buf.len == 0) return; + + self.refill(); + + // Full blocks + while (buf.len >= bytes.len) { + @memcpy(buf[0..bytes.len], bytes); + buf = buf[bytes.len..]; + self.refill(); + } + + // Remaining bytes + if (buf.len > 0) { + @memcpy(buf, bytes[0..buf.len]); + @memset(bytes[0..buf.len], 0); + self.offset = buf.len; + } +} diff --git a/lib/std/Random/Isaac64.zig b/lib/std/Random/Isaac64.zig @@ -0,0 +1,234 @@ +//! ISAAC64 - http://www.burtleburtle.net/bob/rand/isaacafa.html +//! +//! Follows the general idea of the implementation from here with a few shortcuts. +//! https://doc.rust-lang.org/rand/src/rand/prng/isaac64.rs.html + +const std = @import("std"); +const mem = std.mem; +const Isaac64 = @This(); + +r: [256]u64, +m: [256]u64, +a: u64, +b: u64, +c: u64, +i: usize, + +pub fn init(init_s: u64) Isaac64 { + var isaac = Isaac64{ + .r = undefined, + .m = undefined, + .a = undefined, + .b = undefined, + .c = undefined, + .i = undefined, + }; + + // seed == 0 => same result as the unseeded reference implementation + isaac.seed(init_s, 1); + return isaac; +} + +pub fn random(self: *Isaac64) std.Random { + return std.Random.init(self, fill); +} + +fn step(self: *Isaac64, mix: u64, base: usize, comptime m1: usize, comptime m2: usize) void { + const x = self.m[base + m1]; + self.a = mix +% self.m[base + m2]; + + const y = self.a +% self.b +% self.m[@as(usize, @intCast((x >> 3) % self.m.len))]; + self.m[base + m1] = y; + + self.b = x +% self.m[@as(usize, @intCast((y >> 11) % self.m.len))]; + self.r[self.r.len - 1 - base - m1] = self.b; +} + +fn refill(self: *Isaac64) void { + const midpoint = self.r.len / 2; + + self.c +%= 1; + self.b +%= self.c; + + { + var i: usize = 0; + while (i < midpoint) : (i += 4) { + self.step(~(self.a ^ (self.a << 21)), i + 0, 0, midpoint); + self.step(self.a ^ (self.a >> 5), i + 1, 0, midpoint); + self.step(self.a ^ (self.a << 12), i + 2, 0, midpoint); + self.step(self.a ^ (self.a >> 33), i + 3, 0, midpoint); + } + } + + { + var i: usize = 0; + while (i < midpoint) : (i += 4) { + self.step(~(self.a ^ (self.a << 21)), i + 0, midpoint, 0); + self.step(self.a ^ (self.a >> 5), i + 1, midpoint, 0); + self.step(self.a ^ (self.a << 12), i + 2, midpoint, 0); + self.step(self.a ^ (self.a >> 33), i + 3, midpoint, 0); + } + } + + self.i = 0; +} + +fn next(self: *Isaac64) u64 { + if (self.i >= self.r.len) { + self.refill(); + } + + const value = self.r[self.i]; + self.i += 1; + return value; +} + +fn seed(self: *Isaac64, init_s: u64, comptime rounds: usize) void { + // We ignore the multi-pass requirement since we don't currently expose full access to + // seeding the self.m array completely. + @memset(self.m[0..], 0); + self.m[0] = init_s; + + // prescrambled golden ratio constants + var a = [_]u64{ + 0x647c4677a2884b7c, + 0xb9f8b322c73ac862, + 0x8c0ea5053d4712a0, + 0xb29b2e824a595524, + 0x82f053db8355e0ce, + 0x48fe4a0fa5a09315, + 0xae985bf2cbfc89ed, + 0x98f5704f6c44c0ab, + }; + + comptime var i: usize = 0; + inline while (i < rounds) : (i += 1) { + var j: usize = 0; + while (j < self.m.len) : (j += 8) { + comptime var x1: usize = 0; + inline while (x1 < 8) : (x1 += 1) { + a[x1] +%= self.m[j + x1]; + } + + a[0] -%= a[4]; + a[5] ^= a[7] >> 9; + a[7] +%= a[0]; + a[1] -%= a[5]; + a[6] ^= a[0] << 9; + a[0] +%= a[1]; + a[2] -%= a[6]; + a[7] ^= a[1] >> 23; + a[1] +%= a[2]; + a[3] -%= a[7]; + a[0] ^= a[2] << 15; + a[2] +%= a[3]; + a[4] -%= a[0]; + a[1] ^= a[3] >> 14; + a[3] +%= a[4]; + a[5] -%= a[1]; + a[2] ^= a[4] << 20; + a[4] +%= a[5]; + a[6] -%= a[2]; + a[3] ^= a[5] >> 17; + a[5] +%= a[6]; + a[7] -%= a[3]; + a[4] ^= a[6] << 14; + a[6] +%= a[7]; + + comptime var x2: usize = 0; + inline while (x2 < 8) : (x2 += 1) { + self.m[j + x2] = a[x2]; + } + } + } + + @memset(self.r[0..], 0); + self.a = 0; + self.b = 0; + self.c = 0; + self.i = self.r.len; // trigger refill on first value +} + +pub fn fill(self: *Isaac64, buf: []u8) void { + var i: usize = 0; + const aligned_len = buf.len - (buf.len & 7); + + // Fill complete 64-byte segments + while (i < aligned_len) : (i += 8) { + var n = self.next(); + comptime var j: usize = 0; + inline while (j < 8) : (j += 1) { + buf[i + j] = @as(u8, @truncate(n)); + n >>= 8; + } + } + + // Fill trailing, ignoring excess (cut the stream). + if (i != buf.len) { + var n = self.next(); + while (i < buf.len) : (i += 1) { + buf[i] = @as(u8, @truncate(n)); + n >>= 8; + } + } +} + +test "isaac64 sequence" { + var r = Isaac64.init(0); + + // from reference implementation + const seq = [_]u64{ + 0xf67dfba498e4937c, + 0x84a5066a9204f380, + 0xfee34bd5f5514dbb, + 0x4d1664739b8f80d6, + 0x8607459ab52a14aa, + 0x0e78bc5a98529e49, + 0xfe5332822ad13777, + 0x556c27525e33d01a, + 0x08643ca615f3149f, + 0xd0771faf3cb04714, + 0x30e86f68a37b008d, + 0x3074ebc0488a3adf, + 0x270645ea7a2790bc, + 0x5601a0a8d3763c6a, + 0x2f83071f53f325dd, + 0xb9090f3d42d2d2ea, + }; + + for (seq) |s| { + try std.testing.expect(s == r.next()); + } +} + +test "isaac64 fill" { + var r = Isaac64.init(0); + + // from reference implementation + const seq = [_]u64{ + 0xf67dfba498e4937c, + 0x84a5066a9204f380, + 0xfee34bd5f5514dbb, + 0x4d1664739b8f80d6, + 0x8607459ab52a14aa, + 0x0e78bc5a98529e49, + 0xfe5332822ad13777, + 0x556c27525e33d01a, + 0x08643ca615f3149f, + 0xd0771faf3cb04714, + 0x30e86f68a37b008d, + 0x3074ebc0488a3adf, + 0x270645ea7a2790bc, + 0x5601a0a8d3763c6a, + 0x2f83071f53f325dd, + 0xb9090f3d42d2d2ea, + }; + + for (seq) |s| { + var buf0: [8]u8 = undefined; + var buf1: [7]u8 = undefined; + std.mem.writeInt(u64, &buf0, s, .little); + r.fill(&buf1); + try std.testing.expect(std.mem.eql(u8, buf0[0..7], buf1[0..])); + } +} diff --git a/lib/std/Random/Pcg.zig b/lib/std/Random/Pcg.zig @@ -0,0 +1,121 @@ +//! PCG32 - http://www.pcg-random.org/ +//! +//! PRNG + +const std = @import("std"); +const Pcg = @This(); + +const default_multiplier = 6364136223846793005; + +s: u64, +i: u64, + +pub fn init(init_s: u64) Pcg { + var pcg = Pcg{ + .s = undefined, + .i = undefined, + }; + + pcg.seed(init_s); + return pcg; +} + +pub fn random(self: *Pcg) std.Random { + return std.Random.init(self, fill); +} + +fn next(self: *Pcg) u32 { + const l = self.s; + self.s = l *% default_multiplier +% (self.i | 1); + + const xor_s: u32 = @truncate(((l >> 18) ^ l) >> 27); + const rot: u32 = @intCast(l >> 59); + + return (xor_s >> @as(u5, @intCast(rot))) | (xor_s << @as(u5, @intCast((0 -% rot) & 31))); +} + +fn seed(self: *Pcg, init_s: u64) void { + // Pcg requires 128-bits of seed. + var gen = std.Random.SplitMix64.init(init_s); + self.seedTwo(gen.next(), gen.next()); +} + +fn seedTwo(self: *Pcg, init_s: u64, init_i: u64) void { + self.s = 0; + self.i = (init_s << 1) | 1; + self.s = self.s *% default_multiplier +% self.i; + self.s +%= init_i; + self.s = self.s *% default_multiplier +% self.i; +} + +pub fn fill(self: *Pcg, buf: []u8) void { + var i: usize = 0; + const aligned_len = buf.len - (buf.len & 3); + + // Complete 4 byte segments. + while (i < aligned_len) : (i += 4) { + var n = self.next(); + comptime var j: usize = 0; + inline while (j < 4) : (j += 1) { + buf[i + j] = @as(u8, @truncate(n)); + n >>= 8; + } + } + + // Remaining. (cuts the stream) + if (i != buf.len) { + var n = self.next(); + while (i < buf.len) : (i += 1) { + buf[i] = @as(u8, @truncate(n)); + n >>= 8; + } + } +} + +test "pcg sequence" { + var r = Pcg.init(0); + const s0: u64 = 0x9394bf54ce5d79de; + const s1: u64 = 0x84e9c579ef59bbf7; + r.seedTwo(s0, s1); + + const seq = [_]u32{ + 2881561918, + 3063928540, + 1199791034, + 2487695858, + 1479648952, + 3247963454, + }; + + for (seq) |s| { + try std.testing.expect(s == r.next()); + } +} + +test "pcg fill" { + var r = Pcg.init(0); + const s0: u64 = 0x9394bf54ce5d79de; + const s1: u64 = 0x84e9c579ef59bbf7; + r.seedTwo(s0, s1); + + const seq = [_]u32{ + 2881561918, + 3063928540, + 1199791034, + 2487695858, + 1479648952, + 3247963454, + }; + + var i: u32 = 0; + while (i < seq.len) : (i += 2) { + var buf0: [8]u8 = undefined; + std.mem.writeInt(u32, buf0[0..4], seq[i], .little); + std.mem.writeInt(u32, buf0[4..8], seq[i + 1], .little); + + var buf1: [7]u8 = undefined; + r.fill(&buf1); + + try std.testing.expect(std.mem.eql(u8, buf0[0..7], buf1[0..])); + } +} diff --git a/lib/std/Random/RomuTrio.zig b/lib/std/Random/RomuTrio.zig @@ -0,0 +1,131 @@ +// Website: romu-random.org +// Reference paper: http://arxiv.org/abs/2002.11331 +// Beware: this PRNG is trivially predictable. While fast, it should *never* be used for cryptographic purposes. + +const std = @import("std"); +const math = std.math; +const RomuTrio = @This(); + +x_state: u64, +y_state: u64, +z_state: u64, // set to nonzero seed + +pub fn init(init_s: u64) RomuTrio { + var x = RomuTrio{ .x_state = undefined, .y_state = undefined, .z_state = undefined }; + x.seed(init_s); + return x; +} + +pub fn random(self: *RomuTrio) std.Random { + return std.Random.init(self, fill); +} + +fn next(self: *RomuTrio) u64 { + const xp = self.x_state; + const yp = self.y_state; + const zp = self.z_state; + self.x_state = 15241094284759029579 *% zp; + self.y_state = yp -% xp; + self.y_state = std.math.rotl(u64, self.y_state, 12); + self.z_state = zp -% yp; + self.z_state = std.math.rotl(u64, self.z_state, 44); + return xp; +} + +pub fn seedWithBuf(self: *RomuTrio, buf: [24]u8) void { + const seed_buf = @as([3]u64, @bitCast(buf)); + self.x_state = seed_buf[0]; + self.y_state = seed_buf[1]; + self.z_state = seed_buf[2]; +} + +pub fn seed(self: *RomuTrio, init_s: u64) void { + // RomuTrio requires 192-bits of seed. + var gen = std.Random.SplitMix64.init(init_s); + + self.x_state = gen.next(); + self.y_state = gen.next(); + self.z_state = gen.next(); +} + +pub fn fill(self: *RomuTrio, buf: []u8) void { + var i: usize = 0; + const aligned_len = buf.len - (buf.len & 7); + + // Complete 8 byte segments. + while (i < aligned_len) : (i += 8) { + var n = self.next(); + comptime var j: usize = 0; + inline while (j < 8) : (j += 1) { + buf[i + j] = @as(u8, @truncate(n)); + n >>= 8; + } + } + + // Remaining. (cuts the stream) + if (i != buf.len) { + var n = self.next(); + while (i < buf.len) : (i += 1) { + buf[i] = @as(u8, @truncate(n)); + n >>= 8; + } + } +} + +test "RomuTrio sequence" { + // Unfortunately there does not seem to be an official test sequence. + var r = RomuTrio.init(0); + + const seq = [_]u64{ + 16294208416658607535, + 13964609475759908645, + 4703697494102998476, + 3425221541186733346, + 2285772463536419399, + 9454187757529463048, + 13695907680080547496, + 8328236714879408626, + 12323357569716880909, + 12375466223337721820, + }; + + for (seq) |s| { + try std.testing.expectEqual(s, r.next()); + } +} + +test "RomuTrio fill" { + // Unfortunately there does not seem to be an official test sequence. + var r = RomuTrio.init(0); + + const seq = [_]u64{ + 16294208416658607535, + 13964609475759908645, + 4703697494102998476, + 3425221541186733346, + 2285772463536419399, + 9454187757529463048, + 13695907680080547496, + 8328236714879408626, + 12323357569716880909, + 12375466223337721820, + }; + + for (seq) |s| { + var buf0: [8]u8 = undefined; + var buf1: [7]u8 = undefined; + std.mem.writeInt(u64, &buf0, s, .little); + r.fill(&buf1); + try std.testing.expect(std.mem.eql(u8, buf0[0..7], buf1[0..])); + } +} + +test "RomuTrio buf seeding test" { + const buf0 = @as([24]u8, @bitCast([3]u64{ 16294208416658607535, 13964609475759908645, 4703697494102998476 })); + const resulting_state = .{ .x = 16294208416658607535, .y = 13964609475759908645, .z = 4703697494102998476 }; + var r = RomuTrio.init(0); + r.seedWithBuf(buf0); + try std.testing.expect(r.x_state == resulting_state.x); + try std.testing.expect(r.y_state == resulting_state.y); + try std.testing.expect(r.z_state == resulting_state.z); +} diff --git a/lib/std/Random/Sfc64.zig b/lib/std/Random/Sfc64.zig @@ -0,0 +1,131 @@ +//! Sfc64 pseudo-random number generator from Practically Random. +//! Fastest engine of pracrand and smallest footprint. +//! See http://pracrand.sourceforge.net/ + +const std = @import("std"); +const math = std.math; +const Sfc64 = @This(); + +a: u64 = undefined, +b: u64 = undefined, +c: u64 = undefined, +counter: u64 = undefined, + +const Rotation = 24; +const RightShift = 11; +const LeftShift = 3; + +pub fn init(init_s: u64) Sfc64 { + var x = Sfc64{}; + + x.seed(init_s); + return x; +} + +pub fn random(self: *Sfc64) std.Random { + return std.Random.init(self, fill); +} + +fn next(self: *Sfc64) u64 { + const tmp = self.a +% self.b +% self.counter; + self.counter += 1; + self.a = self.b ^ (self.b >> RightShift); + self.b = self.c +% (self.c << LeftShift); + self.c = math.rotl(u64, self.c, Rotation) +% tmp; + return tmp; +} + +fn seed(self: *Sfc64, init_s: u64) void { + self.a = init_s; + self.b = init_s; + self.c = init_s; + self.counter = 1; + var i: u32 = 0; + while (i < 12) : (i += 1) { + _ = self.next(); + } +} + +pub fn fill(self: *Sfc64, buf: []u8) void { + var i: usize = 0; + const aligned_len = buf.len - (buf.len & 7); + + // Complete 8 byte segments. + while (i < aligned_len) : (i += 8) { + var n = self.next(); + comptime var j: usize = 0; + inline while (j < 8) : (j += 1) { + buf[i + j] = @as(u8, @truncate(n)); + n >>= 8; + } + } + + // Remaining. (cuts the stream) + if (i != buf.len) { + var n = self.next(); + while (i < buf.len) : (i += 1) { + buf[i] = @as(u8, @truncate(n)); + n >>= 8; + } + } +} + +test "Sfc64 sequence" { + // Unfortunately there does not seem to be an official test sequence. + var r = Sfc64.init(0); + + const seq = [_]u64{ + 0x3acfa029e3cc6041, + 0xf5b6515bf2ee419c, + 0x1259635894a29b61, + 0xb6ae75395f8ebd6, + 0x225622285ce302e2, + 0x520d28611395cb21, + 0xdb909c818901599d, + 0x8ffd195365216f57, + 0xe8c4ad5e258ac04a, + 0x8f8ef2c89fdb63ca, + 0xf9865b01d98d8e2f, + 0x46555871a65d08ba, + 0x66868677c6298fcd, + 0x2ce15a7e6329f57d, + 0xb2f1833ca91ca79, + 0x4b0890ac9bf453ca, + }; + + for (seq) |s| { + try std.testing.expectEqual(s, r.next()); + } +} + +test "Sfc64 fill" { + // Unfortunately there does not seem to be an official test sequence. + var r = Sfc64.init(0); + + const seq = [_]u64{ + 0x3acfa029e3cc6041, + 0xf5b6515bf2ee419c, + 0x1259635894a29b61, + 0xb6ae75395f8ebd6, + 0x225622285ce302e2, + 0x520d28611395cb21, + 0xdb909c818901599d, + 0x8ffd195365216f57, + 0xe8c4ad5e258ac04a, + 0x8f8ef2c89fdb63ca, + 0xf9865b01d98d8e2f, + 0x46555871a65d08ba, + 0x66868677c6298fcd, + 0x2ce15a7e6329f57d, + 0xb2f1833ca91ca79, + 0x4b0890ac9bf453ca, + }; + + for (seq) |s| { + var buf0: [8]u8 = undefined; + var buf1: [7]u8 = undefined; + std.mem.writeInt(u64, &buf0, s, .little); + r.fill(&buf1); + try std.testing.expect(std.mem.eql(u8, buf0[0..7], buf1[0..])); + } +} diff --git a/lib/std/Random/SplitMix64.zig b/lib/std/Random/SplitMix64.zig @@ -0,0 +1,21 @@ +//! Generator to extend 64-bit seed values into longer sequences. +//! +//! The number of cycles is thus limited to 64-bits regardless of the engine, but this +//! is still plenty for practical purposes. + +const SplitMix64 = @This(); + +s: u64, + +pub fn init(seed: u64) SplitMix64 { + return SplitMix64{ .s = seed }; +} + +pub fn next(self: *SplitMix64) u64 { + self.s +%= 0x9e3779b97f4a7c15; + + var z = self.s; + z = (z ^ (z >> 30)) *% 0xbf58476d1ce4e5b9; + z = (z ^ (z >> 27)) *% 0x94d049bb133111eb; + return z ^ (z >> 31); +} diff --git a/lib/std/Random/Xoroshiro128.zig b/lib/std/Random/Xoroshiro128.zig @@ -0,0 +1,146 @@ +//! Xoroshiro128+ - http://xoroshiro.di.unimi.it/ +//! +//! PRNG + +const std = @import("std"); +const math = std.math; +const Xoroshiro128 = @This(); + +s: [2]u64, + +pub fn init(init_s: u64) Xoroshiro128 { + var x = Xoroshiro128{ .s = undefined }; + + x.seed(init_s); + return x; +} + +pub fn random(self: *Xoroshiro128) std.Random { + return std.Random.init(self, fill); +} + +pub fn next(self: *Xoroshiro128) u64 { + const s0 = self.s[0]; + var s1 = self.s[1]; + const r = s0 +% s1; + + s1 ^= s0; + self.s[0] = math.rotl(u64, s0, @as(u8, 55)) ^ s1 ^ (s1 << 14); + self.s[1] = math.rotl(u64, s1, @as(u8, 36)); + + return r; +} + +// Skip 2^64 places ahead in the sequence +pub fn jump(self: *Xoroshiro128) void { + var s0: u64 = 0; + var s1: u64 = 0; + + const table = [_]u64{ + 0xbeac0467eba5facb, + 0xd86b048b86aa9922, + }; + + inline for (table) |entry| { + var b: usize = 0; + while (b < 64) : (b += 1) { + if ((entry & (@as(u64, 1) << @as(u6, @intCast(b)))) != 0) { + s0 ^= self.s[0]; + s1 ^= self.s[1]; + } + _ = self.next(); + } + } + + self.s[0] = s0; + self.s[1] = s1; +} + +pub fn seed(self: *Xoroshiro128, init_s: u64) void { + // Xoroshiro requires 128-bits of seed. + var gen = std.Random.SplitMix64.init(init_s); + + self.s[0] = gen.next(); + self.s[1] = gen.next(); +} + +pub fn fill(self: *Xoroshiro128, buf: []u8) void { + var i: usize = 0; + const aligned_len = buf.len - (buf.len & 7); + + // Complete 8 byte segments. + while (i < aligned_len) : (i += 8) { + var n = self.next(); + comptime var j: usize = 0; + inline while (j < 8) : (j += 1) { + buf[i + j] = @as(u8, @truncate(n)); + n >>= 8; + } + } + + // Remaining. (cuts the stream) + if (i != buf.len) { + var n = self.next(); + while (i < buf.len) : (i += 1) { + buf[i] = @as(u8, @truncate(n)); + n >>= 8; + } + } +} + +test "xoroshiro sequence" { + var r = Xoroshiro128.init(0); + r.s[0] = 0xaeecf86f7878dd75; + r.s[1] = 0x01cd153642e72622; + + const seq1 = [_]u64{ + 0xb0ba0da5bb600397, + 0x18a08afde614dccc, + 0xa2635b956a31b929, + 0xabe633c971efa045, + 0x9ac19f9706ca3cac, + 0xf62b426578c1e3fb, + }; + + for (seq1) |s| { + try std.testing.expect(s == r.next()); + } + + r.jump(); + + const seq2 = [_]u64{ + 0x95344a13556d3e22, + 0xb4fb32dafa4d00df, + 0xb2011d9ccdcfe2dd, + 0x05679a9b2119b908, + 0xa860a1da7c9cd8a0, + 0x658a96efe3f86550, + }; + + for (seq2) |s| { + try std.testing.expect(s == r.next()); + } +} + +test "xoroshiro fill" { + var r = Xoroshiro128.init(0); + r.s[0] = 0xaeecf86f7878dd75; + r.s[1] = 0x01cd153642e72622; + + const seq = [_]u64{ + 0xb0ba0da5bb600397, + 0x18a08afde614dccc, + 0xa2635b956a31b929, + 0xabe633c971efa045, + 0x9ac19f9706ca3cac, + 0xf62b426578c1e3fb, + }; + + for (seq) |s| { + var buf0: [8]u8 = undefined; + var buf1: [7]u8 = undefined; + std.mem.writeInt(u64, &buf0, s, .little); + r.fill(&buf1); + try std.testing.expect(std.mem.eql(u8, buf0[0..7], buf1[0..])); + } +} diff --git a/lib/std/Random/Xoshiro256.zig b/lib/std/Random/Xoshiro256.zig @@ -0,0 +1,145 @@ +//! Xoshiro256++ - http://xoroshiro.di.unimi.it/ +//! +//! PRNG + +const std = @import("std"); +const math = std.math; +const Xoshiro256 = @This(); + +s: [4]u64, + +pub fn init(init_s: u64) Xoshiro256 { + var x = Xoshiro256{ + .s = undefined, + }; + + x.seed(init_s); + return x; +} + +pub fn random(self: *Xoshiro256) std.Random { + return std.Random.init(self, fill); +} + +pub fn next(self: *Xoshiro256) u64 { + const r = math.rotl(u64, self.s[0] +% self.s[3], 23) +% self.s[0]; + + const t = self.s[1] << 17; + + self.s[2] ^= self.s[0]; + self.s[3] ^= self.s[1]; + self.s[1] ^= self.s[2]; + self.s[0] ^= self.s[3]; + + self.s[2] ^= t; + + self.s[3] = math.rotl(u64, self.s[3], 45); + + return r; +} + +// Skip 2^128 places ahead in the sequence +pub fn jump(self: *Xoshiro256) void { + var s: u256 = 0; + + var table: u256 = 0x39abdc4529b1661ca9582618e03fc9aad5a61266f0c9392c180ec6d33cfd0aba; + + while (table != 0) : (table >>= 1) { + if (@as(u1, @truncate(table)) != 0) { + s ^= @as(u256, @bitCast(self.s)); + } + _ = self.next(); + } + + self.s = @as([4]u64, @bitCast(s)); +} + +pub fn seed(self: *Xoshiro256, init_s: u64) void { + // Xoshiro requires 256-bits of seed. + var gen = std.Random.SplitMix64.init(init_s); + + self.s[0] = gen.next(); + self.s[1] = gen.next(); + self.s[2] = gen.next(); + self.s[3] = gen.next(); +} + +pub fn fill(self: *Xoshiro256, buf: []u8) void { + var i: usize = 0; + const aligned_len = buf.len - (buf.len & 7); + + // Complete 8 byte segments. + while (i < aligned_len) : (i += 8) { + var n = self.next(); + comptime var j: usize = 0; + inline while (j < 8) : (j += 1) { + buf[i + j] = @as(u8, @truncate(n)); + n >>= 8; + } + } + + // Remaining. (cuts the stream) + if (i != buf.len) { + var n = self.next(); + while (i < buf.len) : (i += 1) { + buf[i] = @as(u8, @truncate(n)); + n >>= 8; + } + } +} + +test "xoroshiro sequence" { + if (@import("builtin").zig_backend == .stage2_c) return error.SkipZigTest; + if (@import("builtin").zig_backend == .stage2_x86_64) return error.SkipZigTest; + + var r = Xoshiro256.init(0); + + const seq1 = [_]u64{ + 0x53175d61490b23df, + 0x61da6f3dc380d507, + 0x5c0fdf91ec9a7bfc, + 0x02eebf8c3bbe5e1a, + 0x7eca04ebaf4a5eea, + 0x0543c37757f08d9a, + }; + + for (seq1) |s| { + try std.testing.expect(s == r.next()); + } + + r.jump(); + + const seq2 = [_]u64{ + 0xae1db5c5e27807be, + 0xb584c6a7fd8709fe, + 0xc46a0ee9330fb6e, + 0xdc0c9606f49ed76e, + 0x1f5bb6540f6651fb, + 0x72fa2ca734601488, + }; + + for (seq2) |s| { + try std.testing.expect(s == r.next()); + } +} + +test "xoroshiro fill" { + var r = Xoshiro256.init(0); + + const seq = [_]u64{ + 0x53175d61490b23df, + 0x61da6f3dc380d507, + 0x5c0fdf91ec9a7bfc, + 0x02eebf8c3bbe5e1a, + 0x7eca04ebaf4a5eea, + 0x0543c37757f08d9a, + }; + + for (seq) |s| { + var buf0: [8]u8 = undefined; + var buf1: [7]u8 = undefined; + std.mem.writeInt(u64, &buf0, s, .little); + r.fill(&buf1); + try std.testing.expect(std.mem.eql(u8, buf0[0..7], buf1[0..])); + } +} diff --git a/lib/std/Random/benchmark.zig b/lib/std/Random/benchmark.zig @@ -0,0 +1,217 @@ +// zig run -O ReleaseFast --zig-lib-dir ../.. benchmark.zig + +const std = @import("std"); +const builtin = @import("builtin"); +const time = std.time; +const Timer = time.Timer; +const Random = std.Random; + +const KiB = 1024; +const MiB = 1024 * KiB; +const GiB = 1024 * MiB; + +const Rng = struct { + ty: type, + name: []const u8, + init_u8s: ?[]const u8 = null, + init_u64: ?u64 = null, +}; + +const prngs = [_]Rng{ + Rng{ + .ty = Random.Isaac64, + .name = "isaac64", + .init_u64 = 0, + }, + Rng{ + .ty = Random.Pcg, + .name = "pcg", + .init_u64 = 0, + }, + Rng{ + .ty = Random.RomuTrio, + .name = "romutrio", + .init_u64 = 0, + }, + Rng{ + .ty = Random.Sfc64, + .name = "sfc64", + .init_u64 = 0, + }, + Rng{ + .ty = Random.Xoroshiro128, + .name = "xoroshiro128", + .init_u64 = 0, + }, + Rng{ + .ty = Random.Xoshiro256, + .name = "xoshiro256", + .init_u64 = 0, + }, +}; + +const csprngs = [_]Rng{ + Rng{ + .ty = Random.Ascon, + .name = "ascon", + .init_u8s = &[_]u8{0} ** 32, + }, + Rng{ + .ty = Random.ChaCha, + .name = "chacha", + .init_u8s = &[_]u8{0} ** 32, + }, +}; + +const Result = struct { + throughput: u64, +}; + +const long_block_size: usize = 8 * 8192; +const short_block_size: usize = 8; + +pub fn benchmark(comptime H: anytype, bytes: usize, comptime block_size: usize) !Result { + var rng = blk: { + if (H.init_u8s) |init| { + break :blk H.ty.init(init[0..].*); + } + if (H.init_u64) |init| { + break :blk H.ty.init(init); + } + break :blk H.ty.init(); + }; + + var block: [block_size]u8 = undefined; + + var offset: usize = 0; + var timer = try Timer.start(); + const start = timer.lap(); + while (offset < bytes) : (offset += block.len) { + rng.fill(block[0..]); + } + const end = timer.read(); + + const elapsed_s = @as(f64, @floatFromInt(end - start)) / time.ns_per_s; + const throughput = @as(u64, @intFromFloat(@as(f64, @floatFromInt(bytes)) / elapsed_s)); + + std.debug.assert(rng.random().int(u64) != 0); + + return Result{ + .throughput = throughput, + }; +} + +fn usage() void { + std.debug.print( + \\throughput_test [options] + \\ + \\Options: + \\ --filter [test-name] + \\ --count [int] + \\ --prngs-only + \\ --csprngs-only + \\ --short-only + \\ --long-only + \\ --help + \\ + , .{}); +} + +fn mode(comptime x: comptime_int) comptime_int { + return if (builtin.mode == .Debug) x / 64 else x; +} + +pub fn main() !void { + const stdout = std.io.getStdOut().writer(); + + var buffer: [1024]u8 = undefined; + var fixed = std.heap.FixedBufferAllocator.init(buffer[0..]); + const args = try std.process.argsAlloc(fixed.allocator()); + + var filter: ?[]u8 = ""; + var count: usize = mode(128 * MiB); + var bench_prngs = true; + var bench_csprngs = true; + var bench_long = true; + var bench_short = true; + + var i: usize = 1; + while (i < args.len) : (i += 1) { + if (std.mem.eql(u8, args[i], "--mode")) { + try stdout.print("{}\n", .{builtin.mode}); + return; + } else if (std.mem.eql(u8, args[i], "--filter")) { + i += 1; + if (i == args.len) { + usage(); + std.os.exit(1); + } + + filter = args[i]; + } else if (std.mem.eql(u8, args[i], "--count")) { + i += 1; + if (i == args.len) { + usage(); + std.os.exit(1); + } + + const c = try std.fmt.parseUnsigned(usize, args[i], 10); + count = c * MiB; + } else if (std.mem.eql(u8, args[i], "--csprngs-only")) { + bench_prngs = false; + } else if (std.mem.eql(u8, args[i], "--prngs-only")) { + bench_csprngs = false; + } else if (std.mem.eql(u8, args[i], "--short-only")) { + bench_long = false; + } else if (std.mem.eql(u8, args[i], "--long-only")) { + bench_short = false; + } else if (std.mem.eql(u8, args[i], "--help")) { + usage(); + return; + } else { + usage(); + std.os.exit(1); + } + } + + if (bench_prngs) { + if (bench_long) { + inline for (prngs) |R| { + if (filter == null or std.mem.indexOf(u8, R.name, filter.?) != null) { + try stdout.print("{s} (long outputs)\n", .{R.name}); + const result_long = try benchmark(R, count, long_block_size); + try stdout.print(" {:5} MiB/s\n", .{result_long.throughput / (1 * MiB)}); + } + } + } + if (bench_short) { + inline for (prngs) |R| { + if (filter == null or std.mem.indexOf(u8, R.name, filter.?) != null) { + try stdout.print("{s} (short outputs)\n", .{R.name}); + const result_short = try benchmark(R, count, short_block_size); + try stdout.print(" {:5} MiB/s\n", .{result_short.throughput / (1 * MiB)}); + } + } + } + } + if (bench_csprngs) { + if (bench_long) { + inline for (csprngs) |R| { + if (filter == null or std.mem.indexOf(u8, R.name, filter.?) != null) { + try stdout.print("{s} (cryptographic, long outputs)\n", .{R.name}); + const result_long = try benchmark(R, count, long_block_size); + try stdout.print(" {:5} MiB/s\n", .{result_long.throughput / (1 * MiB)}); + } + } + } + if (bench_short) { + inline for (csprngs) |R| { + if (filter == null or std.mem.indexOf(u8, R.name, filter.?) != null) { + try stdout.print("{s} (cryptographic, short outputs)\n", .{R.name}); + const result_short = try benchmark(R, count, short_block_size); + try stdout.print(" {:5} MiB/s\n", .{result_short.throughput / (1 * MiB)}); + } + } + } + } +} diff --git a/lib/std/Random/test.zig b/lib/std/Random/test.zig @@ -0,0 +1,473 @@ +const std = @import("../std.zig"); +const math = std.math; +const Random = std.Random; +const DefaultPrng = Random.DefaultPrng; +const SplitMix64 = Random.SplitMix64; +const DefaultCsprng = Random.DefaultCsprng; +const expect = std.testing.expect; +const expectEqual = std.testing.expectEqual; + +const SequentialPrng = struct { + const Self = @This(); + next_value: u8, + + pub fn init() Self { + return Self{ + .next_value = 0, + }; + } + + pub fn random(self: *Self) Random { + return Random.init(self, fill); + } + + pub fn fill(self: *Self, buf: []u8) void { + for (buf) |*b| { + b.* = self.next_value; + } + self.next_value +%= 1; + } +}; + +/// Do not use this PRNG! It is meant to be predictable, for the purposes of test reproducibility and coverage. +/// Its output is just a repeat of a user-specified byte pattern. +/// Name is a reference to this comic: https://dilbert.com/strip/2001-10-25 +const Dilbert = struct { + pattern: []const u8 = undefined, + curr_idx: usize = 0, + + pub fn init(pattern: []const u8) !Dilbert { + if (pattern.len == 0) + return error.EmptyPattern; + var self = Dilbert{}; + self.pattern = pattern; + self.curr_idx = 0; + return self; + } + + pub fn random(self: *Dilbert) Random { + return Random.init(self, fill); + } + + pub fn fill(self: *Dilbert, buf: []u8) void { + for (buf) |*byte| { + byte.* = self.pattern[self.curr_idx]; + self.curr_idx = (self.curr_idx + 1) % self.pattern.len; + } + } + + test "Dilbert fill" { + var r = try Dilbert.init("9nine"); + + const seq = [_]u64{ + 0x396E696E65396E69, + 0x6E65396E696E6539, + 0x6E696E65396E696E, + 0x65396E696E65396E, + 0x696E65396E696E65, + }; + + for (seq) |s| { + var buf0: [8]u8 = undefined; + var buf1: [8]u8 = undefined; + std.mem.writeInt(u64, &buf0, s, .big); + r.fill(&buf1); + try std.testing.expect(std.mem.eql(u8, buf0[0..], buf1[0..])); + } + } +}; + +test "Random int" { + try testRandomInt(); + try comptime testRandomInt(); +} +fn testRandomInt() !void { + var rng = SequentialPrng.init(); + const random = rng.random(); + + try expect(random.int(u0) == 0); + + rng.next_value = 0; + try expect(random.int(u1) == 0); + try expect(random.int(u1) == 1); + try expect(random.int(u2) == 2); + try expect(random.int(u2) == 3); + try expect(random.int(u2) == 0); + + rng.next_value = 0xff; + try expect(random.int(u8) == 0xff); + rng.next_value = 0x11; + try expect(random.int(u8) == 0x11); + + rng.next_value = 0xff; + try expect(random.int(u32) == 0xffffffff); + rng.next_value = 0x11; + try expect(random.int(u32) == 0x11111111); + + rng.next_value = 0xff; + try expect(random.int(i32) == -1); + rng.next_value = 0x11; + try expect(random.int(i32) == 0x11111111); + + rng.next_value = 0xff; + try expect(random.int(i8) == -1); + rng.next_value = 0x11; + try expect(random.int(i8) == 0x11); + + rng.next_value = 0xff; + try expect(random.int(u33) == 0x1ffffffff); + rng.next_value = 0xff; + try expect(random.int(i1) == -1); + rng.next_value = 0xff; + try expect(random.int(i2) == -1); + rng.next_value = 0xff; + try expect(random.int(i33) == -1); +} + +test "Random boolean" { + try testRandomBoolean(); + try comptime testRandomBoolean(); +} +fn testRandomBoolean() !void { + var rng = SequentialPrng.init(); + const random = rng.random(); + + try expect(random.boolean() == false); + try expect(random.boolean() == true); + try expect(random.boolean() == false); + try expect(random.boolean() == true); +} + +test "Random enum" { + try testRandomEnumValue(); + try comptime testRandomEnumValue(); +} +fn testRandomEnumValue() !void { + const TestEnum = enum { + First, + Second, + Third, + }; + var rng = SequentialPrng.init(); + const random = rng.random(); + rng.next_value = 0; + try expect(random.enumValue(TestEnum) == TestEnum.First); + try expect(random.enumValue(TestEnum) == TestEnum.First); + try expect(random.enumValue(TestEnum) == TestEnum.First); +} + +test "Random intLessThan" { + @setEvalBranchQuota(10000); + try testRandomIntLessThan(); + try comptime testRandomIntLessThan(); +} +fn testRandomIntLessThan() !void { + var rng = SequentialPrng.init(); + const random = rng.random(); + + rng.next_value = 0xff; + try expect(random.uintLessThan(u8, 4) == 3); + try expect(rng.next_value == 0); + try expect(random.uintLessThan(u8, 4) == 0); + try expect(rng.next_value == 1); + + rng.next_value = 0; + try expect(random.uintLessThan(u64, 32) == 0); + + // trigger the bias rejection code path + rng.next_value = 0; + try expect(random.uintLessThan(u8, 3) == 0); + // verify we incremented twice + try expect(rng.next_value == 2); + + rng.next_value = 0xff; + try expect(random.intRangeLessThan(u8, 0, 0x80) == 0x7f); + rng.next_value = 0xff; + try expect(random.intRangeLessThan(u8, 0x7f, 0xff) == 0xfe); + + rng.next_value = 0xff; + try expect(random.intRangeLessThan(i8, 0, 0x40) == 0x3f); + rng.next_value = 0xff; + try expect(random.intRangeLessThan(i8, -0x40, 0x40) == 0x3f); + rng.next_value = 0xff; + try expect(random.intRangeLessThan(i8, -0x80, 0) == -1); + + rng.next_value = 0xff; + try expect(random.intRangeLessThan(i3, -4, 0) == -1); + rng.next_value = 0xff; + try expect(random.intRangeLessThan(i3, -2, 2) == 1); +} + +test "Random intAtMost" { + @setEvalBranchQuota(10000); + try testRandomIntAtMost(); + try comptime testRandomIntAtMost(); +} +fn testRandomIntAtMost() !void { + var rng = SequentialPrng.init(); + const random = rng.random(); + + rng.next_value = 0xff; + try expect(random.uintAtMost(u8, 3) == 3); + try expect(rng.next_value == 0); + try expect(random.uintAtMost(u8, 3) == 0); + + // trigger the bias rejection code path + rng.next_value = 0; + try expect(random.uintAtMost(u8, 2) == 0); + // verify we incremented twice + try expect(rng.next_value == 2); + + rng.next_value = 0xff; + try expect(random.intRangeAtMost(u8, 0, 0x7f) == 0x7f); + rng.next_value = 0xff; + try expect(random.intRangeAtMost(u8, 0x7f, 0xfe) == 0xfe); + + rng.next_value = 0xff; + try expect(random.intRangeAtMost(i8, 0, 0x3f) == 0x3f); + rng.next_value = 0xff; + try expect(random.intRangeAtMost(i8, -0x40, 0x3f) == 0x3f); + rng.next_value = 0xff; + try expect(random.intRangeAtMost(i8, -0x80, -1) == -1); + + rng.next_value = 0xff; + try expect(random.intRangeAtMost(i3, -4, -1) == -1); + rng.next_value = 0xff; + try expect(random.intRangeAtMost(i3, -2, 1) == 1); + + try expect(random.uintAtMost(u0, 0) == 0); +} + +test "Random Biased" { + var prng = DefaultPrng.init(0); + const random = prng.random(); + // Not thoroughly checking the logic here. + // Just want to execute all the paths with different types. + + try expect(random.uintLessThanBiased(u1, 1) == 0); + try expect(random.uintLessThanBiased(u32, 10) < 10); + try expect(random.uintLessThanBiased(u64, 20) < 20); + + try expect(random.uintAtMostBiased(u0, 0) == 0); + try expect(random.uintAtMostBiased(u1, 0) <= 0); + try expect(random.uintAtMostBiased(u32, 10) <= 10); + try expect(random.uintAtMostBiased(u64, 20) <= 20); + + try expect(random.intRangeLessThanBiased(u1, 0, 1) == 0); + try expect(random.intRangeLessThanBiased(i1, -1, 0) == -1); + try expect(random.intRangeLessThanBiased(u32, 10, 20) >= 10); + try expect(random.intRangeLessThanBiased(i32, 10, 20) >= 10); + try expect(random.intRangeLessThanBiased(u64, 20, 40) >= 20); + try expect(random.intRangeLessThanBiased(i64, 20, 40) >= 20); + + // uncomment for broken module error: + //expect(random.intRangeAtMostBiased(u0, 0, 0) == 0); + try expect(random.intRangeAtMostBiased(u1, 0, 1) >= 0); + try expect(random.intRangeAtMostBiased(i1, -1, 0) >= -1); + try expect(random.intRangeAtMostBiased(u32, 10, 20) >= 10); + try expect(random.intRangeAtMostBiased(i32, 10, 20) >= 10); + try expect(random.intRangeAtMostBiased(u64, 20, 40) >= 20); + try expect(random.intRangeAtMostBiased(i64, 20, 40) >= 20); +} + +test "splitmix64 sequence" { + var r = SplitMix64.init(0xaeecf86f7878dd75); + + const seq = [_]u64{ + 0x5dbd39db0178eb44, + 0xa9900fb66b397da3, + 0x5c1a28b1aeebcf5c, + 0x64a963238f776912, + 0xc6d4177b21d1c0ab, + 0xb2cbdbdb5ea35394, + }; + + for (seq) |s| { + try expect(s == r.next()); + } +} + +// Actual Random helper function tests, pcg engine is assumed correct. +test "Random float correctness" { + var prng = DefaultPrng.init(0); + const random = prng.random(); + + var i: usize = 0; + while (i < 1000) : (i += 1) { + const val1 = random.float(f32); + try expect(val1 >= 0.0); + try expect(val1 < 1.0); + + const val2 = random.float(f64); + try expect(val2 >= 0.0); + try expect(val2 < 1.0); + } +} + +// Check the "astronomically unlikely" code paths. +test "Random float coverage" { + var prng = try Dilbert.init(&[_]u8{0}); + const random = prng.random(); + + const rand_f64 = random.float(f64); + const rand_f32 = random.float(f32); + + try expect(rand_f32 == 0.0); + try expect(rand_f64 == 0.0); +} + +test "Random float chi-square goodness of fit" { + const num_numbers = 100000; + const num_buckets = 1000; + + var f32_hist = std.AutoHashMap(u32, u32).init(std.testing.allocator); + defer f32_hist.deinit(); + var f64_hist = std.AutoHashMap(u64, u32).init(std.testing.allocator); + defer f64_hist.deinit(); + + var prng = DefaultPrng.init(0); + const random = prng.random(); + + var i: usize = 0; + while (i < num_numbers) : (i += 1) { + const rand_f32 = random.float(f32); + const rand_f64 = random.float(f64); + const f32_put = try f32_hist.getOrPut(@as(u32, @intFromFloat(rand_f32 * @as(f32, @floatFromInt(num_buckets))))); + if (f32_put.found_existing) { + f32_put.value_ptr.* += 1; + } else { + f32_put.value_ptr.* = 1; + } + const f64_put = try f64_hist.getOrPut(@as(u32, @intFromFloat(rand_f64 * @as(f64, @floatFromInt(num_buckets))))); + if (f64_put.found_existing) { + f64_put.value_ptr.* += 1; + } else { + f64_put.value_ptr.* = 1; + } + } + + var f32_total_variance: f64 = 0; + var f64_total_variance: f64 = 0; + + { + var j: u32 = 0; + while (j < num_buckets) : (j += 1) { + const count = @as(f64, @floatFromInt((if (f32_hist.get(j)) |v| v else 0))); + const expected = @as(f64, @floatFromInt(num_numbers)) / @as(f64, @floatFromInt(num_buckets)); + const delta = count - expected; + const variance = (delta * delta) / expected; + f32_total_variance += variance; + } + } + + { + var j: u64 = 0; + while (j < num_buckets) : (j += 1) { + const count = @as(f64, @floatFromInt((if (f64_hist.get(j)) |v| v else 0))); + const expected = @as(f64, @floatFromInt(num_numbers)) / @as(f64, @floatFromInt(num_buckets)); + const delta = count - expected; + const variance = (delta * delta) / expected; + f64_total_variance += variance; + } + } + + // Accept p-values >= 0.05. + // Critical value is calculated by opening a Python interpreter and running: + // scipy.stats.chi2.isf(0.05, num_buckets - 1) + const critical_value = 1073.6426506574246; + try expect(f32_total_variance < critical_value); + try expect(f64_total_variance < critical_value); +} + +test "Random shuffle" { + var prng = DefaultPrng.init(0); + const random = prng.random(); + + var seq = [_]u8{ 0, 1, 2, 3, 4 }; + var seen = [_]bool{false} ** 5; + + var i: usize = 0; + while (i < 1000) : (i += 1) { + random.shuffle(u8, seq[0..]); + seen[seq[0]] = true; + try expect(sumArray(seq[0..]) == 10); + } + + // we should see every entry at the head at least once + for (seen) |e| { + try expect(e == true); + } +} + +fn sumArray(s: []const u8) u32 { + var r: u32 = 0; + for (s) |e| + r += e; + return r; +} + +test "Random range" { + var prng = DefaultPrng.init(0); + const random = prng.random(); + + try testRange(random, -4, 3); + try testRange(random, -4, -1); + try testRange(random, 10, 14); + try testRange(random, -0x80, 0x7f); +} + +fn testRange(r: Random, start: i8, end: i8) !void { + try testRangeBias(r, start, end, true); + try testRangeBias(r, start, end, false); +} +fn testRangeBias(r: Random, start: i8, end: i8, biased: bool) !void { + const count = @as(usize, @intCast(@as(i32, end) - @as(i32, start))); + var values_buffer = [_]bool{false} ** 0x100; + const values = values_buffer[0..count]; + var i: usize = 0; + while (i < count) { + const value: i32 = if (biased) r.intRangeLessThanBiased(i8, start, end) else r.intRangeLessThan(i8, start, end); + const index = @as(usize, @intCast(value - start)); + if (!values[index]) { + i += 1; + values[index] = true; + } + } +} + +test "CSPRNG" { + var secret_seed: [DefaultCsprng.secret_seed_length]u8 = undefined; + std.crypto.random.bytes(&secret_seed); + var csprng = DefaultCsprng.init(secret_seed); + const random = csprng.random(); + const a = random.int(u64); + const b = random.int(u64); + const c = random.int(u64); + try expect(a ^ b ^ c != 0); +} + +test "Random weightedIndex" { + // Make sure weightedIndex works for various integers and floats + inline for (.{ u64, i4, f32, f64 }) |T| { + var prng = DefaultPrng.init(0); + const random = prng.random(); + + const proportions = [_]T{ 2, 1, 1, 2 }; + var counts = [_]f64{ 0, 0, 0, 0 }; + + const n_trials: u64 = 10_000; + var i: usize = 0; + while (i < n_trials) : (i += 1) { + const pick = random.weightedIndex(T, &proportions); + counts[pick] += 1; + } + + // We expect the first and last counts to be roughly 2x the second and third + const approxEqRel = std.math.approxEqRel; + // Define "roughly" to be within 10% + const tolerance = 0.1; + try std.testing.expect(approxEqRel(f64, counts[0], counts[1] * 2, tolerance)); + try std.testing.expect(approxEqRel(f64, counts[1], counts[2], tolerance)); + try std.testing.expect(approxEqRel(f64, counts[2] * 2, counts[3], tolerance)); + } +} diff --git a/lib/std/Random/ziggurat.zig b/lib/std/Random/ziggurat.zig @@ -0,0 +1,170 @@ +//! Implements [ZIGNOR][1] (Jurgen A. Doornik, 2005, Nuffield College, Oxford). +//! +//! [1]: https://www.doornik.com/research/ziggurat.pdf +//! +//! rust/rand used as a reference; +//! +//! NOTE: This seems interesting but reference code is a bit hard to grok: +//! https://sbarral.github.io/etf. + +const std = @import("../std.zig"); +const builtin = @import("builtin"); +const math = std.math; +const Random = std.Random; + +pub fn next_f64(random: Random, comptime tables: ZigTable) f64 { + while (true) { + // We manually construct a float from parts as we can avoid an extra random lookup here by + // using the unused exponent for the lookup table entry. + const bits = random.int(u64); + const i = @as(usize, @as(u8, @truncate(bits))); + + const u = blk: { + if (tables.is_symmetric) { + // Generate a value in the range [2, 4) and scale into [-1, 1) + const repr = ((0x3ff + 1) << 52) | (bits >> 12); + break :blk @as(f64, @bitCast(repr)) - 3.0; + } else { + // Generate a value in the range [1, 2) and scale into (0, 1) + const repr = (0x3ff << 52) | (bits >> 12); + break :blk @as(f64, @bitCast(repr)) - (1.0 - math.floatEps(f64) / 2.0); + } + }; + + const x = u * tables.x[i]; + const test_x = if (tables.is_symmetric) @abs(x) else x; + + // equivalent to |u| < tables.x[i+1] / tables.x[i] (or u < tables.x[i+1] / tables.x[i]) + if (test_x < tables.x[i + 1]) { + return x; + } + + if (i == 0) { + return tables.zero_case(random, u); + } + + // equivalent to f1 + DRanU() * (f0 - f1) < 1 + if (tables.f[i + 1] + (tables.f[i] - tables.f[i + 1]) * random.float(f64) < tables.pdf(x)) { + return x; + } + } +} + +pub const ZigTable = struct { + r: f64, + x: [257]f64, + f: [257]f64, + + // probability density function used as a fallback + pdf: fn (f64) f64, + // whether the distribution is symmetric + is_symmetric: bool, + // fallback calculation in the case we are in the 0 block + zero_case: fn (Random, f64) f64, +}; + +// zigNorInit +pub fn ZigTableGen( + comptime is_symmetric: bool, + comptime r: f64, + comptime v: f64, + comptime f: fn (f64) f64, + comptime f_inv: fn (f64) f64, + comptime zero_case: fn (Random, f64) f64, +) ZigTable { + var tables: ZigTable = undefined; + + tables.is_symmetric = is_symmetric; + tables.r = r; + tables.pdf = f; + tables.zero_case = zero_case; + + tables.x[0] = v / f(r); + tables.x[1] = r; + + for (tables.x[2..256], 0..) |*entry, i| { + const last = tables.x[2 + i - 1]; + entry.* = f_inv(v / last + f(last)); + } + tables.x[256] = 0; + + for (tables.f[0..], 0..) |*entry, i| { + entry.* = f(tables.x[i]); + } + + return tables; +} + +// N(0, 1) +pub const NormDist = blk: { + @setEvalBranchQuota(30000); + break :blk ZigTableGen(true, norm_r, norm_v, norm_f, norm_f_inv, norm_zero_case); +}; + +pub const norm_r = 3.6541528853610088; +pub const norm_v = 0.00492867323399; + +pub fn norm_f(x: f64) f64 { + return @exp(-x * x / 2.0); +} +pub fn norm_f_inv(y: f64) f64 { + return @sqrt(-2.0 * @log(y)); +} +pub fn norm_zero_case(random: Random, u: f64) f64 { + var x: f64 = 1; + var y: f64 = 0; + + while (-2.0 * y < x * x) { + x = @log(random.float(f64)) / norm_r; + y = @log(random.float(f64)); + } + + if (u < 0) { + return x - norm_r; + } else { + return norm_r - x; + } +} + +test "normal dist sanity" { + var prng = Random.DefaultPrng.init(0); + const random = prng.random(); + + var i: usize = 0; + while (i < 1000) : (i += 1) { + _ = random.floatNorm(f64); + } +} + +// Exp(1) +pub const ExpDist = blk: { + @setEvalBranchQuota(30000); + break :blk ZigTableGen(false, exp_r, exp_v, exp_f, exp_f_inv, exp_zero_case); +}; + +pub const exp_r = 7.69711747013104972; +pub const exp_v = 0.0039496598225815571993; + +pub fn exp_f(x: f64) f64 { + return @exp(-x); +} +pub fn exp_f_inv(y: f64) f64 { + return -@log(y); +} +pub fn exp_zero_case(random: Random, _: f64) f64 { + return exp_r - @log(random.float(f64)); +} + +test "exp dist smoke test" { + var prng = Random.DefaultPrng.init(0); + const random = prng.random(); + + var i: usize = 0; + while (i < 1000) : (i += 1) { + _ = random.floatExp(f64); + } +} + +test { + _ = NormDist; +} diff --git a/lib/std/Thread/RwLock.zig b/lib/std/Thread/RwLock.zig @@ -328,7 +328,7 @@ test "RwLock - concurrent access" { } fn writer(self: *Self, thread_idx: usize) !void { - var prng = std.rand.DefaultPrng.init(thread_idx); + var prng = std.Random.DefaultPrng.init(thread_idx); var rnd = prng.random(); while (true) { diff --git a/lib/std/crypto/benchmark.zig b/lib/std/crypto/benchmark.zig @@ -10,7 +10,7 @@ const crypto = std.crypto; const KiB = 1024; const MiB = 1024 * KiB; -var prng = std.rand.DefaultPrng.init(0); +var prng = std.Random.DefaultPrng.init(0); const random = prng.random(); const Crypto = struct { diff --git a/lib/std/crypto/kyber_d00.zig b/lib/std/crypto/kyber_d00.zig @@ -110,7 +110,7 @@ const assert = std.debug.assert; const crypto = std.crypto; const math = std.math; const mem = std.mem; -const RndGen = std.rand.DefaultPrng; +const RndGen = std.Random.DefaultPrng; const sha3 = crypto.hash.sha3; // Q is the parameter q ≡ 3329 = 2¹¹ + 2¹⁰ + 2⁸ + 1. diff --git a/lib/std/crypto/tlcsprng.zig b/lib/std/crypto/tlcsprng.zig @@ -10,7 +10,7 @@ const os = std.os; /// We use this as a layer of indirection because global const pointers cannot /// point to thread-local variables. -pub const interface = std.rand.Random{ +pub const interface = std.Random{ .ptr = undefined, .fillFn = tlsCsprngFill, }; @@ -43,7 +43,7 @@ const maybe_have_wipe_on_fork = builtin.os.isAtLeast(.linux, .{ }) orelse true; const is_haiku = builtin.os.tag == .haiku; -const Rng = std.rand.DefaultCsprng; +const Rng = std.Random.DefaultCsprng; const Context = struct { init_state: enum(u8) { uninitialized = 0, initialized, failed }, diff --git a/lib/std/hash/benchmark.zig b/lib/std/hash/benchmark.zig @@ -10,7 +10,7 @@ const KiB = 1024; const MiB = 1024 * KiB; const GiB = 1024 * MiB; -var prng = std.rand.DefaultPrng.init(0); +var prng = std.Random.DefaultPrng.init(0); const random = prng.random(); const Hash = struct { diff --git a/lib/std/hash_map.zig b/lib/std/hash_map.zig @@ -1884,7 +1884,7 @@ test "std.hash_map put and remove loop in random order" { while (i < size) : (i += 1) { try keys.append(i); } - var prng = std.rand.DefaultPrng.init(0); + var prng = std.Random.DefaultPrng.init(0); const random = prng.random(); while (i < iterations) : (i += 1) { @@ -1916,7 +1916,7 @@ test "std.hash_map remove one million elements in random order" { keys.append(i) catch unreachable; } - var prng = std.rand.DefaultPrng.init(0); + var prng = std.Random.DefaultPrng.init(0); const random = prng.random(); random.shuffle(u32, keys.items); diff --git a/lib/std/heap/arena_allocator.zig b/lib/std/heap/arena_allocator.zig @@ -250,7 +250,7 @@ test "ArenaAllocator (reset with preheating)" { var arena_allocator = ArenaAllocator.init(std.testing.allocator); defer arena_allocator.deinit(); // provides some variance in the allocated data - var rng_src = std.rand.DefaultPrng.init(19930913); + var rng_src = std.Random.DefaultPrng.init(19930913); const random = rng_src.random(); var rounds: usize = 25; while (rounds > 0) { diff --git a/lib/std/io/test.zig b/lib/std/io/test.zig @@ -1,6 +1,6 @@ const std = @import("std"); const io = std.io; -const DefaultPrng = std.rand.DefaultPrng; +const DefaultPrng = std.Random.DefaultPrng; const expect = std.testing.expect; const expectEqual = std.testing.expectEqual; const expectError = std.testing.expectError; diff --git a/lib/std/math/big/rational.zig b/lib/std/math/big/rational.zig @@ -594,7 +594,7 @@ test "big.rational toFloat" { test "big.rational set/to Float round-trip" { var a = try Rational.init(testing.allocator); defer a.deinit(); - var prng = std.rand.DefaultPrng.init(0x5EED); + var prng = std.Random.DefaultPrng.init(0x5EED); const random = prng.random(); var i: usize = 0; while (i < 512) : (i += 1) { diff --git a/lib/std/mem.zig b/lib/std/mem.zig @@ -4423,7 +4423,7 @@ test "read/write(Var)PackedInt" { const foreign_endian: Endian = if (native_endian == .big) .little else .big; const expect = std.testing.expect; - var prng = std.rand.DefaultPrng.init(1234); + var prng = std.Random.DefaultPrng.init(1234); const random = prng.random(); @setEvalBranchQuota(10_000); diff --git a/lib/std/priority_dequeue.zig b/lib/std/priority_dequeue.zig @@ -866,7 +866,7 @@ test "std.PriorityDequeue: shrinkAndFree" { } test "std.PriorityDequeue: fuzz testing min" { - var prng = std.rand.DefaultPrng.init(0x12345678); + var prng = std.Random.DefaultPrng.init(0x12345678); const random = prng.random(); const test_case_count = 100; @@ -878,7 +878,7 @@ test "std.PriorityDequeue: fuzz testing min" { } } -fn fuzzTestMin(rng: std.rand.Random, comptime queue_size: usize) !void { +fn fuzzTestMin(rng: std.Random, comptime queue_size: usize) !void { const allocator = testing.allocator; const items = try generateRandomSlice(allocator, rng, queue_size); @@ -895,7 +895,7 @@ fn fuzzTestMin(rng: std.rand.Random, comptime queue_size: usize) !void { } test "std.PriorityDequeue: fuzz testing max" { - var prng = std.rand.DefaultPrng.init(0x87654321); + var prng = std.Random.DefaultPrng.init(0x87654321); const random = prng.random(); const test_case_count = 100; @@ -907,7 +907,7 @@ test "std.PriorityDequeue: fuzz testing max" { } } -fn fuzzTestMax(rng: std.rand.Random, queue_size: usize) !void { +fn fuzzTestMax(rng: std.Random, queue_size: usize) !void { const allocator = testing.allocator; const items = try generateRandomSlice(allocator, rng, queue_size); @@ -924,7 +924,7 @@ fn fuzzTestMax(rng: std.rand.Random, queue_size: usize) !void { } test "std.PriorityDequeue: fuzz testing min and max" { - var prng = std.rand.DefaultPrng.init(0x87654321); + var prng = std.Random.DefaultPrng.init(0x87654321); const random = prng.random(); const test_case_count = 100; @@ -936,7 +936,7 @@ test "std.PriorityDequeue: fuzz testing min and max" { } } -fn fuzzTestMinMax(rng: std.rand.Random, queue_size: usize) !void { +fn fuzzTestMinMax(rng: std.Random, queue_size: usize) !void { const allocator = testing.allocator; const items = try generateRandomSlice(allocator, rng, queue_size); @@ -963,7 +963,7 @@ fn fuzzTestMinMax(rng: std.rand.Random, queue_size: usize) !void { } } -fn generateRandomSlice(allocator: std.mem.Allocator, rng: std.rand.Random, size: usize) ![]u32 { +fn generateRandomSlice(allocator: std.mem.Allocator, rng: std.Random, size: usize) ![]u32 { var array = std.ArrayList(u32).init(allocator); try array.ensureTotalCapacity(size); diff --git a/lib/std/rand.zig b/lib/std/rand.zig @@ -1,460 +0,0 @@ -//! The engines provided here should be initialized from an external source. -//! For a thread-local cryptographically secure pseudo random number generator, -//! use `std.crypto.random`. -//! Be sure to use a CSPRNG when required, otherwise using a normal PRNG will -//! be faster and use substantially less stack space. - -const std = @import("std.zig"); -const builtin = @import("builtin"); -const assert = std.debug.assert; -const mem = std.mem; -const math = std.math; -const maxInt = std.math.maxInt; - -/// Fast unbiased random numbers. -pub const DefaultPrng = Xoshiro256; - -/// Cryptographically secure random numbers. -pub const DefaultCsprng = ChaCha; - -pub const Ascon = @import("rand/Ascon.zig"); -pub const ChaCha = @import("rand/ChaCha.zig"); - -pub const Isaac64 = @import("rand/Isaac64.zig"); -pub const Pcg = @import("rand/Pcg.zig"); -pub const Xoroshiro128 = @import("rand/Xoroshiro128.zig"); -pub const Xoshiro256 = @import("rand/Xoshiro256.zig"); -pub const Sfc64 = @import("rand/Sfc64.zig"); -pub const RomuTrio = @import("rand/RomuTrio.zig"); -pub const ziggurat = @import("rand/ziggurat.zig"); - -pub const Random = struct { - ptr: *anyopaque, - fillFn: *const fn (ptr: *anyopaque, buf: []u8) void, - - pub fn init(pointer: anytype, comptime fillFn: fn (ptr: @TypeOf(pointer), buf: []u8) void) Random { - const Ptr = @TypeOf(pointer); - assert(@typeInfo(Ptr) == .Pointer); // Must be a pointer - assert(@typeInfo(Ptr).Pointer.size == .One); // Must be a single-item pointer - assert(@typeInfo(@typeInfo(Ptr).Pointer.child) == .Struct); // Must point to a struct - const gen = struct { - fn fill(ptr: *anyopaque, buf: []u8) void { - const self: Ptr = @ptrCast(@alignCast(ptr)); - fillFn(self, buf); - } - }; - - return .{ - .ptr = pointer, - .fillFn = gen.fill, - }; - } - - /// Read random bytes into the specified buffer until full. - pub fn bytes(r: Random, buf: []u8) void { - r.fillFn(r.ptr, buf); - } - - pub fn boolean(r: Random) bool { - return r.int(u1) != 0; - } - - /// Returns a random value from an enum, evenly distributed. - /// - /// Note that this will not yield consistent results across all targets - /// due to dependence on the representation of `usize` as an index. - /// See `enumValueWithIndex` for further commentary. - pub inline fn enumValue(r: Random, comptime EnumType: type) EnumType { - return r.enumValueWithIndex(EnumType, usize); - } - - /// Returns a random value from an enum, evenly distributed. - /// - /// An index into an array of all named values is generated using the - /// specified `Index` type to determine the return value. - /// This allows for results to be independent of `usize` representation. - /// - /// Prefer `enumValue` if this isn't important. - /// - /// See `uintLessThan`, which this function uses in most cases, - /// for commentary on the runtime of this function. - pub fn enumValueWithIndex(r: Random, comptime EnumType: type, comptime Index: type) EnumType { - comptime assert(@typeInfo(EnumType) == .Enum); - - // We won't use int -> enum casting because enum elements can have - // arbitrary values. Instead we'll randomly pick one of the type's values. - const values = comptime std.enums.values(EnumType); - comptime assert(values.len > 0); // can't return anything - comptime assert(maxInt(Index) >= values.len - 1); // can't access all values - comptime if (values.len == 1) return values[0]; - - const index = if (comptime values.len - 1 == maxInt(Index)) - r.int(Index) - else - r.uintLessThan(Index, values.len); - - const MinInt = MinArrayIndex(Index); - return values[@as(MinInt, @intCast(index))]; - } - - /// Returns a random int `i` such that `minInt(T) <= i <= maxInt(T)`. - /// `i` is evenly distributed. - pub fn int(r: Random, comptime T: type) T { - const bits = @typeInfo(T).Int.bits; - const UnsignedT = std.meta.Int(.unsigned, bits); - const ceil_bytes = comptime std.math.divCeil(u16, bits, 8) catch unreachable; - const ByteAlignedT = std.meta.Int(.unsigned, ceil_bytes * 8); - - var rand_bytes: [ceil_bytes]u8 = undefined; - r.bytes(&rand_bytes); - - // use LE instead of native endian for better portability maybe? - // TODO: endian portability is pointless if the underlying prng isn't endian portable. - // TODO: document the endian portability of this library. - const byte_aligned_result = mem.readInt(ByteAlignedT, &rand_bytes, .little); - const unsigned_result: UnsignedT = @truncate(byte_aligned_result); - return @bitCast(unsigned_result); - } - - /// Constant-time implementation off `uintLessThan`. - /// The results of this function may be biased. - pub fn uintLessThanBiased(r: Random, comptime T: type, less_than: T) T { - comptime assert(@typeInfo(T).Int.signedness == .unsigned); - assert(0 < less_than); - return limitRangeBiased(T, r.int(T), less_than); - } - - /// Returns an evenly distributed random unsigned integer `0 <= i < less_than`. - /// This function assumes that the underlying `fillFn` produces evenly distributed values. - /// Within this assumption, the runtime of this function is exponentially distributed. - /// If `fillFn` were backed by a true random generator, - /// the runtime of this function would technically be unbounded. - /// However, if `fillFn` is backed by any evenly distributed pseudo random number generator, - /// this function is guaranteed to return. - /// If you need deterministic runtime bounds, use `uintLessThanBiased`. - pub fn uintLessThan(r: Random, comptime T: type, less_than: T) T { - comptime assert(@typeInfo(T).Int.signedness == .unsigned); - const bits = @typeInfo(T).Int.bits; - assert(0 < less_than); - - // adapted from: - // http://www.pcg-random.org/posts/bounded-rands.html - // "Lemire's (with an extra tweak from me)" - var x = r.int(T); - var m = math.mulWide(T, x, less_than); - var l: T = @truncate(m); - if (l < less_than) { - var t = -%less_than; - - if (t >= less_than) { - t -= less_than; - if (t >= less_than) { - t %= less_than; - } - } - while (l < t) { - x = r.int(T); - m = math.mulWide(T, x, less_than); - l = @truncate(m); - } - } - return @intCast(m >> bits); - } - - /// Constant-time implementation off `uintAtMost`. - /// The results of this function may be biased. - pub fn uintAtMostBiased(r: Random, comptime T: type, at_most: T) T { - assert(@typeInfo(T).Int.signedness == .unsigned); - if (at_most == maxInt(T)) { - // have the full range - return r.int(T); - } - return r.uintLessThanBiased(T, at_most + 1); - } - - /// Returns an evenly distributed random unsigned integer `0 <= i <= at_most`. - /// See `uintLessThan`, which this function uses in most cases, - /// for commentary on the runtime of this function. - pub fn uintAtMost(r: Random, comptime T: type, at_most: T) T { - assert(@typeInfo(T).Int.signedness == .unsigned); - if (at_most == maxInt(T)) { - // have the full range - return r.int(T); - } - return r.uintLessThan(T, at_most + 1); - } - - /// Constant-time implementation off `intRangeLessThan`. - /// The results of this function may be biased. - pub fn intRangeLessThanBiased(r: Random, comptime T: type, at_least: T, less_than: T) T { - assert(at_least < less_than); - const info = @typeInfo(T).Int; - if (info.signedness == .signed) { - // Two's complement makes this math pretty easy. - const UnsignedT = std.meta.Int(.unsigned, info.bits); - const lo: UnsignedT = @bitCast(at_least); - const hi: UnsignedT = @bitCast(less_than); - const result = lo +% r.uintLessThanBiased(UnsignedT, hi -% lo); - return @bitCast(result); - } else { - // The signed implementation would work fine, but we can use stricter arithmetic operators here. - return at_least + r.uintLessThanBiased(T, less_than - at_least); - } - } - - /// Returns an evenly distributed random integer `at_least <= i < less_than`. - /// See `uintLessThan`, which this function uses in most cases, - /// for commentary on the runtime of this function. - pub fn intRangeLessThan(r: Random, comptime T: type, at_least: T, less_than: T) T { - assert(at_least < less_than); - const info = @typeInfo(T).Int; - if (info.signedness == .signed) { - // Two's complement makes this math pretty easy. - const UnsignedT = std.meta.Int(.unsigned, info.bits); - const lo: UnsignedT = @bitCast(at_least); - const hi: UnsignedT = @bitCast(less_than); - const result = lo +% r.uintLessThan(UnsignedT, hi -% lo); - return @bitCast(result); - } else { - // The signed implementation would work fine, but we can use stricter arithmetic operators here. - return at_least + r.uintLessThan(T, less_than - at_least); - } - } - - /// Constant-time implementation off `intRangeAtMostBiased`. - /// The results of this function may be biased. - pub fn intRangeAtMostBiased(r: Random, comptime T: type, at_least: T, at_most: T) T { - assert(at_least <= at_most); - const info = @typeInfo(T).Int; - if (info.signedness == .signed) { - // Two's complement makes this math pretty easy. - const UnsignedT = std.meta.Int(.unsigned, info.bits); - const lo: UnsignedT = @bitCast(at_least); - const hi: UnsignedT = @bitCast(at_most); - const result = lo +% r.uintAtMostBiased(UnsignedT, hi -% lo); - return @bitCast(result); - } else { - // The signed implementation would work fine, but we can use stricter arithmetic operators here. - return at_least + r.uintAtMostBiased(T, at_most - at_least); - } - } - - /// Returns an evenly distributed random integer `at_least <= i <= at_most`. - /// See `uintLessThan`, which this function uses in most cases, - /// for commentary on the runtime of this function. - pub fn intRangeAtMost(r: Random, comptime T: type, at_least: T, at_most: T) T { - assert(at_least <= at_most); - const info = @typeInfo(T).Int; - if (info.signedness == .signed) { - // Two's complement makes this math pretty easy. - const UnsignedT = std.meta.Int(.unsigned, info.bits); - const lo: UnsignedT = @bitCast(at_least); - const hi: UnsignedT = @bitCast(at_most); - const result = lo +% r.uintAtMost(UnsignedT, hi -% lo); - return @bitCast(result); - } else { - // The signed implementation would work fine, but we can use stricter arithmetic operators here. - return at_least + r.uintAtMost(T, at_most - at_least); - } - } - - /// Return a floating point value evenly distributed in the range [0, 1). - pub fn float(r: Random, comptime T: type) T { - // Generate a uniformly random value for the mantissa. - // Then generate an exponentially biased random value for the exponent. - // This covers every possible value in the range. - switch (T) { - f32 => { - // Use 23 random bits for the mantissa, and the rest for the exponent. - // If all 41 bits are zero, generate additional random bits, until a - // set bit is found, or 126 bits have been generated. - const rand = r.int(u64); - var rand_lz = @clz(rand); - if (rand_lz >= 41) { - // TODO: when #5177 or #489 is implemented, - // tell the compiler it is unlikely (1/2^41) to reach this point. - // (Same for the if branch and the f64 calculations below.) - rand_lz = 41 + @clz(r.int(u64)); - if (rand_lz == 41 + 64) { - // It is astronomically unlikely to reach this point. - rand_lz += @clz(r.int(u32) | 0x7FF); - } - } - const mantissa: u23 = @truncate(rand); - const exponent = @as(u32, 126 - rand_lz) << 23; - return @bitCast(exponent | mantissa); - }, - f64 => { - // Use 52 random bits for the mantissa, and the rest for the exponent. - // If all 12 bits are zero, generate additional random bits, until a - // set bit is found, or 1022 bits have been generated. - const rand = r.int(u64); - var rand_lz: u64 = @clz(rand); - if (rand_lz >= 12) { - rand_lz = 12; - while (true) { - // It is astronomically unlikely for this loop to execute more than once. - const addl_rand_lz = @clz(r.int(u64)); - rand_lz += addl_rand_lz; - if (addl_rand_lz != 64) { - break; - } - if (rand_lz >= 1022) { - rand_lz = 1022; - break; - } - } - } - const mantissa = rand & 0xFFFFFFFFFFFFF; - const exponent = (1022 - rand_lz) << 52; - return @bitCast(exponent | mantissa); - }, - else => @compileError("unknown floating point type"), - } - } - - /// Return a floating point value normally distributed with mean = 0, stddev = 1. - /// - /// To use different parameters, use: floatNorm(...) * desiredStddev + desiredMean. - pub fn floatNorm(r: Random, comptime T: type) T { - const value = ziggurat.next_f64(r, ziggurat.NormDist); - switch (T) { - f32 => return @floatCast(value), - f64 => return value, - else => @compileError("unknown floating point type"), - } - } - - /// Return an exponentially distributed float with a rate parameter of 1. - /// - /// To use a different rate parameter, use: floatExp(...) / desiredRate. - pub fn floatExp(r: Random, comptime T: type) T { - const value = ziggurat.next_f64(r, ziggurat.ExpDist); - switch (T) { - f32 => return @floatCast(value), - f64 => return value, - else => @compileError("unknown floating point type"), - } - } - - /// Shuffle a slice into a random order. - /// - /// Note that this will not yield consistent results across all targets - /// due to dependence on the representation of `usize` as an index. - /// See `shuffleWithIndex` for further commentary. - pub inline fn shuffle(r: Random, comptime T: type, buf: []T) void { - r.shuffleWithIndex(T, buf, usize); - } - - /// Shuffle a slice into a random order, using an index of a - /// specified type to maintain distribution across targets. - /// Asserts the index type can represent `buf.len`. - /// - /// Indexes into the slice are generated using the specified `Index` - /// type, which determines distribution properties. This allows for - /// results to be independent of `usize` representation. - /// - /// Prefer `shuffle` if this isn't important. - /// - /// See `intRangeLessThan`, which this function uses, - /// for commentary on the runtime of this function. - pub fn shuffleWithIndex(r: Random, comptime T: type, buf: []T, comptime Index: type) void { - const MinInt = MinArrayIndex(Index); - if (buf.len < 2) { - return; - } - - // `i <= j < max <= maxInt(MinInt)` - const max: MinInt = @intCast(buf.len); - var i: MinInt = 0; - while (i < max - 1) : (i += 1) { - const j: MinInt = @intCast(r.intRangeLessThan(Index, i, max)); - mem.swap(T, &buf[i], &buf[j]); - } - } - - /// Randomly selects an index into `proportions`, where the likelihood of each - /// index is weighted by that proportion. - /// It is more likely for the index of the last proportion to be returned - /// than the index of the first proportion in the slice, and vice versa. - /// - /// This is useful for selecting an item from a slice where weights are not equal. - /// `T` must be a numeric type capable of holding the sum of `proportions`. - pub fn weightedIndex(r: std.rand.Random, comptime T: type, proportions: []const T) usize { - // This implementation works by summing the proportions and picking a - // random point in [0, sum). We then loop over the proportions, - // accumulating until our accumulator is greater than the random point. - - const sum = s: { - var sum: T = 0; - for (proportions) |v| sum += v; - break :s sum; - }; - - const point = switch (@typeInfo(T)) { - .Int => |int_info| switch (int_info.signedness) { - .signed => r.intRangeLessThan(T, 0, sum), - .unsigned => r.uintLessThan(T, sum), - }, - // take care that imprecision doesn't lead to a value slightly greater than sum - .Float => @min(r.float(T) * sum, sum - std.math.floatEps(T)), - else => @compileError("weightedIndex does not support proportions of type " ++ - @typeName(T)), - }; - - assert(point < sum); - - var accumulator: T = 0; - for (proportions, 0..) |p, index| { - accumulator += p; - if (point < accumulator) return index; - } else unreachable; - } - - /// Returns the smallest of `Index` and `usize`. - fn MinArrayIndex(comptime Index: type) type { - const index_info = @typeInfo(Index).Int; - assert(index_info.signedness == .unsigned); - return if (index_info.bits >= @typeInfo(usize).Int.bits) usize else Index; - } -}; - -/// Convert a random integer 0 <= random_int <= maxValue(T), -/// into an integer 0 <= result < less_than. -/// This function introduces a minor bias. -pub fn limitRangeBiased(comptime T: type, random_int: T, less_than: T) T { - comptime assert(@typeInfo(T).Int.signedness == .unsigned); - const bits = @typeInfo(T).Int.bits; - - // adapted from: - // http://www.pcg-random.org/posts/bounded-rands.html - // "Integer Multiplication (Biased)" - const m = math.mulWide(T, random_int, less_than); - return @intCast(m >> bits); -} - -// Generator to extend 64-bit seed values into longer sequences. -// -// The number of cycles is thus limited to 64-bits regardless of the engine, but this -// is still plenty for practical purposes. -pub const SplitMix64 = struct { - s: u64, - - pub fn init(seed: u64) SplitMix64 { - return SplitMix64{ .s = seed }; - } - - pub fn next(self: *SplitMix64) u64 { - self.s +%= 0x9e3779b97f4a7c15; - - var z = self.s; - z = (z ^ (z >> 30)) *% 0xbf58476d1ce4e5b9; - z = (z ^ (z >> 27)) *% 0x94d049bb133111eb; - return z ^ (z >> 31); - } -}; - -test { - std.testing.refAllDecls(@This()); - _ = @import("rand/test.zig"); -} diff --git a/lib/std/rand/Ascon.zig b/lib/std/rand/Ascon.zig @@ -1,59 +0,0 @@ -//! CSPRNG based on the Reverie construction, a permutation-based PRNG -//! with forward security, instantiated with the Ascon(128,12,8) permutation. -//! -//! Compared to ChaCha, this PRNG has a much smaller state, and can be -//! a better choice for constrained environments. -//! -//! References: -//! - A Robust and Sponge-Like PRNG with Improved Efficiency https://eprint.iacr.org/2016/886.pdf -//! - Ascon https://ascon.iaik.tugraz.at/files/asconv12-nist.pdf - -const std = @import("std"); -const mem = std.mem; -const Random = std.rand.Random; -const Self = @This(); - -const Ascon = std.crypto.core.Ascon(.little); - -state: Ascon, - -const rate = 16; -pub const secret_seed_length = 32; - -/// The seed must be uniform, secret and `secret_seed_length` bytes long. -pub fn init(secret_seed: [secret_seed_length]u8) Self { - var self = Self{ .state = Ascon.initXof() }; - self.addEntropy(&secret_seed); - return self; -} - -/// Inserts entropy to refresh the internal state. -pub fn addEntropy(self: *Self, bytes: []const u8) void { - comptime std.debug.assert(secret_seed_length % rate == 0); - var i: usize = 0; - while (i + rate < bytes.len) : (i += rate) { - self.state.addBytes(bytes[i..][0..rate]); - self.state.permuteR(8); - } - if (i != bytes.len) self.state.addBytes(bytes[i..]); - self.state.permute(); -} - -/// Returns a `std.rand.Random` structure backed by the current RNG. -pub fn random(self: *Self) Random { - return Random.init(self, fill); -} - -/// Fills the buffer with random bytes. -pub fn fill(self: *Self, buf: []u8) void { - var i: usize = 0; - while (true) { - const left = buf.len - i; - const n = @min(left, rate); - self.state.extractBytes(buf[i..][0..n]); - if (left == 0) break; - self.state.permuteR(8); - i += n; - } - self.state.permuteRatchet(6, rate); -} diff --git a/lib/std/rand/ChaCha.zig b/lib/std/rand/ChaCha.zig @@ -1,98 +0,0 @@ -//! CSPRNG based on the ChaCha8 stream cipher, with forward security. -//! -//! References: -//! - Fast-key-erasure random-number generators https://blog.cr.yp.to/20170723-random.html - -const std = @import("std"); -const mem = std.mem; -const Random = std.rand.Random; -const Self = @This(); - -const Cipher = std.crypto.stream.chacha.ChaCha8IETF; - -const State = [8 * Cipher.block_length]u8; - -state: State, -offset: usize, - -const nonce = [_]u8{0} ** Cipher.nonce_length; - -pub const secret_seed_length = Cipher.key_length; - -/// The seed must be uniform, secret and `secret_seed_length` bytes long. -pub fn init(secret_seed: [secret_seed_length]u8) Self { - var self = Self{ .state = undefined, .offset = 0 }; - Cipher.stream(&self.state, 0, secret_seed, nonce); - return self; -} - -/// Inserts entropy to refresh the internal state. -pub fn addEntropy(self: *Self, bytes: []const u8) void { - var i: usize = 0; - while (i + Cipher.key_length <= bytes.len) : (i += Cipher.key_length) { - Cipher.xor( - self.state[0..Cipher.key_length], - self.state[0..Cipher.key_length], - 0, - bytes[i..][0..Cipher.key_length].*, - nonce, - ); - } - if (i < bytes.len) { - var k = [_]u8{0} ** Cipher.key_length; - const src = bytes[i..]; - @memcpy(k[0..src.len], src); - Cipher.xor( - self.state[0..Cipher.key_length], - self.state[0..Cipher.key_length], - 0, - k, - nonce, - ); - } - self.refill(); -} - -/// Returns a `std.rand.Random` structure backed by the current RNG. -pub fn random(self: *Self) Random { - return Random.init(self, fill); -} - -// Refills the buffer with random bytes, overwriting the previous key. -fn refill(self: *Self) void { - Cipher.stream(&self.state, 0, self.state[0..Cipher.key_length].*, nonce); - self.offset = 0; -} - -/// Fills the buffer with random bytes. -pub fn fill(self: *Self, buf_: []u8) void { - const bytes = self.state[Cipher.key_length..]; - var buf = buf_; - - const avail = bytes.len - self.offset; - if (avail > 0) { - // Bytes from the current block - const n = @min(avail, buf.len); - @memcpy(buf[0..n], bytes[self.offset..][0..n]); - @memset(bytes[self.offset..][0..n], 0); - buf = buf[n..]; - self.offset += n; - } - if (buf.len == 0) return; - - self.refill(); - - // Full blocks - while (buf.len >= bytes.len) { - @memcpy(buf[0..bytes.len], bytes); - buf = buf[bytes.len..]; - self.refill(); - } - - // Remaining bytes - if (buf.len > 0) { - @memcpy(buf, bytes[0..buf.len]); - @memset(bytes[0..buf.len], 0); - self.offset = buf.len; - } -} diff --git a/lib/std/rand/Isaac64.zig b/lib/std/rand/Isaac64.zig @@ -1,235 +0,0 @@ -//! ISAAC64 - http://www.burtleburtle.net/bob/rand/isaacafa.html -//! -//! Follows the general idea of the implementation from here with a few shortcuts. -//! https://doc.rust-lang.org/rand/src/rand/prng/isaac64.rs.html - -const std = @import("std"); -const Random = std.rand.Random; -const mem = std.mem; -const Isaac64 = @This(); - -r: [256]u64, -m: [256]u64, -a: u64, -b: u64, -c: u64, -i: usize, - -pub fn init(init_s: u64) Isaac64 { - var isaac = Isaac64{ - .r = undefined, - .m = undefined, - .a = undefined, - .b = undefined, - .c = undefined, - .i = undefined, - }; - - // seed == 0 => same result as the unseeded reference implementation - isaac.seed(init_s, 1); - return isaac; -} - -pub fn random(self: *Isaac64) Random { - return Random.init(self, fill); -} - -fn step(self: *Isaac64, mix: u64, base: usize, comptime m1: usize, comptime m2: usize) void { - const x = self.m[base + m1]; - self.a = mix +% self.m[base + m2]; - - const y = self.a +% self.b +% self.m[@as(usize, @intCast((x >> 3) % self.m.len))]; - self.m[base + m1] = y; - - self.b = x +% self.m[@as(usize, @intCast((y >> 11) % self.m.len))]; - self.r[self.r.len - 1 - base - m1] = self.b; -} - -fn refill(self: *Isaac64) void { - const midpoint = self.r.len / 2; - - self.c +%= 1; - self.b +%= self.c; - - { - var i: usize = 0; - while (i < midpoint) : (i += 4) { - self.step(~(self.a ^ (self.a << 21)), i + 0, 0, midpoint); - self.step(self.a ^ (self.a >> 5), i + 1, 0, midpoint); - self.step(self.a ^ (self.a << 12), i + 2, 0, midpoint); - self.step(self.a ^ (self.a >> 33), i + 3, 0, midpoint); - } - } - - { - var i: usize = 0; - while (i < midpoint) : (i += 4) { - self.step(~(self.a ^ (self.a << 21)), i + 0, midpoint, 0); - self.step(self.a ^ (self.a >> 5), i + 1, midpoint, 0); - self.step(self.a ^ (self.a << 12), i + 2, midpoint, 0); - self.step(self.a ^ (self.a >> 33), i + 3, midpoint, 0); - } - } - - self.i = 0; -} - -fn next(self: *Isaac64) u64 { - if (self.i >= self.r.len) { - self.refill(); - } - - const value = self.r[self.i]; - self.i += 1; - return value; -} - -fn seed(self: *Isaac64, init_s: u64, comptime rounds: usize) void { - // We ignore the multi-pass requirement since we don't currently expose full access to - // seeding the self.m array completely. - @memset(self.m[0..], 0); - self.m[0] = init_s; - - // prescrambled golden ratio constants - var a = [_]u64{ - 0x647c4677a2884b7c, - 0xb9f8b322c73ac862, - 0x8c0ea5053d4712a0, - 0xb29b2e824a595524, - 0x82f053db8355e0ce, - 0x48fe4a0fa5a09315, - 0xae985bf2cbfc89ed, - 0x98f5704f6c44c0ab, - }; - - comptime var i: usize = 0; - inline while (i < rounds) : (i += 1) { - var j: usize = 0; - while (j < self.m.len) : (j += 8) { - comptime var x1: usize = 0; - inline while (x1 < 8) : (x1 += 1) { - a[x1] +%= self.m[j + x1]; - } - - a[0] -%= a[4]; - a[5] ^= a[7] >> 9; - a[7] +%= a[0]; - a[1] -%= a[5]; - a[6] ^= a[0] << 9; - a[0] +%= a[1]; - a[2] -%= a[6]; - a[7] ^= a[1] >> 23; - a[1] +%= a[2]; - a[3] -%= a[7]; - a[0] ^= a[2] << 15; - a[2] +%= a[3]; - a[4] -%= a[0]; - a[1] ^= a[3] >> 14; - a[3] +%= a[4]; - a[5] -%= a[1]; - a[2] ^= a[4] << 20; - a[4] +%= a[5]; - a[6] -%= a[2]; - a[3] ^= a[5] >> 17; - a[5] +%= a[6]; - a[7] -%= a[3]; - a[4] ^= a[6] << 14; - a[6] +%= a[7]; - - comptime var x2: usize = 0; - inline while (x2 < 8) : (x2 += 1) { - self.m[j + x2] = a[x2]; - } - } - } - - @memset(self.r[0..], 0); - self.a = 0; - self.b = 0; - self.c = 0; - self.i = self.r.len; // trigger refill on first value -} - -pub fn fill(self: *Isaac64, buf: []u8) void { - var i: usize = 0; - const aligned_len = buf.len - (buf.len & 7); - - // Fill complete 64-byte segments - while (i < aligned_len) : (i += 8) { - var n = self.next(); - comptime var j: usize = 0; - inline while (j < 8) : (j += 1) { - buf[i + j] = @as(u8, @truncate(n)); - n >>= 8; - } - } - - // Fill trailing, ignoring excess (cut the stream). - if (i != buf.len) { - var n = self.next(); - while (i < buf.len) : (i += 1) { - buf[i] = @as(u8, @truncate(n)); - n >>= 8; - } - } -} - -test "isaac64 sequence" { - var r = Isaac64.init(0); - - // from reference implementation - const seq = [_]u64{ - 0xf67dfba498e4937c, - 0x84a5066a9204f380, - 0xfee34bd5f5514dbb, - 0x4d1664739b8f80d6, - 0x8607459ab52a14aa, - 0x0e78bc5a98529e49, - 0xfe5332822ad13777, - 0x556c27525e33d01a, - 0x08643ca615f3149f, - 0xd0771faf3cb04714, - 0x30e86f68a37b008d, - 0x3074ebc0488a3adf, - 0x270645ea7a2790bc, - 0x5601a0a8d3763c6a, - 0x2f83071f53f325dd, - 0xb9090f3d42d2d2ea, - }; - - for (seq) |s| { - try std.testing.expect(s == r.next()); - } -} - -test "isaac64 fill" { - var r = Isaac64.init(0); - - // from reference implementation - const seq = [_]u64{ - 0xf67dfba498e4937c, - 0x84a5066a9204f380, - 0xfee34bd5f5514dbb, - 0x4d1664739b8f80d6, - 0x8607459ab52a14aa, - 0x0e78bc5a98529e49, - 0xfe5332822ad13777, - 0x556c27525e33d01a, - 0x08643ca615f3149f, - 0xd0771faf3cb04714, - 0x30e86f68a37b008d, - 0x3074ebc0488a3adf, - 0x270645ea7a2790bc, - 0x5601a0a8d3763c6a, - 0x2f83071f53f325dd, - 0xb9090f3d42d2d2ea, - }; - - for (seq) |s| { - var buf0: [8]u8 = undefined; - var buf1: [7]u8 = undefined; - std.mem.writeInt(u64, &buf0, s, .little); - r.fill(&buf1); - try std.testing.expect(std.mem.eql(u8, buf0[0..7], buf1[0..])); - } -} diff --git a/lib/std/rand/Pcg.zig b/lib/std/rand/Pcg.zig @@ -1,122 +0,0 @@ -//! PCG32 - http://www.pcg-random.org/ -//! -//! PRNG - -const std = @import("std"); -const Random = std.rand.Random; -const Pcg = @This(); - -const default_multiplier = 6364136223846793005; - -s: u64, -i: u64, - -pub fn init(init_s: u64) Pcg { - var pcg = Pcg{ - .s = undefined, - .i = undefined, - }; - - pcg.seed(init_s); - return pcg; -} - -pub fn random(self: *Pcg) Random { - return Random.init(self, fill); -} - -fn next(self: *Pcg) u32 { - const l = self.s; - self.s = l *% default_multiplier +% (self.i | 1); - - const xor_s: u32 = @truncate(((l >> 18) ^ l) >> 27); - const rot: u32 = @intCast(l >> 59); - - return (xor_s >> @as(u5, @intCast(rot))) | (xor_s << @as(u5, @intCast((0 -% rot) & 31))); -} - -fn seed(self: *Pcg, init_s: u64) void { - // Pcg requires 128-bits of seed. - var gen = std.rand.SplitMix64.init(init_s); - self.seedTwo(gen.next(), gen.next()); -} - -fn seedTwo(self: *Pcg, init_s: u64, init_i: u64) void { - self.s = 0; - self.i = (init_s << 1) | 1; - self.s = self.s *% default_multiplier +% self.i; - self.s +%= init_i; - self.s = self.s *% default_multiplier +% self.i; -} - -pub fn fill(self: *Pcg, buf: []u8) void { - var i: usize = 0; - const aligned_len = buf.len - (buf.len & 3); - - // Complete 4 byte segments. - while (i < aligned_len) : (i += 4) { - var n = self.next(); - comptime var j: usize = 0; - inline while (j < 4) : (j += 1) { - buf[i + j] = @as(u8, @truncate(n)); - n >>= 8; - } - } - - // Remaining. (cuts the stream) - if (i != buf.len) { - var n = self.next(); - while (i < buf.len) : (i += 1) { - buf[i] = @as(u8, @truncate(n)); - n >>= 8; - } - } -} - -test "pcg sequence" { - var r = Pcg.init(0); - const s0: u64 = 0x9394bf54ce5d79de; - const s1: u64 = 0x84e9c579ef59bbf7; - r.seedTwo(s0, s1); - - const seq = [_]u32{ - 2881561918, - 3063928540, - 1199791034, - 2487695858, - 1479648952, - 3247963454, - }; - - for (seq) |s| { - try std.testing.expect(s == r.next()); - } -} - -test "pcg fill" { - var r = Pcg.init(0); - const s0: u64 = 0x9394bf54ce5d79de; - const s1: u64 = 0x84e9c579ef59bbf7; - r.seedTwo(s0, s1); - - const seq = [_]u32{ - 2881561918, - 3063928540, - 1199791034, - 2487695858, - 1479648952, - 3247963454, - }; - - var i: u32 = 0; - while (i < seq.len) : (i += 2) { - var buf0: [8]u8 = undefined; - std.mem.writeInt(u32, buf0[0..4], seq[i], .little); - std.mem.writeInt(u32, buf0[4..8], seq[i + 1], .little); - - var buf1: [7]u8 = undefined; - r.fill(&buf1); - - try std.testing.expect(std.mem.eql(u8, buf0[0..7], buf1[0..])); - } -} diff --git a/lib/std/rand/RomuTrio.zig b/lib/std/rand/RomuTrio.zig @@ -1,132 +0,0 @@ -// Website: romu-random.org -// Reference paper: http://arxiv.org/abs/2002.11331 -// Beware: this PRNG is trivially predictable. While fast, it should *never* be used for cryptographic purposes. - -const std = @import("std"); -const Random = std.rand.Random; -const math = std.math; -const RomuTrio = @This(); - -x_state: u64, -y_state: u64, -z_state: u64, // set to nonzero seed - -pub fn init(init_s: u64) RomuTrio { - var x = RomuTrio{ .x_state = undefined, .y_state = undefined, .z_state = undefined }; - x.seed(init_s); - return x; -} - -pub fn random(self: *RomuTrio) Random { - return Random.init(self, fill); -} - -fn next(self: *RomuTrio) u64 { - const xp = self.x_state; - const yp = self.y_state; - const zp = self.z_state; - self.x_state = 15241094284759029579 *% zp; - self.y_state = yp -% xp; - self.y_state = std.math.rotl(u64, self.y_state, 12); - self.z_state = zp -% yp; - self.z_state = std.math.rotl(u64, self.z_state, 44); - return xp; -} - -pub fn seedWithBuf(self: *RomuTrio, buf: [24]u8) void { - const seed_buf = @as([3]u64, @bitCast(buf)); - self.x_state = seed_buf[0]; - self.y_state = seed_buf[1]; - self.z_state = seed_buf[2]; -} - -pub fn seed(self: *RomuTrio, init_s: u64) void { - // RomuTrio requires 192-bits of seed. - var gen = std.rand.SplitMix64.init(init_s); - - self.x_state = gen.next(); - self.y_state = gen.next(); - self.z_state = gen.next(); -} - -pub fn fill(self: *RomuTrio, buf: []u8) void { - var i: usize = 0; - const aligned_len = buf.len - (buf.len & 7); - - // Complete 8 byte segments. - while (i < aligned_len) : (i += 8) { - var n = self.next(); - comptime var j: usize = 0; - inline while (j < 8) : (j += 1) { - buf[i + j] = @as(u8, @truncate(n)); - n >>= 8; - } - } - - // Remaining. (cuts the stream) - if (i != buf.len) { - var n = self.next(); - while (i < buf.len) : (i += 1) { - buf[i] = @as(u8, @truncate(n)); - n >>= 8; - } - } -} - -test "RomuTrio sequence" { - // Unfortunately there does not seem to be an official test sequence. - var r = RomuTrio.init(0); - - const seq = [_]u64{ - 16294208416658607535, - 13964609475759908645, - 4703697494102998476, - 3425221541186733346, - 2285772463536419399, - 9454187757529463048, - 13695907680080547496, - 8328236714879408626, - 12323357569716880909, - 12375466223337721820, - }; - - for (seq) |s| { - try std.testing.expectEqual(s, r.next()); - } -} - -test "RomuTrio fill" { - // Unfortunately there does not seem to be an official test sequence. - var r = RomuTrio.init(0); - - const seq = [_]u64{ - 16294208416658607535, - 13964609475759908645, - 4703697494102998476, - 3425221541186733346, - 2285772463536419399, - 9454187757529463048, - 13695907680080547496, - 8328236714879408626, - 12323357569716880909, - 12375466223337721820, - }; - - for (seq) |s| { - var buf0: [8]u8 = undefined; - var buf1: [7]u8 = undefined; - std.mem.writeInt(u64, &buf0, s, .little); - r.fill(&buf1); - try std.testing.expect(std.mem.eql(u8, buf0[0..7], buf1[0..])); - } -} - -test "RomuTrio buf seeding test" { - const buf0 = @as([24]u8, @bitCast([3]u64{ 16294208416658607535, 13964609475759908645, 4703697494102998476 })); - const resulting_state = .{ .x = 16294208416658607535, .y = 13964609475759908645, .z = 4703697494102998476 }; - var r = RomuTrio.init(0); - r.seedWithBuf(buf0); - try std.testing.expect(r.x_state == resulting_state.x); - try std.testing.expect(r.y_state == resulting_state.y); - try std.testing.expect(r.z_state == resulting_state.z); -} diff --git a/lib/std/rand/Sfc64.zig b/lib/std/rand/Sfc64.zig @@ -1,132 +0,0 @@ -//! Sfc64 pseudo-random number generator from Practically Random. -//! Fastest engine of pracrand and smallest footprint. -//! See http://pracrand.sourceforge.net/ - -const std = @import("std"); -const Random = std.rand.Random; -const math = std.math; -const Sfc64 = @This(); - -a: u64 = undefined, -b: u64 = undefined, -c: u64 = undefined, -counter: u64 = undefined, - -const Rotation = 24; -const RightShift = 11; -const LeftShift = 3; - -pub fn init(init_s: u64) Sfc64 { - var x = Sfc64{}; - - x.seed(init_s); - return x; -} - -pub fn random(self: *Sfc64) Random { - return Random.init(self, fill); -} - -fn next(self: *Sfc64) u64 { - const tmp = self.a +% self.b +% self.counter; - self.counter += 1; - self.a = self.b ^ (self.b >> RightShift); - self.b = self.c +% (self.c << LeftShift); - self.c = math.rotl(u64, self.c, Rotation) +% tmp; - return tmp; -} - -fn seed(self: *Sfc64, init_s: u64) void { - self.a = init_s; - self.b = init_s; - self.c = init_s; - self.counter = 1; - var i: u32 = 0; - while (i < 12) : (i += 1) { - _ = self.next(); - } -} - -pub fn fill(self: *Sfc64, buf: []u8) void { - var i: usize = 0; - const aligned_len = buf.len - (buf.len & 7); - - // Complete 8 byte segments. - while (i < aligned_len) : (i += 8) { - var n = self.next(); - comptime var j: usize = 0; - inline while (j < 8) : (j += 1) { - buf[i + j] = @as(u8, @truncate(n)); - n >>= 8; - } - } - - // Remaining. (cuts the stream) - if (i != buf.len) { - var n = self.next(); - while (i < buf.len) : (i += 1) { - buf[i] = @as(u8, @truncate(n)); - n >>= 8; - } - } -} - -test "Sfc64 sequence" { - // Unfortunately there does not seem to be an official test sequence. - var r = Sfc64.init(0); - - const seq = [_]u64{ - 0x3acfa029e3cc6041, - 0xf5b6515bf2ee419c, - 0x1259635894a29b61, - 0xb6ae75395f8ebd6, - 0x225622285ce302e2, - 0x520d28611395cb21, - 0xdb909c818901599d, - 0x8ffd195365216f57, - 0xe8c4ad5e258ac04a, - 0x8f8ef2c89fdb63ca, - 0xf9865b01d98d8e2f, - 0x46555871a65d08ba, - 0x66868677c6298fcd, - 0x2ce15a7e6329f57d, - 0xb2f1833ca91ca79, - 0x4b0890ac9bf453ca, - }; - - for (seq) |s| { - try std.testing.expectEqual(s, r.next()); - } -} - -test "Sfc64 fill" { - // Unfortunately there does not seem to be an official test sequence. - var r = Sfc64.init(0); - - const seq = [_]u64{ - 0x3acfa029e3cc6041, - 0xf5b6515bf2ee419c, - 0x1259635894a29b61, - 0xb6ae75395f8ebd6, - 0x225622285ce302e2, - 0x520d28611395cb21, - 0xdb909c818901599d, - 0x8ffd195365216f57, - 0xe8c4ad5e258ac04a, - 0x8f8ef2c89fdb63ca, - 0xf9865b01d98d8e2f, - 0x46555871a65d08ba, - 0x66868677c6298fcd, - 0x2ce15a7e6329f57d, - 0xb2f1833ca91ca79, - 0x4b0890ac9bf453ca, - }; - - for (seq) |s| { - var buf0: [8]u8 = undefined; - var buf1: [7]u8 = undefined; - std.mem.writeInt(u64, &buf0, s, .little); - r.fill(&buf1); - try std.testing.expect(std.mem.eql(u8, buf0[0..7], buf1[0..])); - } -} diff --git a/lib/std/rand/Xoroshiro128.zig b/lib/std/rand/Xoroshiro128.zig @@ -1,147 +0,0 @@ -//! Xoroshiro128+ - http://xoroshiro.di.unimi.it/ -//! -//! PRNG - -const std = @import("std"); -const Random = std.rand.Random; -const math = std.math; -const Xoroshiro128 = @This(); - -s: [2]u64, - -pub fn init(init_s: u64) Xoroshiro128 { - var x = Xoroshiro128{ .s = undefined }; - - x.seed(init_s); - return x; -} - -pub fn random(self: *Xoroshiro128) Random { - return Random.init(self, fill); -} - -pub fn next(self: *Xoroshiro128) u64 { - const s0 = self.s[0]; - var s1 = self.s[1]; - const r = s0 +% s1; - - s1 ^= s0; - self.s[0] = math.rotl(u64, s0, @as(u8, 55)) ^ s1 ^ (s1 << 14); - self.s[1] = math.rotl(u64, s1, @as(u8, 36)); - - return r; -} - -// Skip 2^64 places ahead in the sequence -pub fn jump(self: *Xoroshiro128) void { - var s0: u64 = 0; - var s1: u64 = 0; - - const table = [_]u64{ - 0xbeac0467eba5facb, - 0xd86b048b86aa9922, - }; - - inline for (table) |entry| { - var b: usize = 0; - while (b < 64) : (b += 1) { - if ((entry & (@as(u64, 1) << @as(u6, @intCast(b)))) != 0) { - s0 ^= self.s[0]; - s1 ^= self.s[1]; - } - _ = self.next(); - } - } - - self.s[0] = s0; - self.s[1] = s1; -} - -pub fn seed(self: *Xoroshiro128, init_s: u64) void { - // Xoroshiro requires 128-bits of seed. - var gen = std.rand.SplitMix64.init(init_s); - - self.s[0] = gen.next(); - self.s[1] = gen.next(); -} - -pub fn fill(self: *Xoroshiro128, buf: []u8) void { - var i: usize = 0; - const aligned_len = buf.len - (buf.len & 7); - - // Complete 8 byte segments. - while (i < aligned_len) : (i += 8) { - var n = self.next(); - comptime var j: usize = 0; - inline while (j < 8) : (j += 1) { - buf[i + j] = @as(u8, @truncate(n)); - n >>= 8; - } - } - - // Remaining. (cuts the stream) - if (i != buf.len) { - var n = self.next(); - while (i < buf.len) : (i += 1) { - buf[i] = @as(u8, @truncate(n)); - n >>= 8; - } - } -} - -test "xoroshiro sequence" { - var r = Xoroshiro128.init(0); - r.s[0] = 0xaeecf86f7878dd75; - r.s[1] = 0x01cd153642e72622; - - const seq1 = [_]u64{ - 0xb0ba0da5bb600397, - 0x18a08afde614dccc, - 0xa2635b956a31b929, - 0xabe633c971efa045, - 0x9ac19f9706ca3cac, - 0xf62b426578c1e3fb, - }; - - for (seq1) |s| { - try std.testing.expect(s == r.next()); - } - - r.jump(); - - const seq2 = [_]u64{ - 0x95344a13556d3e22, - 0xb4fb32dafa4d00df, - 0xb2011d9ccdcfe2dd, - 0x05679a9b2119b908, - 0xa860a1da7c9cd8a0, - 0x658a96efe3f86550, - }; - - for (seq2) |s| { - try std.testing.expect(s == r.next()); - } -} - -test "xoroshiro fill" { - var r = Xoroshiro128.init(0); - r.s[0] = 0xaeecf86f7878dd75; - r.s[1] = 0x01cd153642e72622; - - const seq = [_]u64{ - 0xb0ba0da5bb600397, - 0x18a08afde614dccc, - 0xa2635b956a31b929, - 0xabe633c971efa045, - 0x9ac19f9706ca3cac, - 0xf62b426578c1e3fb, - }; - - for (seq) |s| { - var buf0: [8]u8 = undefined; - var buf1: [7]u8 = undefined; - std.mem.writeInt(u64, &buf0, s, .little); - r.fill(&buf1); - try std.testing.expect(std.mem.eql(u8, buf0[0..7], buf1[0..])); - } -} diff --git a/lib/std/rand/Xoshiro256.zig b/lib/std/rand/Xoshiro256.zig @@ -1,146 +0,0 @@ -//! Xoshiro256++ - http://xoroshiro.di.unimi.it/ -//! -//! PRNG - -const std = @import("std"); -const Random = std.rand.Random; -const math = std.math; -const Xoshiro256 = @This(); - -s: [4]u64, - -pub fn init(init_s: u64) Xoshiro256 { - var x = Xoshiro256{ - .s = undefined, - }; - - x.seed(init_s); - return x; -} - -pub fn random(self: *Xoshiro256) Random { - return Random.init(self, fill); -} - -pub fn next(self: *Xoshiro256) u64 { - const r = math.rotl(u64, self.s[0] +% self.s[3], 23) +% self.s[0]; - - const t = self.s[1] << 17; - - self.s[2] ^= self.s[0]; - self.s[3] ^= self.s[1]; - self.s[1] ^= self.s[2]; - self.s[0] ^= self.s[3]; - - self.s[2] ^= t; - - self.s[3] = math.rotl(u64, self.s[3], 45); - - return r; -} - -// Skip 2^128 places ahead in the sequence -pub fn jump(self: *Xoshiro256) void { - var s: u256 = 0; - - var table: u256 = 0x39abdc4529b1661ca9582618e03fc9aad5a61266f0c9392c180ec6d33cfd0aba; - - while (table != 0) : (table >>= 1) { - if (@as(u1, @truncate(table)) != 0) { - s ^= @as(u256, @bitCast(self.s)); - } - _ = self.next(); - } - - self.s = @as([4]u64, @bitCast(s)); -} - -pub fn seed(self: *Xoshiro256, init_s: u64) void { - // Xoshiro requires 256-bits of seed. - var gen = std.rand.SplitMix64.init(init_s); - - self.s[0] = gen.next(); - self.s[1] = gen.next(); - self.s[2] = gen.next(); - self.s[3] = gen.next(); -} - -pub fn fill(self: *Xoshiro256, buf: []u8) void { - var i: usize = 0; - const aligned_len = buf.len - (buf.len & 7); - - // Complete 8 byte segments. - while (i < aligned_len) : (i += 8) { - var n = self.next(); - comptime var j: usize = 0; - inline while (j < 8) : (j += 1) { - buf[i + j] = @as(u8, @truncate(n)); - n >>= 8; - } - } - - // Remaining. (cuts the stream) - if (i != buf.len) { - var n = self.next(); - while (i < buf.len) : (i += 1) { - buf[i] = @as(u8, @truncate(n)); - n >>= 8; - } - } -} - -test "xoroshiro sequence" { - if (@import("builtin").zig_backend == .stage2_c) return error.SkipZigTest; - if (@import("builtin").zig_backend == .stage2_x86_64) return error.SkipZigTest; - - var r = Xoshiro256.init(0); - - const seq1 = [_]u64{ - 0x53175d61490b23df, - 0x61da6f3dc380d507, - 0x5c0fdf91ec9a7bfc, - 0x02eebf8c3bbe5e1a, - 0x7eca04ebaf4a5eea, - 0x0543c37757f08d9a, - }; - - for (seq1) |s| { - try std.testing.expect(s == r.next()); - } - - r.jump(); - - const seq2 = [_]u64{ - 0xae1db5c5e27807be, - 0xb584c6a7fd8709fe, - 0xc46a0ee9330fb6e, - 0xdc0c9606f49ed76e, - 0x1f5bb6540f6651fb, - 0x72fa2ca734601488, - }; - - for (seq2) |s| { - try std.testing.expect(s == r.next()); - } -} - -test "xoroshiro fill" { - var r = Xoshiro256.init(0); - - const seq = [_]u64{ - 0x53175d61490b23df, - 0x61da6f3dc380d507, - 0x5c0fdf91ec9a7bfc, - 0x02eebf8c3bbe5e1a, - 0x7eca04ebaf4a5eea, - 0x0543c37757f08d9a, - }; - - for (seq) |s| { - var buf0: [8]u8 = undefined; - var buf1: [7]u8 = undefined; - std.mem.writeInt(u64, &buf0, s, .little); - r.fill(&buf1); - try std.testing.expect(std.mem.eql(u8, buf0[0..7], buf1[0..])); - } -} diff --git a/lib/std/rand/benchmark.zig b/lib/std/rand/benchmark.zig @@ -1,217 +0,0 @@ -// zig run -O ReleaseFast --zig-lib-dir ../.. benchmark.zig - -const std = @import("std"); -const builtin = @import("builtin"); -const time = std.time; -const Timer = time.Timer; -const rand = std.rand; - -const KiB = 1024; -const MiB = 1024 * KiB; -const GiB = 1024 * MiB; - -const Rng = struct { - ty: type, - name: []const u8, - init_u8s: ?[]const u8 = null, - init_u64: ?u64 = null, -}; - -const prngs = [_]Rng{ - Rng{ - .ty = rand.Isaac64, - .name = "isaac64", - .init_u64 = 0, - }, - Rng{ - .ty = rand.Pcg, - .name = "pcg", - .init_u64 = 0, - }, - Rng{ - .ty = rand.RomuTrio, - .name = "romutrio", - .init_u64 = 0, - }, - Rng{ - .ty = std.rand.Sfc64, - .name = "sfc64", - .init_u64 = 0, - }, - Rng{ - .ty = std.rand.Xoroshiro128, - .name = "xoroshiro128", - .init_u64 = 0, - }, - Rng{ - .ty = std.rand.Xoshiro256, - .name = "xoshiro256", - .init_u64 = 0, - }, -}; - -const csprngs = [_]Rng{ - Rng{ - .ty = rand.Ascon, - .name = "ascon", - .init_u8s = &[_]u8{0} ** 32, - }, - Rng{ - .ty = rand.ChaCha, - .name = "chacha", - .init_u8s = &[_]u8{0} ** 32, - }, -}; - -const Result = struct { - throughput: u64, -}; - -const long_block_size: usize = 8 * 8192; -const short_block_size: usize = 8; - -pub fn benchmark(comptime H: anytype, bytes: usize, comptime block_size: usize) !Result { - var rng = blk: { - if (H.init_u8s) |init| { - break :blk H.ty.init(init[0..].*); - } - if (H.init_u64) |init| { - break :blk H.ty.init(init); - } - break :blk H.ty.init(); - }; - - var block: [block_size]u8 = undefined; - - var offset: usize = 0; - var timer = try Timer.start(); - const start = timer.lap(); - while (offset < bytes) : (offset += block.len) { - rng.fill(block[0..]); - } - const end = timer.read(); - - const elapsed_s = @as(f64, @floatFromInt(end - start)) / time.ns_per_s; - const throughput = @as(u64, @intFromFloat(@as(f64, @floatFromInt(bytes)) / elapsed_s)); - - std.debug.assert(rng.random().int(u64) != 0); - - return Result{ - .throughput = throughput, - }; -} - -fn usage() void { - std.debug.print( - \\throughput_test [options] - \\ - \\Options: - \\ --filter [test-name] - \\ --count [int] - \\ --prngs-only - \\ --csprngs-only - \\ --short-only - \\ --long-only - \\ --help - \\ - , .{}); -} - -fn mode(comptime x: comptime_int) comptime_int { - return if (builtin.mode == .Debug) x / 64 else x; -} - -pub fn main() !void { - const stdout = std.io.getStdOut().writer(); - - var buffer: [1024]u8 = undefined; - var fixed = std.heap.FixedBufferAllocator.init(buffer[0..]); - const args = try std.process.argsAlloc(fixed.allocator()); - - var filter: ?[]u8 = ""; - var count: usize = mode(128 * MiB); - var bench_prngs = true; - var bench_csprngs = true; - var bench_long = true; - var bench_short = true; - - var i: usize = 1; - while (i < args.len) : (i += 1) { - if (std.mem.eql(u8, args[i], "--mode")) { - try stdout.print("{}\n", .{builtin.mode}); - return; - } else if (std.mem.eql(u8, args[i], "--filter")) { - i += 1; - if (i == args.len) { - usage(); - std.os.exit(1); - } - - filter = args[i]; - } else if (std.mem.eql(u8, args[i], "--count")) { - i += 1; - if (i == args.len) { - usage(); - std.os.exit(1); - } - - const c = try std.fmt.parseUnsigned(usize, args[i], 10); - count = c * MiB; - } else if (std.mem.eql(u8, args[i], "--csprngs-only")) { - bench_prngs = false; - } else if (std.mem.eql(u8, args[i], "--prngs-only")) { - bench_csprngs = false; - } else if (std.mem.eql(u8, args[i], "--short-only")) { - bench_long = false; - } else if (std.mem.eql(u8, args[i], "--long-only")) { - bench_short = false; - } else if (std.mem.eql(u8, args[i], "--help")) { - usage(); - return; - } else { - usage(); - std.os.exit(1); - } - } - - if (bench_prngs) { - if (bench_long) { - inline for (prngs) |R| { - if (filter == null or std.mem.indexOf(u8, R.name, filter.?) != null) { - try stdout.print("{s} (long outputs)\n", .{R.name}); - const result_long = try benchmark(R, count, long_block_size); - try stdout.print(" {:5} MiB/s\n", .{result_long.throughput / (1 * MiB)}); - } - } - } - if (bench_short) { - inline for (prngs) |R| { - if (filter == null or std.mem.indexOf(u8, R.name, filter.?) != null) { - try stdout.print("{s} (short outputs)\n", .{R.name}); - const result_short = try benchmark(R, count, short_block_size); - try stdout.print(" {:5} MiB/s\n", .{result_short.throughput / (1 * MiB)}); - } - } - } - } - if (bench_csprngs) { - if (bench_long) { - inline for (csprngs) |R| { - if (filter == null or std.mem.indexOf(u8, R.name, filter.?) != null) { - try stdout.print("{s} (cryptographic, long outputs)\n", .{R.name}); - const result_long = try benchmark(R, count, long_block_size); - try stdout.print(" {:5} MiB/s\n", .{result_long.throughput / (1 * MiB)}); - } - } - } - if (bench_short) { - inline for (csprngs) |R| { - if (filter == null or std.mem.indexOf(u8, R.name, filter.?) != null) { - try stdout.print("{s} (cryptographic, short outputs)\n", .{R.name}); - const result_short = try benchmark(R, count, short_block_size); - try stdout.print(" {:5} MiB/s\n", .{result_short.throughput / (1 * MiB)}); - } - } - } - } -} diff --git a/lib/std/rand/test.zig b/lib/std/rand/test.zig @@ -1,473 +0,0 @@ -const std = @import("../std.zig"); -const math = std.math; -const DefaultPrng = std.rand.DefaultPrng; -const Random = std.rand.Random; -const SplitMix64 = std.rand.SplitMix64; -const DefaultCsprng = std.rand.DefaultCsprng; -const expect = std.testing.expect; -const expectEqual = std.testing.expectEqual; - -const SequentialPrng = struct { - const Self = @This(); - next_value: u8, - - pub fn init() Self { - return Self{ - .next_value = 0, - }; - } - - pub fn random(self: *Self) Random { - return Random.init(self, fill); - } - - pub fn fill(self: *Self, buf: []u8) void { - for (buf) |*b| { - b.* = self.next_value; - } - self.next_value +%= 1; - } -}; - -/// Do not use this PRNG! It is meant to be predictable, for the purposes of test reproducibility and coverage. -/// Its output is just a repeat of a user-specified byte pattern. -/// Name is a reference to this comic: https://dilbert.com/strip/2001-10-25 -const Dilbert = struct { - pattern: []const u8 = undefined, - curr_idx: usize = 0, - - pub fn init(pattern: []const u8) !Dilbert { - if (pattern.len == 0) - return error.EmptyPattern; - var self = Dilbert{}; - self.pattern = pattern; - self.curr_idx = 0; - return self; - } - - pub fn random(self: *Dilbert) Random { - return Random.init(self, fill); - } - - pub fn fill(self: *Dilbert, buf: []u8) void { - for (buf) |*byte| { - byte.* = self.pattern[self.curr_idx]; - self.curr_idx = (self.curr_idx + 1) % self.pattern.len; - } - } - - test "Dilbert fill" { - var r = try Dilbert.init("9nine"); - - const seq = [_]u64{ - 0x396E696E65396E69, - 0x6E65396E696E6539, - 0x6E696E65396E696E, - 0x65396E696E65396E, - 0x696E65396E696E65, - }; - - for (seq) |s| { - var buf0: [8]u8 = undefined; - var buf1: [8]u8 = undefined; - std.mem.writeInt(u64, &buf0, s, .big); - r.fill(&buf1); - try std.testing.expect(std.mem.eql(u8, buf0[0..], buf1[0..])); - } - } -}; - -test "Random int" { - try testRandomInt(); - try comptime testRandomInt(); -} -fn testRandomInt() !void { - var rng = SequentialPrng.init(); - const random = rng.random(); - - try expect(random.int(u0) == 0); - - rng.next_value = 0; - try expect(random.int(u1) == 0); - try expect(random.int(u1) == 1); - try expect(random.int(u2) == 2); - try expect(random.int(u2) == 3); - try expect(random.int(u2) == 0); - - rng.next_value = 0xff; - try expect(random.int(u8) == 0xff); - rng.next_value = 0x11; - try expect(random.int(u8) == 0x11); - - rng.next_value = 0xff; - try expect(random.int(u32) == 0xffffffff); - rng.next_value = 0x11; - try expect(random.int(u32) == 0x11111111); - - rng.next_value = 0xff; - try expect(random.int(i32) == -1); - rng.next_value = 0x11; - try expect(random.int(i32) == 0x11111111); - - rng.next_value = 0xff; - try expect(random.int(i8) == -1); - rng.next_value = 0x11; - try expect(random.int(i8) == 0x11); - - rng.next_value = 0xff; - try expect(random.int(u33) == 0x1ffffffff); - rng.next_value = 0xff; - try expect(random.int(i1) == -1); - rng.next_value = 0xff; - try expect(random.int(i2) == -1); - rng.next_value = 0xff; - try expect(random.int(i33) == -1); -} - -test "Random boolean" { - try testRandomBoolean(); - try comptime testRandomBoolean(); -} -fn testRandomBoolean() !void { - var rng = SequentialPrng.init(); - const random = rng.random(); - - try expect(random.boolean() == false); - try expect(random.boolean() == true); - try expect(random.boolean() == false); - try expect(random.boolean() == true); -} - -test "Random enum" { - try testRandomEnumValue(); - try comptime testRandomEnumValue(); -} -fn testRandomEnumValue() !void { - const TestEnum = enum { - First, - Second, - Third, - }; - var rng = SequentialPrng.init(); - const random = rng.random(); - rng.next_value = 0; - try expect(random.enumValue(TestEnum) == TestEnum.First); - try expect(random.enumValue(TestEnum) == TestEnum.First); - try expect(random.enumValue(TestEnum) == TestEnum.First); -} - -test "Random intLessThan" { - @setEvalBranchQuota(10000); - try testRandomIntLessThan(); - try comptime testRandomIntLessThan(); -} -fn testRandomIntLessThan() !void { - var rng = SequentialPrng.init(); - const random = rng.random(); - - rng.next_value = 0xff; - try expect(random.uintLessThan(u8, 4) == 3); - try expect(rng.next_value == 0); - try expect(random.uintLessThan(u8, 4) == 0); - try expect(rng.next_value == 1); - - rng.next_value = 0; - try expect(random.uintLessThan(u64, 32) == 0); - - // trigger the bias rejection code path - rng.next_value = 0; - try expect(random.uintLessThan(u8, 3) == 0); - // verify we incremented twice - try expect(rng.next_value == 2); - - rng.next_value = 0xff; - try expect(random.intRangeLessThan(u8, 0, 0x80) == 0x7f); - rng.next_value = 0xff; - try expect(random.intRangeLessThan(u8, 0x7f, 0xff) == 0xfe); - - rng.next_value = 0xff; - try expect(random.intRangeLessThan(i8, 0, 0x40) == 0x3f); - rng.next_value = 0xff; - try expect(random.intRangeLessThan(i8, -0x40, 0x40) == 0x3f); - rng.next_value = 0xff; - try expect(random.intRangeLessThan(i8, -0x80, 0) == -1); - - rng.next_value = 0xff; - try expect(random.intRangeLessThan(i3, -4, 0) == -1); - rng.next_value = 0xff; - try expect(random.intRangeLessThan(i3, -2, 2) == 1); -} - -test "Random intAtMost" { - @setEvalBranchQuota(10000); - try testRandomIntAtMost(); - try comptime testRandomIntAtMost(); -} -fn testRandomIntAtMost() !void { - var rng = SequentialPrng.init(); - const random = rng.random(); - - rng.next_value = 0xff; - try expect(random.uintAtMost(u8, 3) == 3); - try expect(rng.next_value == 0); - try expect(random.uintAtMost(u8, 3) == 0); - - // trigger the bias rejection code path - rng.next_value = 0; - try expect(random.uintAtMost(u8, 2) == 0); - // verify we incremented twice - try expect(rng.next_value == 2); - - rng.next_value = 0xff; - try expect(random.intRangeAtMost(u8, 0, 0x7f) == 0x7f); - rng.next_value = 0xff; - try expect(random.intRangeAtMost(u8, 0x7f, 0xfe) == 0xfe); - - rng.next_value = 0xff; - try expect(random.intRangeAtMost(i8, 0, 0x3f) == 0x3f); - rng.next_value = 0xff; - try expect(random.intRangeAtMost(i8, -0x40, 0x3f) == 0x3f); - rng.next_value = 0xff; - try expect(random.intRangeAtMost(i8, -0x80, -1) == -1); - - rng.next_value = 0xff; - try expect(random.intRangeAtMost(i3, -4, -1) == -1); - rng.next_value = 0xff; - try expect(random.intRangeAtMost(i3, -2, 1) == 1); - - try expect(random.uintAtMost(u0, 0) == 0); -} - -test "Random Biased" { - var prng = DefaultPrng.init(0); - const random = prng.random(); - // Not thoroughly checking the logic here. - // Just want to execute all the paths with different types. - - try expect(random.uintLessThanBiased(u1, 1) == 0); - try expect(random.uintLessThanBiased(u32, 10) < 10); - try expect(random.uintLessThanBiased(u64, 20) < 20); - - try expect(random.uintAtMostBiased(u0, 0) == 0); - try expect(random.uintAtMostBiased(u1, 0) <= 0); - try expect(random.uintAtMostBiased(u32, 10) <= 10); - try expect(random.uintAtMostBiased(u64, 20) <= 20); - - try expect(random.intRangeLessThanBiased(u1, 0, 1) == 0); - try expect(random.intRangeLessThanBiased(i1, -1, 0) == -1); - try expect(random.intRangeLessThanBiased(u32, 10, 20) >= 10); - try expect(random.intRangeLessThanBiased(i32, 10, 20) >= 10); - try expect(random.intRangeLessThanBiased(u64, 20, 40) >= 20); - try expect(random.intRangeLessThanBiased(i64, 20, 40) >= 20); - - // uncomment for broken module error: - //expect(random.intRangeAtMostBiased(u0, 0, 0) == 0); - try expect(random.intRangeAtMostBiased(u1, 0, 1) >= 0); - try expect(random.intRangeAtMostBiased(i1, -1, 0) >= -1); - try expect(random.intRangeAtMostBiased(u32, 10, 20) >= 10); - try expect(random.intRangeAtMostBiased(i32, 10, 20) >= 10); - try expect(random.intRangeAtMostBiased(u64, 20, 40) >= 20); - try expect(random.intRangeAtMostBiased(i64, 20, 40) >= 20); -} - -test "splitmix64 sequence" { - var r = SplitMix64.init(0xaeecf86f7878dd75); - - const seq = [_]u64{ - 0x5dbd39db0178eb44, - 0xa9900fb66b397da3, - 0x5c1a28b1aeebcf5c, - 0x64a963238f776912, - 0xc6d4177b21d1c0ab, - 0xb2cbdbdb5ea35394, - }; - - for (seq) |s| { - try expect(s == r.next()); - } -} - -// Actual Random helper function tests, pcg engine is assumed correct. -test "Random float correctness" { - var prng = DefaultPrng.init(0); - const random = prng.random(); - - var i: usize = 0; - while (i < 1000) : (i += 1) { - const val1 = random.float(f32); - try expect(val1 >= 0.0); - try expect(val1 < 1.0); - - const val2 = random.float(f64); - try expect(val2 >= 0.0); - try expect(val2 < 1.0); - } -} - -// Check the "astronomically unlikely" code paths. -test "Random float coverage" { - var prng = try Dilbert.init(&[_]u8{0}); - const random = prng.random(); - - const rand_f64 = random.float(f64); - const rand_f32 = random.float(f32); - - try expect(rand_f32 == 0.0); - try expect(rand_f64 == 0.0); -} - -test "Random float chi-square goodness of fit" { - const num_numbers = 100000; - const num_buckets = 1000; - - var f32_hist = std.AutoHashMap(u32, u32).init(std.testing.allocator); - defer f32_hist.deinit(); - var f64_hist = std.AutoHashMap(u64, u32).init(std.testing.allocator); - defer f64_hist.deinit(); - - var prng = DefaultPrng.init(0); - const random = prng.random(); - - var i: usize = 0; - while (i < num_numbers) : (i += 1) { - const rand_f32 = random.float(f32); - const rand_f64 = random.float(f64); - const f32_put = try f32_hist.getOrPut(@as(u32, @intFromFloat(rand_f32 * @as(f32, @floatFromInt(num_buckets))))); - if (f32_put.found_existing) { - f32_put.value_ptr.* += 1; - } else { - f32_put.value_ptr.* = 1; - } - const f64_put = try f64_hist.getOrPut(@as(u32, @intFromFloat(rand_f64 * @as(f64, @floatFromInt(num_buckets))))); - if (f64_put.found_existing) { - f64_put.value_ptr.* += 1; - } else { - f64_put.value_ptr.* = 1; - } - } - - var f32_total_variance: f64 = 0; - var f64_total_variance: f64 = 0; - - { - var j: u32 = 0; - while (j < num_buckets) : (j += 1) { - const count = @as(f64, @floatFromInt((if (f32_hist.get(j)) |v| v else 0))); - const expected = @as(f64, @floatFromInt(num_numbers)) / @as(f64, @floatFromInt(num_buckets)); - const delta = count - expected; - const variance = (delta * delta) / expected; - f32_total_variance += variance; - } - } - - { - var j: u64 = 0; - while (j < num_buckets) : (j += 1) { - const count = @as(f64, @floatFromInt((if (f64_hist.get(j)) |v| v else 0))); - const expected = @as(f64, @floatFromInt(num_numbers)) / @as(f64, @floatFromInt(num_buckets)); - const delta = count - expected; - const variance = (delta * delta) / expected; - f64_total_variance += variance; - } - } - - // Accept p-values >= 0.05. - // Critical value is calculated by opening a Python interpreter and running: - // scipy.stats.chi2.isf(0.05, num_buckets - 1) - const critical_value = 1073.6426506574246; - try expect(f32_total_variance < critical_value); - try expect(f64_total_variance < critical_value); -} - -test "Random shuffle" { - var prng = DefaultPrng.init(0); - const random = prng.random(); - - var seq = [_]u8{ 0, 1, 2, 3, 4 }; - var seen = [_]bool{false} ** 5; - - var i: usize = 0; - while (i < 1000) : (i += 1) { - random.shuffle(u8, seq[0..]); - seen[seq[0]] = true; - try expect(sumArray(seq[0..]) == 10); - } - - // we should see every entry at the head at least once - for (seen) |e| { - try expect(e == true); - } -} - -fn sumArray(s: []const u8) u32 { - var r: u32 = 0; - for (s) |e| - r += e; - return r; -} - -test "Random range" { - var prng = DefaultPrng.init(0); - const random = prng.random(); - - try testRange(random, -4, 3); - try testRange(random, -4, -1); - try testRange(random, 10, 14); - try testRange(random, -0x80, 0x7f); -} - -fn testRange(r: Random, start: i8, end: i8) !void { - try testRangeBias(r, start, end, true); - try testRangeBias(r, start, end, false); -} -fn testRangeBias(r: Random, start: i8, end: i8, biased: bool) !void { - const count = @as(usize, @intCast(@as(i32, end) - @as(i32, start))); - var values_buffer = [_]bool{false} ** 0x100; - const values = values_buffer[0..count]; - var i: usize = 0; - while (i < count) { - const value: i32 = if (biased) r.intRangeLessThanBiased(i8, start, end) else r.intRangeLessThan(i8, start, end); - const index = @as(usize, @intCast(value - start)); - if (!values[index]) { - i += 1; - values[index] = true; - } - } -} - -test "CSPRNG" { - var secret_seed: [DefaultCsprng.secret_seed_length]u8 = undefined; - std.crypto.random.bytes(&secret_seed); - var csprng = DefaultCsprng.init(secret_seed); - const random = csprng.random(); - const a = random.int(u64); - const b = random.int(u64); - const c = random.int(u64); - try expect(a ^ b ^ c != 0); -} - -test "Random weightedIndex" { - // Make sure weightedIndex works for various integers and floats - inline for (.{ u64, i4, f32, f64 }) |T| { - var prng = DefaultPrng.init(0); - const random = prng.random(); - - const proportions = [_]T{ 2, 1, 1, 2 }; - var counts = [_]f64{ 0, 0, 0, 0 }; - - const n_trials: u64 = 10_000; - var i: usize = 0; - while (i < n_trials) : (i += 1) { - const pick = random.weightedIndex(T, &proportions); - counts[pick] += 1; - } - - // We expect the first and last counts to be roughly 2x the second and third - const approxEqRel = std.math.approxEqRel; - // Define "roughly" to be within 10% - const tolerance = 0.1; - try std.testing.expect(approxEqRel(f64, counts[0], counts[1] * 2, tolerance)); - try std.testing.expect(approxEqRel(f64, counts[1], counts[2], tolerance)); - try std.testing.expect(approxEqRel(f64, counts[2] * 2, counts[3], tolerance)); - } -} diff --git a/lib/std/rand/ziggurat.zig b/lib/std/rand/ziggurat.zig @@ -1,170 +0,0 @@ -//! Implements [ZIGNOR][1] (Jurgen A. Doornik, 2005, Nuffield College, Oxford). -//! -//! [1]: https://www.doornik.com/research/ziggurat.pdf -//! -//! rust/rand used as a reference; -//! -//! NOTE: This seems interesting but reference code is a bit hard to grok: -//! https://sbarral.github.io/etf. - -const std = @import("../std.zig"); -const builtin = @import("builtin"); -const math = std.math; -const Random = std.rand.Random; - -pub fn next_f64(random: Random, comptime tables: ZigTable) f64 { - while (true) { - // We manually construct a float from parts as we can avoid an extra random lookup here by - // using the unused exponent for the lookup table entry. - const bits = random.int(u64); - const i = @as(usize, @as(u8, @truncate(bits))); - - const u = blk: { - if (tables.is_symmetric) { - // Generate a value in the range [2, 4) and scale into [-1, 1) - const repr = ((0x3ff + 1) << 52) | (bits >> 12); - break :blk @as(f64, @bitCast(repr)) - 3.0; - } else { - // Generate a value in the range [1, 2) and scale into (0, 1) - const repr = (0x3ff << 52) | (bits >> 12); - break :blk @as(f64, @bitCast(repr)) - (1.0 - math.floatEps(f64) / 2.0); - } - }; - - const x = u * tables.x[i]; - const test_x = if (tables.is_symmetric) @abs(x) else x; - - // equivalent to |u| < tables.x[i+1] / tables.x[i] (or u < tables.x[i+1] / tables.x[i]) - if (test_x < tables.x[i + 1]) { - return x; - } - - if (i == 0) { - return tables.zero_case(random, u); - } - - // equivalent to f1 + DRanU() * (f0 - f1) < 1 - if (tables.f[i + 1] + (tables.f[i] - tables.f[i + 1]) * random.float(f64) < tables.pdf(x)) { - return x; - } - } -} - -pub const ZigTable = struct { - r: f64, - x: [257]f64, - f: [257]f64, - - // probability density function used as a fallback - pdf: fn (f64) f64, - // whether the distribution is symmetric - is_symmetric: bool, - // fallback calculation in the case we are in the 0 block - zero_case: fn (Random, f64) f64, -}; - -// zigNorInit -pub fn ZigTableGen( - comptime is_symmetric: bool, - comptime r: f64, - comptime v: f64, - comptime f: fn (f64) f64, - comptime f_inv: fn (f64) f64, - comptime zero_case: fn (Random, f64) f64, -) ZigTable { - var tables: ZigTable = undefined; - - tables.is_symmetric = is_symmetric; - tables.r = r; - tables.pdf = f; - tables.zero_case = zero_case; - - tables.x[0] = v / f(r); - tables.x[1] = r; - - for (tables.x[2..256], 0..) |*entry, i| { - const last = tables.x[2 + i - 1]; - entry.* = f_inv(v / last + f(last)); - } - tables.x[256] = 0; - - for (tables.f[0..], 0..) |*entry, i| { - entry.* = f(tables.x[i]); - } - - return tables; -} - -// N(0, 1) -pub const NormDist = blk: { - @setEvalBranchQuota(30000); - break :blk ZigTableGen(true, norm_r, norm_v, norm_f, norm_f_inv, norm_zero_case); -}; - -pub const norm_r = 3.6541528853610088; -pub const norm_v = 0.00492867323399; - -pub fn norm_f(x: f64) f64 { - return @exp(-x * x / 2.0); -} -pub fn norm_f_inv(y: f64) f64 { - return @sqrt(-2.0 * @log(y)); -} -pub fn norm_zero_case(random: Random, u: f64) f64 { - var x: f64 = 1; - var y: f64 = 0; - - while (-2.0 * y < x * x) { - x = @log(random.float(f64)) / norm_r; - y = @log(random.float(f64)); - } - - if (u < 0) { - return x - norm_r; - } else { - return norm_r - x; - } -} - -test "normal dist sanity" { - var prng = std.rand.DefaultPrng.init(0); - const random = prng.random(); - - var i: usize = 0; - while (i < 1000) : (i += 1) { - _ = random.floatNorm(f64); - } -} - -// Exp(1) -pub const ExpDist = blk: { - @setEvalBranchQuota(30000); - break :blk ZigTableGen(false, exp_r, exp_v, exp_f, exp_f_inv, exp_zero_case); -}; - -pub const exp_r = 7.69711747013104972; -pub const exp_v = 0.0039496598225815571993; - -pub fn exp_f(x: f64) f64 { - return @exp(-x); -} -pub fn exp_f_inv(y: f64) f64 { - return -@log(y); -} -pub fn exp_zero_case(random: Random, _: f64) f64 { - return exp_r - @log(random.float(f64)); -} - -test "exp dist smoke test" { - var prng = std.rand.DefaultPrng.init(0); - const random = prng.random(); - - var i: usize = 0; - while (i < 1000) : (i += 1) { - _ = random.floatExp(f64); - } -} - -test { - _ = NormDist; -} diff --git a/lib/std/sort.zig b/lib/std/sort.zig @@ -379,7 +379,7 @@ test "sort with context in the middle of a slice" { } test "sort fuzz testing" { - var prng = std.rand.DefaultPrng.init(0x12345678); + var prng = std.Random.DefaultPrng.init(0x12345678); const random = prng.random(); const test_case_count = 10; diff --git a/lib/std/std.zig b/lib/std/std.zig @@ -36,6 +36,7 @@ pub const PackedIntSliceEndian = @import("packed_int_array.zig").PackedIntSliceE pub const PriorityQueue = @import("priority_queue.zig").PriorityQueue; pub const PriorityDequeue = @import("priority_dequeue.zig").PriorityDequeue; pub const Progress = @import("Progress.zig"); +pub const Random = @import("Random.zig"); pub const RingBuffer = @import("RingBuffer.zig"); pub const SegmentedList = @import("segmented_list.zig").SegmentedList; pub const SemanticVersion = @import("SemanticVersion.zig"); @@ -156,8 +157,8 @@ pub const pdb = @import("pdb.zig"); /// and spawning of child processes. pub const process = @import("process.zig"); -/// Fast pseudo-random number generators (i.e. not cryptographically secure). -pub const rand = @import("rand.zig"); +/// Deprecated: use `Random` instead. +pub const rand = Random; /// Sorting. pub const sort = @import("sort.zig"); diff --git a/lib/std/treap.zig b/lib/std/treap.zig @@ -18,7 +18,7 @@ pub fn Treap(comptime Key: type, comptime compareFn: anytype) type { /// A customized pseudo random number generator for the treap. /// This just helps reducing the memory size of the treap itself - /// as std.rand.DefaultPrng requires larger state (while producing better entropy for randomness to be fair). + /// as std.Random.DefaultPrng requires larger state (while producing better entropy for randomness to be fair). const Prng = struct { xorshift: usize = 0, @@ -305,7 +305,7 @@ pub fn Treap(comptime Key: type, comptime compareFn: anytype) type { // https://lemire.me/blog/2017/09/18/visiting-all-values-in-an-array-exactly-once-in-random-order/ fn SliceIterRandomOrder(comptime T: type) type { return struct { - rng: std.rand.Random, + rng: std.Random, slice: []T, index: usize = undefined, offset: usize = undefined, @@ -313,7 +313,7 @@ fn SliceIterRandomOrder(comptime T: type) type { const Self = @This(); - pub fn init(slice: []T, rng: std.rand.Random) Self { + pub fn init(slice: []T, rng: std.Random) Self { return Self{ .rng = rng, .slice = slice, @@ -353,7 +353,7 @@ test "std.Treap: insert, find, replace, remove" { var treap = TestTreap{}; var nodes: [10]TestNode = undefined; - var prng = std.rand.DefaultPrng.init(0xdeadbeef); + var prng = std.Random.DefaultPrng.init(0xdeadbeef); var iter = SliceIterRandomOrder(TestNode).init(&nodes, prng.random()); // insert check diff --git a/src/reduce.zig b/src/reduce.zig @@ -136,7 +136,7 @@ pub fn main(gpa: Allocator, arena: Allocator, args: []const []const u8) !void { var more_fixups: Ast.Fixups = .{}; defer more_fixups.deinit(gpa); - var rng = std.rand.DefaultPrng.init(seed); + var rng = std.Random.DefaultPrng.init(seed); // 1. Walk the AST of the source file looking for independent // reductions and collecting them all into an array list. @@ -274,7 +274,7 @@ pub fn main(gpa: Allocator, arena: Allocator, args: []const []const u8) !void { return std.process.cleanExit(); } -fn sortTransformations(transformations: []Walk.Transformation, rng: std.rand.Random) void { +fn sortTransformations(transformations: []Walk.Transformation, rng: std.Random) void { rng.shuffle(Walk.Transformation, transformations); // Stable sort based on priority to keep randomness as the secondary sort. // TODO: introduce transformation priorities diff --git a/src/resinator/compile.zig b/src/resinator/compile.zig @@ -3356,7 +3356,7 @@ test "StringTable" { } break :ids buf; }; - var prng = std.rand.DefaultPrng.init(0); + var prng = std.Random.DefaultPrng.init(0); var random = prng.random(); random.shuffle(u16, &ids);