1//! The engines provided here should be initialized from an external source.
2//! Be sure to use a CSPRNG when required, otherwise using a normal PRNG will
3//! be faster and use substantially less stack space.
4const Random = @This();
5
6const std = @import("std.zig");
7const math = std.math;
8const mem = std.mem;
9const assert = std.debug.assert;
10const maxInt = std.math.maxInt;
11
12/// Fast unbiased random numbers.
13pub const DefaultPrng = Xoshiro256;
14
15/// Cryptographically secure random numbers.
16pub const DefaultCsprng = ChaCha;
17
18pub const Ascon = @import("Random/Ascon.zig");
19pub const ChaCha = @import("Random/ChaCha.zig");
20
21pub const Isaac64 = @import("Random/Isaac64.zig");
22pub const Pcg = @import("Random/Pcg.zig");
23pub const Xoroshiro128 = @import("Random/Xoroshiro128.zig");
24pub const Xoshiro256 = @import("Random/Xoshiro256.zig");
25pub const Sfc64 = @import("Random/Sfc64.zig");
26pub const RomuTrio = @import("Random/RomuTrio.zig");
27pub const SplitMix64 = @import("Random/SplitMix64.zig");
28pub const ziggurat = @import("Random/ziggurat.zig");
29pub const lcg = @import("Random/lcg.zig");
30
31/// Any comparison of this field may result in illegal behavior, since it may be set to
32/// `undefined` in cases where the random implementation does not have any associated
33/// state.
34ptr: *anyopaque,
35fillFn: *const fn (ptr: *anyopaque, buf: []u8) void,
36
37pub const IoSource = struct {
38 io: std.Io,
39
40 pub fn interface(this: *const @This()) std.Random {
41 return .{
42 .ptr = @constCast(this),
43 .fillFn = fill,
44 };
45 }
46
47 fn fill(ptr: *anyopaque, buffer: []u8) void {
48 const this: *const @This() = @ptrCast(@alignCast(ptr));
49 this.io.random(buffer);
50 }
51};
52
53pub fn init(pointer: anytype, comptime fillFn: fn (ptr: @TypeOf(pointer), buf: []u8) void) Random {
54 const Ptr = @TypeOf(pointer);
55 assert(@typeInfo(Ptr) == .pointer); // Must be a pointer
56 assert(@typeInfo(Ptr).pointer.size == .one); // Must be a single-item pointer
57 assert(@typeInfo(@typeInfo(Ptr).pointer.child) == .@"struct"); // Must point to a struct
58 const gen = struct {
59 fn fill(ptr: *anyopaque, buf: []u8) void {
60 const self: Ptr = @ptrCast(@alignCast(ptr));
61 fillFn(self, buf);
62 }
63 };
64
65 return .{
66 .ptr = pointer,
67 .fillFn = gen.fill,
68 };
69}
70
71/// Read random bytes into the specified buffer until full.
72pub fn bytes(r: Random, buf: []u8) void {
73 r.fillFn(r.ptr, buf);
74}
75
76pub fn array(r: Random, comptime E: type, comptime N: usize) [N]E {
77 var result: [N]E = undefined;
78 bytes(r, &result);
79 return result;
80}
81
82pub fn boolean(r: Random) bool {
83 return r.int(u1) != 0;
84}
85
86/// Returns a random value from an enum, evenly distributed.
87///
88/// Note that this will not yield consistent results across all targets
89/// due to dependence on the representation of `usize` as an index.
90/// See `enumValueWithIndex` for further commentary.
91pub inline fn enumValue(r: Random, comptime EnumType: type) EnumType {
92 return r.enumValueWithIndex(EnumType, usize);
93}
94
95/// Returns a random value from an enum, evenly distributed.
96///
97/// An index into an array of all named values is generated using the
98/// specified `Index` type to determine the return value.
99/// This allows for results to be independent of `usize` representation.
100///
101/// Prefer `enumValue` if this isn't important.
102///
103/// See `uintLessThan`, which this function uses in most cases,
104/// for commentary on the runtime of this function.
105pub fn enumValueWithIndex(r: Random, comptime EnumType: type, comptime Index: type) EnumType {
106 comptime assert(@typeInfo(EnumType) == .@"enum");
107
108 // We won't use int -> enum casting because enum elements can have
109 // arbitrary values. Instead we'll randomly pick one of the type's values.
110 const values = comptime std.enums.values(EnumType);
111 comptime assert(values.len > 0); // can't return anything
112 comptime assert(maxInt(Index) >= values.len - 1); // can't access all values
113 if (values.len == 1) return values[0];
114
115 const index = if (comptime values.len - 1 == maxInt(Index))
116 r.int(Index)
117 else
118 r.uintLessThan(Index, values.len);
119
120 const MinInt = MinArrayIndex(Index);
121 return values[@as(MinInt, @intCast(index))];
122}
123
124/// Returns a random int `i` such that `minInt(T) <= i <= maxInt(T)`.
125/// `i` is evenly distributed.
126pub fn int(r: Random, comptime T: type) T {
127 const bits = @typeInfo(T).int.bits;
128 const UnsignedT = @Int(.unsigned, bits);
129 const ceil_bytes = @divCeil(bits, 8);
130 const ByteAlignedT = @Int(.unsigned, ceil_bytes * 8);
131
132 var rand_bytes: [ceil_bytes]u8 = undefined;
133 r.bytes(&rand_bytes);
134
135 // use LE instead of native endian for better portability maybe?
136 // TODO: endian portability is pointless if the underlying prng isn't endian portable.
137 // TODO: document the endian portability of this library.
138 const byte_aligned_result = mem.readInt(ByteAlignedT, &rand_bytes, .little);
139 const unsigned_result: UnsignedT = @truncate(byte_aligned_result);
140 return @bitCast(unsigned_result);
141}
142
143/// Constant-time implementation off `uintLessThan`.
144/// The results of this function may be biased.
145pub fn uintLessThanBiased(r: Random, comptime T: type, less_than: T) T {
146 comptime assert(@typeInfo(T).int.signedness == .unsigned);
147 assert(0 < less_than);
148 return limitRangeBiased(T, r.int(T), less_than);
149}
150
151/// Returns an evenly distributed random unsigned integer `0 <= i < less_than`.
152/// This function assumes that the underlying `fillFn` produces evenly distributed values.
153/// Within this assumption, the runtime of this function is exponentially distributed.
154/// If `fillFn` were backed by a true random generator,
155/// the runtime of this function would technically be unbounded.
156/// However, if `fillFn` is backed by any evenly distributed pseudo random number generator,
157/// this function is guaranteed to return.
158/// If you need deterministic runtime bounds, use `uintLessThanBiased`.
159pub fn uintLessThan(r: Random, comptime T: type, less_than: T) T {
160 comptime assert(@typeInfo(T).int.signedness == .unsigned);
161 const bits = @typeInfo(T).int.bits;
162 assert(0 < less_than);
163
164 // adapted from:
165 // http://www.pcg-random.org/posts/bounded-rands.html
166 // "Lemire's (with an extra tweak from me)"
167 var x = r.int(T);
168 var m = math.mulWide(T, x, less_than);
169 var l: T = @truncate(m);
170 if (l < less_than) {
171 var t = -%less_than;
172
173 if (t >= less_than) {
174 t -= less_than;
175 if (t >= less_than) {
176 t %= less_than;
177 }
178 }
179 while (l < t) {
180 x = r.int(T);
181 m = math.mulWide(T, x, less_than);
182 l = @truncate(m);
183 }
184 }
185 return @intCast(m >> bits);
186}
187
188/// Constant-time implementation off `uintAtMost`.
189/// The results of this function may be biased.
190pub fn uintAtMostBiased(r: Random, comptime T: type, at_most: T) T {
191 assert(@typeInfo(T).int.signedness == .unsigned);
192 if (at_most == maxInt(T)) {
193 // have the full range
194 return r.int(T);
195 }
196 return r.uintLessThanBiased(T, at_most + 1);
197}
198
199/// Returns an evenly distributed random unsigned integer `0 <= i <= at_most`.
200/// See `uintLessThan`, which this function uses in most cases,
201/// for commentary on the runtime of this function.
202pub fn uintAtMost(r: Random, comptime T: type, at_most: T) T {
203 assert(@typeInfo(T).int.signedness == .unsigned);
204 if (at_most == maxInt(T)) {
205 // have the full range
206 return r.int(T);
207 }
208 return r.uintLessThan(T, at_most + 1);
209}
210
211/// Constant-time implementation off `intRangeLessThan`.
212/// The results of this function may be biased.
213pub fn intRangeLessThanBiased(r: Random, comptime T: type, at_least: T, less_than: T) T {
214 assert(at_least < less_than);
215 const info = @typeInfo(T).int;
216 if (info.signedness == .signed) {
217 // Two's complement makes this math pretty easy.
218 const UnsignedT = @Int(.unsigned, info.bits);
219 const lo: UnsignedT = @bitCast(at_least);
220 const hi: UnsignedT = @bitCast(less_than);
221 const result = lo +% r.uintLessThanBiased(UnsignedT, hi -% lo);
222 return @bitCast(result);
223 } else {
224 // The signed implementation would work fine, but we can use stricter arithmetic operators here.
225 return at_least + r.uintLessThanBiased(T, less_than - at_least);
226 }
227}
228
229/// Returns an evenly distributed random integer `at_least <= i < less_than`.
230/// See `uintLessThan`, which this function uses in most cases,
231/// for commentary on the runtime of this function.
232pub fn intRangeLessThan(r: Random, comptime T: type, at_least: T, less_than: T) T {
233 assert(at_least < less_than);
234 const info = @typeInfo(T).int;
235 if (info.signedness == .signed) {
236 // Two's complement makes this math pretty easy.
237 const UnsignedT = @Int(.unsigned, info.bits);
238 const lo: UnsignedT = @bitCast(at_least);
239 const hi: UnsignedT = @bitCast(less_than);
240 const result = lo +% r.uintLessThan(UnsignedT, hi -% lo);
241 return @bitCast(result);
242 } else {
243 // The signed implementation would work fine, but we can use stricter arithmetic operators here.
244 return at_least + r.uintLessThan(T, less_than - at_least);
245 }
246}
247
248/// Constant-time implementation off `intRangeAtMostBiased`.
249/// The results of this function may be biased.
250pub fn intRangeAtMostBiased(r: Random, comptime T: type, at_least: T, at_most: T) T {
251 assert(at_least <= at_most);
252 const info = @typeInfo(T).int;
253 if (info.signedness == .signed) {
254 // Two's complement makes this math pretty easy.
255 const UnsignedT = @Int(.unsigned, info.bits);
256 const lo: UnsignedT = @bitCast(at_least);
257 const hi: UnsignedT = @bitCast(at_most);
258 const result = lo +% r.uintAtMostBiased(UnsignedT, hi -% lo);
259 return @bitCast(result);
260 } else {
261 // The signed implementation would work fine, but we can use stricter arithmetic operators here.
262 return at_least + r.uintAtMostBiased(T, at_most - at_least);
263 }
264}
265
266/// Returns an evenly distributed random integer `at_least <= i <= at_most`.
267/// See `uintLessThan`, which this function uses in most cases,
268/// for commentary on the runtime of this function.
269pub fn intRangeAtMost(r: Random, comptime T: type, at_least: T, at_most: T) T {
270 assert(at_least <= at_most);
271 const info = @typeInfo(T).int;
272 if (info.signedness == .signed) {
273 // Two's complement makes this math pretty easy.
274 const UnsignedT = @Int(.unsigned, info.bits);
275 const lo: UnsignedT = @bitCast(at_least);
276 const hi: UnsignedT = @bitCast(at_most);
277 const result = lo +% r.uintAtMost(UnsignedT, hi -% lo);
278 return @bitCast(result);
279 } else {
280 // The signed implementation would work fine, but we can use stricter arithmetic operators here.
281 return at_least + r.uintAtMost(T, at_most - at_least);
282 }
283}
284
285/// Return a floating point value evenly distributed in the range [0, 1).
286pub fn float(r: Random, comptime T: type) T {
287 // Generate a uniformly random value for the mantissa.
288 // Then generate an exponentially biased random value for the exponent.
289 // This covers every possible value in the range.
290 switch (T) {
291 f32 => {
292 // Use 23 random bits for the mantissa, and the rest for the exponent.
293 // If all 41 bits are zero, generate additional random bits, until a
294 // set bit is found, or 126 bits have been generated.
295 const rand = r.int(u64);
296 var rand_lz = @clz(rand);
297 if (rand_lz >= 41) {
298 @branchHint(.unlikely);
299 rand_lz = 41 + @clz(r.int(u64));
300 if (rand_lz == 41 + 64) {
301 @branchHint(.unlikely);
302 // It is astronomically unlikely to reach this point.
303 rand_lz += @clz(r.int(u32) | 0x7FF);
304 }
305 }
306 const mantissa: u23 = @truncate(rand);
307 const exponent = @as(u32, 126 - rand_lz) << 23;
308 return @bitCast(exponent | mantissa);
309 },
310 f64 => {
311 // Use 52 random bits for the mantissa, and the rest for the exponent.
312 // If all 12 bits are zero, generate additional random bits, until a
313 // set bit is found, or 1022 bits have been generated.
314 const rand = r.int(u64);
315 var rand_lz: u64 = @clz(rand);
316 if (rand_lz >= 12) {
317 rand_lz = 12;
318 while (true) {
319 // It is astronomically unlikely for this loop to execute more than once.
320 const addl_rand_lz = @clz(r.int(u64));
321 rand_lz += addl_rand_lz;
322 if (addl_rand_lz != 64) {
323 @branchHint(.likely);
324 break;
325 }
326 if (rand_lz >= 1022) {
327 rand_lz = 1022;
328 break;
329 }
330 }
331 }
332 const mantissa = rand & 0xFFFFFFFFFFFFF;
333 const exponent = (1022 - rand_lz) << 52;
334 return @bitCast(exponent | mantissa);
335 },
336 else => @compileError("unknown floating point type"),
337 }
338}
339
340/// Return a floating point value normally distributed with mean = 0, stddev = 1.
341///
342/// To use different parameters, use: floatNorm(...) * desiredStddev + desiredMean.
343pub fn floatNorm(r: Random, comptime T: type) T {
344 const value = ziggurat.next_f64(r, ziggurat.NormDist);
345 switch (T) {
346 f32 => return @floatCast(value),
347 f64 => return value,
348 else => @compileError("unknown floating point type"),
349 }
350}
351
352/// Return an exponentially distributed float with a rate parameter of 1.
353///
354/// To use a different rate parameter, use: floatExp(...) / desiredRate.
355pub fn floatExp(r: Random, comptime T: type) T {
356 const value = ziggurat.next_f64(r, ziggurat.ExpDist);
357 switch (T) {
358 f32 => return @floatCast(value),
359 f64 => return value,
360 else => @compileError("unknown floating point type"),
361 }
362}
363
364/// Shuffle a slice into a random order.
365///
366/// Note that this will not yield consistent results across all targets
367/// due to dependence on the representation of `usize` as an index.
368/// See `shuffleWithIndex` for further commentary.
369pub inline fn shuffle(r: Random, comptime T: type, buf: []T) void {
370 r.shuffleWithIndex(T, buf, usize);
371}
372
373/// Shuffle a slice into a random order, using an index of a
374/// specified type to maintain distribution across targets.
375/// Asserts the index type can represent `buf.len`.
376///
377/// Indexes into the slice are generated using the specified `Index`
378/// type, which determines distribution properties. This allows for
379/// results to be independent of `usize` representation.
380///
381/// Prefer `shuffle` if this isn't important.
382///
383/// See `intRangeLessThan`, which this function uses,
384/// for commentary on the runtime of this function.
385pub fn shuffleWithIndex(r: Random, comptime T: type, buf: []T, comptime Index: type) void {
386 const MinInt = MinArrayIndex(Index);
387 if (buf.len < 2) {
388 return;
389 }
390
391 // `i <= j < max <= maxInt(MinInt)`
392 const max: MinInt = @intCast(buf.len);
393 var i: MinInt = 0;
394 while (i < max - 1) : (i += 1) {
395 const j: MinInt = @intCast(r.intRangeLessThan(Index, i, max));
396 mem.swap(T, &buf[i], &buf[j]);
397 }
398}
399
400/// Randomly selects an index into `proportions`, where the likelihood of each
401/// index is weighted by that proportion.
402/// It is more likely for the index of the last proportion to be returned
403/// than the index of the first proportion in the slice, and vice versa.
404///
405/// This is useful for selecting an item from a slice where weights are not equal.
406/// `T` must be a numeric type capable of holding the sum of `proportions`.
407pub fn weightedIndex(r: Random, comptime T: type, proportions: []const T) usize {
408 // This implementation works by summing the proportions and picking a
409 // random point in [0, sum). We then loop over the proportions,
410 // accumulating until our accumulator is greater than the random point.
411
412 const sum = s: {
413 var sum: T = 0;
414 for (proportions) |v| sum += v;
415 break :s sum;
416 };
417
418 const point = switch (@typeInfo(T)) {
419 .int => |int_info| switch (int_info.signedness) {
420 .signed => r.intRangeLessThan(T, 0, sum),
421 .unsigned => r.uintLessThan(T, sum),
422 },
423 // take care that imprecision doesn't lead to a value slightly greater than sum
424 .float => @min(r.float(T) * sum, sum - std.math.floatEps(T)),
425 else => @compileError("weightedIndex does not support proportions of type " ++
426 @typeName(T)),
427 };
428
429 assert(point < sum);
430
431 var accumulator: T = 0;
432 for (proportions, 0..) |p, index| {
433 accumulator += p;
434 if (point < accumulator) return index;
435 } else unreachable;
436}
437
438/// Convert a random integer 0 <= random_int <= maxValue(T),
439/// into an integer 0 <= result < less_than.
440/// This function introduces a minor bias.
441pub fn limitRangeBiased(comptime T: type, random_int: T, less_than: T) T {
442 comptime assert(@typeInfo(T).int.signedness == .unsigned);
443 const bits = @typeInfo(T).int.bits;
444
445 // adapted from:
446 // http://www.pcg-random.org/posts/bounded-rands.html
447 // "Integer Multiplication (Biased)"
448 const m = math.mulWide(T, random_int, less_than);
449 return @intCast(m >> bits);
450}
451
452/// Returns the smallest of `Index` and `usize`.
453fn MinArrayIndex(comptime Index: type) type {
454 const index_info = @typeInfo(Index).int;
455 assert(index_info.signedness == .unsigned);
456 return if (index_info.bits >= @typeInfo(usize).int.bits) usize else Index;
457}
458
459test {
460 std.testing.refAllDecls(@This());
461 _ = @import("Random/test.zig");
462}