authorgravatar for thejoshwolfe@gmail.comJosh Wolfe <thejoshwolfe@gmail.com> 2018-11-21 17:33:37-05:00
committergravatar for thejoshwolfe@gmail.comJosh Wolfe <thejoshwolfe@gmail.com> 2018-11-24 22:25:21-05:00
log1924ffa67d90f21852680b7c9a09df07fc218ebe
tree28606b05dffdd654ee48e3e04608ba1a10c2e854
parent4d747d452fba0987332244c777dfd3a29d5a20bd

better debiased random range implementation

1 files changed, 46 insertions(+), 32 deletions(-)

std/rand/index.zig+46-32
...@@ -69,22 +69,36 @@ pub const Random = struct {...@@ -69,22 +69,36 @@ pub const Random = struct {
69 pub fn uintLessThan(r: *Random, comptime T: type, less_than: T) T {69 pub fn uintLessThan(r: *Random, comptime T: type, less_than: T) T {
70 assert(T.is_signed == false);70 assert(T.is_signed == false);
71 assert(0 < less_than);71 assert(0 < less_than);
7272 // Small is typically u32
73 const last_group_size_minus_one: T = maxInt(T) % less_than;73 const Small = @IntType(false, @divTrunc(T.bit_count + 31, 32) * 32);
74 if (last_group_size_minus_one == less_than - 1) {74 // Large is typically u64
75 // less_than is a power of two.75 const Large = @IntType(false, Small.bit_count * 2);
76 assert(math.floorPowerOfTwo(T, less_than) == less_than);76
77 // There is no retry zone. The optimal retry_zone_start would be maxInt(T) + 1.77 // adapted from:
78 return r.int(T) % less_than;78 // http://www.pcg-random.org/posts/bounded-rands.html
79 }79 // "Lemire's (with an extra tweak from me)"
80 const retry_zone_start = maxInt(T) - last_group_size_minus_one;80 var x: Small = r.int(Small);
8181 var m: Large = Large(x) * Large(less_than);
82 while (true) {82 var l: Small = @truncate(Small, m);
83 const rand_val = r.int(T);83 if (l < less_than) {
84 if (rand_val < retry_zone_start) {84 // TODO: workaround for https://github.com/ziglang/zig/issues/1770
85 return rand_val % less_than;85 // should be:
86 // var t: Small = -%less_than;
87 var t: Small = @bitCast(Small, -%@bitCast(@IntType(true, Small.bit_count), Small(less_than)));
88
89 if (t >= less_than) {
90 t -= less_than;
91 if (t >= less_than) {
92 t %= less_than;
93 }
94 }
95 while (l < t) {
96 x = r.int(Small);
97 m = Large(x) * Large(less_than);
98 l = @truncate(Small, m);
86 }99 }
87 }100 }
101 return @intCast(T, m >> Small.bit_count);
88 }102 }
89103
90 /// Returns an evenly distributed random unsigned integer `0 <= i <= at_most`.104 /// Returns an evenly distributed random unsigned integer `0 <= i <= at_most`.
...@@ -294,10 +308,19 @@ fn testRandomIntLessThan() void {...@@ -294,10 +308,19 @@ fn testRandomIntLessThan() void {
294 var r = SequentialPrng.init();308 var r = SequentialPrng.init();
295 r.next_value = 0xff;309 r.next_value = 0xff;
296 assert(r.random.uintLessThan(u8, 4) == 3);310 assert(r.random.uintLessThan(u8, 4) == 3);
297 r.next_value = 0xff;311 assert(r.next_value == 0);
298 assert(r.random.uintLessThan(u8, 3) == 0);312 assert(r.random.uintLessThan(u8, 4) == 0);
299 assert(r.next_value == 1);313 assert(r.next_value == 1);
300314
315 r.next_value = 0;
316 assert(r.random.uintLessThan(u64, 32) == 0);
317
318 // trigger the bias rejection code path
319 r.next_value = 0;
320 assert(r.random.uintLessThan(u8, 3) == 0);
321 // verify we incremented twice
322 assert(r.next_value == 2);
323
301 r.next_value = 0xff;324 r.next_value = 0xff;
302 assert(r.random.intRangeLessThan(u8, 0, 0x80) == 0x7f);325 assert(r.random.intRangeLessThan(u8, 0, 0x80) == 0x7f);
303 r.next_value = 0xff;326 r.next_value = 0xff;
...@@ -310,17 +333,10 @@ fn testRandomIntLessThan() void {...@@ -310,17 +333,10 @@ fn testRandomIntLessThan() void {
310 r.next_value = 0xff;333 r.next_value = 0xff;
311 assert(r.random.intRangeLessThan(i8, -0x80, 0) == -1);334 assert(r.random.intRangeLessThan(i8, -0x80, 0) == -1);
312335
313 r.next_value = 0xff;
314 assert(r.random.intRangeLessThan(i64, -0x8000000000000000, 0) == -1);
315 r.next_value = 0xff;336 r.next_value = 0xff;
316 assert(r.random.intRangeLessThan(i3, -4, 0) == -1);337 assert(r.random.intRangeLessThan(i3, -4, 0) == -1);
317 r.next_value = 0xff;338 r.next_value = 0xff;
318 assert(r.random.intRangeLessThan(i3, -2, 2) == 1);339 assert(r.random.intRangeLessThan(i3, -2, 2) == 1);
319
320 // test retrying and eventually getting a good value
321 // start just out of bounds
322 r.next_value = 0x81;
323 assert(r.random.uintLessThan(u8, 0x81) == 0);
324}340}
325341
326test "Random intAtMost" {342test "Random intAtMost" {
...@@ -332,9 +348,14 @@ fn testRandomIntAtMost() void {...@@ -332,9 +348,14 @@ fn testRandomIntAtMost() void {
332 var r = SequentialPrng.init();348 var r = SequentialPrng.init();
333 r.next_value = 0xff;349 r.next_value = 0xff;
334 assert(r.random.uintAtMost(u8, 3) == 3);350 assert(r.random.uintAtMost(u8, 3) == 3);
335 r.next_value = 0xff;351 assert(r.next_value == 0);
352 assert(r.random.uintAtMost(u8, 3) == 0);
353
354 // trigger the bias rejection code path
355 r.next_value = 0;
336 assert(r.random.uintAtMost(u8, 2) == 0);356 assert(r.random.uintAtMost(u8, 2) == 0);
337 assert(r.next_value == 1);357 // verify we incremented twice
358 assert(r.next_value == 2);
338359
339 r.next_value = 0xff;360 r.next_value = 0xff;
340 assert(r.random.intRangeAtMost(u8, 0, 0x7f) == 0x7f);361 assert(r.random.intRangeAtMost(u8, 0, 0x7f) == 0x7f);
...@@ -348,17 +369,10 @@ fn testRandomIntAtMost() void {...@@ -348,17 +369,10 @@ fn testRandomIntAtMost() void {
348 r.next_value = 0xff;369 r.next_value = 0xff;
349 assert(r.random.intRangeAtMost(i8, -0x80, -1) == -1);370 assert(r.random.intRangeAtMost(i8, -0x80, -1) == -1);
350371
351 r.next_value = 0xff;
352 assert(r.random.intRangeAtMost(i64, -0x8000000000000000, -1) == -1);
353 r.next_value = 0xff;372 r.next_value = 0xff;
354 assert(r.random.intRangeAtMost(i3, -4, -1) == -1);373 assert(r.random.intRangeAtMost(i3, -4, -1) == -1);
355 r.next_value = 0xff;374 r.next_value = 0xff;
356 assert(r.random.intRangeAtMost(i3, -2, 1) == 1);375 assert(r.random.intRangeAtMost(i3, -2, 1) == 1);
357
358 // test retrying and eventually getting a good value
359 // start just out of bounds
360 r.next_value = 0x81;
361 assert(r.random.uintAtMost(u8, 0x80) == 0);
362}376}
363377
364// Generator to extend 64-bit seed values into longer sequences.378// Generator to extend 64-bit seed values into longer sequences.