authorgravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2018-11-25 11:44:08-05:00
committergravatar for noreply@github.comGitHub <noreply@github.com> 2018-11-25 11:44:08-05:00
loga07490a4bc501334727e70b9dca387d03c24fe2e
tree0986fb8114e0597c8db6811244a17f2df7ada4c1
parent4d747d452fba0987332244c777dfd3a29d5a20bd
parent9ae5200bd2851aab704b786f53a1495a9f58049e
signaturebadge-question-mark Signed by PGP key 4AEE18F83AFDEB23

Merge pull request #1783 from ziglang/rand-range

Use better rand range implementations

1 files changed, 155 insertions(+), 40 deletions(-)

std/rand/index.zig+155-40
...@@ -57,6 +57,18 @@ pub const Random = struct {...@@ -57,6 +57,18 @@ pub const Random = struct {
57 return @bitCast(T, unsigned_result);57 return @bitCast(T, unsigned_result);
58 }58 }
5959
60 /// Constant-time implementation off ::uintLessThan.
61 /// The results of this function may be biased.
62 pub fn uintLessThanBiased(r: *Random, comptime T: type, less_than: T) T {
63 comptime assert(T.is_signed == false);
64 comptime assert(T.bit_count <= 64); // TODO: workaround: LLVM ERROR: Unsupported library call operation!
65 assert(0 < less_than);
66 if (T.bit_count <= 32) {
67 return @intCast(T, limitRangeBiased(u32, r.int(u32), less_than));
68 } else {
69 return @intCast(T, limitRangeBiased(u64, r.int(u64), less_than));
70 }
71 }
60 /// Returns an evenly distributed random unsigned integer `0 <= i < less_than`.72 /// Returns an evenly distributed random unsigned integer `0 <= i < less_than`.
61 /// This function assumes that the underlying ::fillFn produces evenly distributed values.73 /// This function assumes that the underlying ::fillFn produces evenly distributed values.
62 /// Within this assumption, the runtime of this function is exponentially distributed.74 /// Within this assumption, the runtime of this function is exponentially distributed.
...@@ -64,29 +76,53 @@ pub const Random = struct {...@@ -64,29 +76,53 @@ pub const Random = struct {
64 /// the runtime of this function would technically be unbounded.76 /// the runtime of this function would technically be unbounded.
65 /// However, if ::fillFn is backed by any evenly distributed pseudo random number generator,77 /// However, if ::fillFn is backed by any evenly distributed pseudo random number generator,
66 /// this function is guaranteed to return.78 /// this function is guaranteed to return.
67 /// If you need deterministic runtime bounds, consider instead using `r.int(T) % less_than`,79 /// If you need deterministic runtime bounds, use `::uintLessThanBiased`.
68 /// which will usually be biased toward smaller values.
69 pub fn uintLessThan(r: *Random, comptime T: type, less_than: T) T {80 pub fn uintLessThan(r: *Random, comptime T: type, less_than: T) T {
70 assert(T.is_signed == false);81 comptime assert(T.is_signed == false);
82 comptime assert(T.bit_count <= 64); // TODO: workaround: LLVM ERROR: Unsupported library call operation!
71 assert(0 < less_than);83 assert(0 < less_than);
7284 // Small is typically u32
73 const last_group_size_minus_one: T = maxInt(T) % less_than;85 const Small = @IntType(false, @divTrunc(T.bit_count + 31, 32) * 32);
74 if (last_group_size_minus_one == less_than - 1) {86 // Large is typically u64
75 // less_than is a power of two.87 const Large = @IntType(false, Small.bit_count * 2);
76 assert(math.floorPowerOfTwo(T, less_than) == less_than);88
77 // There is no retry zone. The optimal retry_zone_start would be maxInt(T) + 1.89 // adapted from:
78 return r.int(T) % less_than;90 // http://www.pcg-random.org/posts/bounded-rands.html
79 }91 // "Lemire's (with an extra tweak from me)"
80 const retry_zone_start = maxInt(T) - last_group_size_minus_one;92 var x: Small = r.int(Small);
8193 var m: Large = Large(x) * Large(less_than);
82 while (true) {94 var l: Small = @truncate(Small, m);
83 const rand_val = r.int(T);95 if (l < less_than) {
84 if (rand_val < retry_zone_start) {96 // TODO: workaround for https://github.com/ziglang/zig/issues/1770
85 return rand_val % less_than;97 // should be:
98 // var t: Small = -%less_than;
99 var t: Small = @bitCast(Small, -%@bitCast(@IntType(true, Small.bit_count), Small(less_than)));
100
101 if (t >= less_than) {
102 t -= less_than;
103 if (t >= less_than) {
104 t %= less_than;
105 }
106 }
107 while (l < t) {
108 x = r.int(Small);
109 m = Large(x) * Large(less_than);
110 l = @truncate(Small, m);
86 }111 }
87 }112 }
113 return @intCast(T, m >> Small.bit_count);
88 }114 }
89115
116 /// Constant-time implementation off ::uintAtMost.
117 /// The results of this function may be biased.
118 pub fn uintAtMostBiased(r: *Random, comptime T: type, at_most: T) T {
119 assert(T.is_signed == false);
120 if (at_most == maxInt(T)) {
121 // have the full range
122 return r.int(T);
123 }
124 return r.uintLessThanBiased(T, at_most + 1);
125 }
90 /// Returns an evenly distributed random unsigned integer `0 <= i <= at_most`.126 /// Returns an evenly distributed random unsigned integer `0 <= i <= at_most`.
91 /// See ::uintLessThan, which this function uses in most cases,127 /// See ::uintLessThan, which this function uses in most cases,
92 /// for commentary on the runtime of this function.128 /// for commentary on the runtime of this function.
...@@ -99,6 +135,22 @@ pub const Random = struct {...@@ -99,6 +135,22 @@ pub const Random = struct {
99 return r.uintLessThan(T, at_most + 1);135 return r.uintLessThan(T, at_most + 1);
100 }136 }
101137
138 /// Constant-time implementation off ::intRangeLessThan.
139 /// The results of this function may be biased.
140 pub fn intRangeLessThanBiased(r: *Random, comptime T: type, at_least: T, less_than: T) T {
141 assert(at_least < less_than);
142 if (T.is_signed) {
143 // Two's complement makes this math pretty easy.
144 const UnsignedT = @IntType(false, T.bit_count);
145 const lo = @bitCast(UnsignedT, at_least);
146 const hi = @bitCast(UnsignedT, less_than);
147 const result = lo +% r.uintLessThanBiased(UnsignedT, hi -% lo);
148 return @bitCast(T, result);
149 } else {
150 // The signed implementation would work fine, but we can use stricter arithmetic operators here.
151 return at_least + r.uintLessThanBiased(T, less_than - at_least);
152 }
153 }
102 /// Returns an evenly distributed random integer `at_least <= i < less_than`.154 /// Returns an evenly distributed random integer `at_least <= i < less_than`.
103 /// See ::uintLessThan, which this function uses in most cases,155 /// See ::uintLessThan, which this function uses in most cases,
104 /// for commentary on the runtime of this function.156 /// for commentary on the runtime of this function.
...@@ -117,6 +169,22 @@ pub const Random = struct {...@@ -117,6 +169,22 @@ pub const Random = struct {
117 }169 }
118 }170 }
119171
172 /// Constant-time implementation off ::intRangeAtMostBiased.
173 /// The results of this function may be biased.
174 pub fn intRangeAtMostBiased(r: *Random, comptime T: type, at_least: T, at_most: T) T {
175 assert(at_least <= at_most);
176 if (T.is_signed) {
177 // Two's complement makes this math pretty easy.
178 const UnsignedT = @IntType(false, T.bit_count);
179 const lo = @bitCast(UnsignedT, at_least);
180 const hi = @bitCast(UnsignedT, at_most);
181 const result = lo +% r.uintAtMostBiased(UnsignedT, hi -% lo);
182 return @bitCast(T, result);
183 } else {
184 // The signed implementation would work fine, but we can use stricter arithmetic operators here.
185 return at_least + r.uintAtMostBiased(T, at_most - at_least);
186 }
187 }
120 /// Returns an evenly distributed random integer `at_least <= i <= at_most`.188 /// Returns an evenly distributed random integer `at_least <= i <= at_most`.
121 /// See ::uintLessThan, which this function uses in most cases,189 /// See ::uintLessThan, which this function uses in most cases,
122 /// for commentary on the runtime of this function.190 /// for commentary on the runtime of this function.
...@@ -135,15 +203,11 @@ pub const Random = struct {...@@ -135,15 +203,11 @@ pub const Random = struct {
135 }203 }
136 }204 }
137205
138 /// Return a random integer/boolean type.
139 /// TODO: deprecated. use ::boolean or ::int instead.206 /// TODO: deprecated. use ::boolean or ::int instead.
140 pub fn scalar(r: *Random, comptime T: type) T {207 pub fn scalar(r: *Random, comptime T: type) T {
141 if (T == bool) return r.boolean();208 return if (T == bool) r.boolean() else r.int(T);
142 return r.int(T);
143 }209 }
144210
145 /// Return a random integer with even distribution between `start`
146 /// inclusive and `end` exclusive. `start` must be less than `end`.
147 /// TODO: deprecated. renamed to ::intRangeLessThan211 /// TODO: deprecated. renamed to ::intRangeLessThan
148 pub fn range(r: *Random, comptime T: type, start: T, end: T) T {212 pub fn range(r: *Random, comptime T: type, start: T, end: T) T {
149 return r.intRangeLessThan(T, start, end);213 return r.intRangeLessThan(T, start, end);
...@@ -206,6 +270,20 @@ pub const Random = struct {...@@ -206,6 +270,20 @@ pub const Random = struct {
206 }270 }
207};271};
208272
273/// Convert a random integer 0 <= random_int <= maxValue(T),
274/// into an integer 0 <= result < less_than.
275/// This function introduces a minor bias.
276pub fn limitRangeBiased(comptime T: type, random_int: T, less_than: T) T {
277 comptime assert(T.is_signed == false);
278 const T2 = @IntType(false, T.bit_count * 2);
279
280 // adapted from:
281 // http://www.pcg-random.org/posts/bounded-rands.html
282 // "Integer Multiplication (Biased)"
283 var m: T2 = T2(random_int) * T2(less_than);
284 return @intCast(T, m >> T.bit_count);
285}
286
209const SequentialPrng = struct {287const SequentialPrng = struct {
210 const Self = @This();288 const Self = @This();
211 random: Random,289 random: Random,
...@@ -294,10 +372,19 @@ fn testRandomIntLessThan() void {...@@ -294,10 +372,19 @@ fn testRandomIntLessThan() void {
294 var r = SequentialPrng.init();372 var r = SequentialPrng.init();
295 r.next_value = 0xff;373 r.next_value = 0xff;
296 assert(r.random.uintLessThan(u8, 4) == 3);374 assert(r.random.uintLessThan(u8, 4) == 3);
297 r.next_value = 0xff;375 assert(r.next_value == 0);
298 assert(r.random.uintLessThan(u8, 3) == 0);376 assert(r.random.uintLessThan(u8, 4) == 0);
299 assert(r.next_value == 1);377 assert(r.next_value == 1);
300378
379 r.next_value = 0;
380 assert(r.random.uintLessThan(u64, 32) == 0);
381
382 // trigger the bias rejection code path
383 r.next_value = 0;
384 assert(r.random.uintLessThan(u8, 3) == 0);
385 // verify we incremented twice
386 assert(r.next_value == 2);
387
301 r.next_value = 0xff;388 r.next_value = 0xff;
302 assert(r.random.intRangeLessThan(u8, 0, 0x80) == 0x7f);389 assert(r.random.intRangeLessThan(u8, 0, 0x80) == 0x7f);
303 r.next_value = 0xff;390 r.next_value = 0xff;
...@@ -310,17 +397,10 @@ fn testRandomIntLessThan() void {...@@ -310,17 +397,10 @@ fn testRandomIntLessThan() void {
310 r.next_value = 0xff;397 r.next_value = 0xff;
311 assert(r.random.intRangeLessThan(i8, -0x80, 0) == -1);398 assert(r.random.intRangeLessThan(i8, -0x80, 0) == -1);
312399
313 r.next_value = 0xff;
314 assert(r.random.intRangeLessThan(i64, -0x8000000000000000, 0) == -1);
315 r.next_value = 0xff;400 r.next_value = 0xff;
316 assert(r.random.intRangeLessThan(i3, -4, 0) == -1);401 assert(r.random.intRangeLessThan(i3, -4, 0) == -1);
317 r.next_value = 0xff;402 r.next_value = 0xff;
318 assert(r.random.intRangeLessThan(i3, -2, 2) == 1);403 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}404}
325405
326test "Random intAtMost" {406test "Random intAtMost" {
...@@ -332,9 +412,14 @@ fn testRandomIntAtMost() void {...@@ -332,9 +412,14 @@ fn testRandomIntAtMost() void {
332 var r = SequentialPrng.init();412 var r = SequentialPrng.init();
333 r.next_value = 0xff;413 r.next_value = 0xff;
334 assert(r.random.uintAtMost(u8, 3) == 3);414 assert(r.random.uintAtMost(u8, 3) == 3);
335 r.next_value = 0xff;415 assert(r.next_value == 0);
416 assert(r.random.uintAtMost(u8, 3) == 0);
417
418 // trigger the bias rejection code path
419 r.next_value = 0;
336 assert(r.random.uintAtMost(u8, 2) == 0);420 assert(r.random.uintAtMost(u8, 2) == 0);
337 assert(r.next_value == 1);421 // verify we incremented twice
422 assert(r.next_value == 2);
338423
339 r.next_value = 0xff;424 r.next_value = 0xff;
340 assert(r.random.intRangeAtMost(u8, 0, 0x7f) == 0x7f);425 assert(r.random.intRangeAtMost(u8, 0, 0x7f) == 0x7f);
...@@ -348,17 +433,43 @@ fn testRandomIntAtMost() void {...@@ -348,17 +433,43 @@ fn testRandomIntAtMost() void {
348 r.next_value = 0xff;433 r.next_value = 0xff;
349 assert(r.random.intRangeAtMost(i8, -0x80, -1) == -1);434 assert(r.random.intRangeAtMost(i8, -0x80, -1) == -1);
350435
351 r.next_value = 0xff;
352 assert(r.random.intRangeAtMost(i64, -0x8000000000000000, -1) == -1);
353 r.next_value = 0xff;436 r.next_value = 0xff;
354 assert(r.random.intRangeAtMost(i3, -4, -1) == -1);437 assert(r.random.intRangeAtMost(i3, -4, -1) == -1);
355 r.next_value = 0xff;438 r.next_value = 0xff;
356 assert(r.random.intRangeAtMost(i3, -2, 1) == 1);439 assert(r.random.intRangeAtMost(i3, -2, 1) == 1);
357440
358 // test retrying and eventually getting a good value441 assert(r.random.uintAtMost(u0, 0) == 0);
359 // start just out of bounds442}
360 r.next_value = 0x81;443
361 assert(r.random.uintAtMost(u8, 0x80) == 0);444test "Random Biased" {
445 var r = DefaultPrng.init(0);
446 // Not thoroughly checking the logic here.
447 // Just want to execute all the paths with different types.
448
449 assert(r.random.uintLessThanBiased(u1, 1) == 0);
450 assert(r.random.uintLessThanBiased(u32, 10) < 10);
451 assert(r.random.uintLessThanBiased(u64, 20) < 20);
452
453 assert(r.random.uintAtMostBiased(u0, 0) == 0);
454 assert(r.random.uintAtMostBiased(u1, 0) <= 0);
455 assert(r.random.uintAtMostBiased(u32, 10) <= 10);
456 assert(r.random.uintAtMostBiased(u64, 20) <= 20);
457
458 assert(r.random.intRangeLessThanBiased(u1, 0, 1) == 0);
459 assert(r.random.intRangeLessThanBiased(i1, -1, 0) == -1);
460 assert(r.random.intRangeLessThanBiased(u32, 10, 20) >= 10);
461 assert(r.random.intRangeLessThanBiased(i32, 10, 20) >= 10);
462 assert(r.random.intRangeLessThanBiased(u64, 20, 40) >= 20);
463 assert(r.random.intRangeLessThanBiased(i64, 20, 40) >= 20);
464
465 // uncomment for broken module error:
466 //assert(r.random.intRangeAtMostBiased(u0, 0, 0) == 0);
467 assert(r.random.intRangeAtMostBiased(u1, 0, 1) >= 0);
468 assert(r.random.intRangeAtMostBiased(i1, -1, 0) >= -1);
469 assert(r.random.intRangeAtMostBiased(u32, 10, 20) >= 10);
470 assert(r.random.intRangeAtMostBiased(i32, 10, 20) >= 10);
471 assert(r.random.intRangeAtMostBiased(u64, 20, 40) >= 20);
472 assert(r.random.intRangeAtMostBiased(i64, 20, 40) >= 20);
362}473}
363474
364// Generator to extend 64-bit seed values into longer sequences.475// Generator to extend 64-bit seed values into longer sequences.
...@@ -870,12 +981,16 @@ test "Random range" {...@@ -870,12 +981,16 @@ test "Random range" {
870}981}
871982
872fn testRange(r: *Random, start: i8, end: i8) void {983fn testRange(r: *Random, start: i8, end: i8) void {
984 testRangeBias(r, start, end, true);
985 testRangeBias(r, start, end, false);
986}
987fn testRangeBias(r: *Random, start: i8, end: i8, biased: bool) void {
873 const count = @intCast(usize, i32(end) - i32(start));988 const count = @intCast(usize, i32(end) - i32(start));
874 var values_buffer = []bool{false} ** 0x100;989 var values_buffer = []bool{false} ** 0x100;
875 const values = values_buffer[0..count];990 const values = values_buffer[0..count];
876 var i: usize = 0;991 var i: usize = 0;
877 while (i < count) {992 while (i < count) {
878 const value: i32 = r.intRangeLessThan(i8, start, end);993 const value: i32 = if (biased) r.intRangeLessThanBiased(i8, start, end) else r.intRangeLessThan(i8, start, end);
879 const index = @intCast(usize, value - start);994 const index = @intCast(usize, value - start);
880 if (!values[index]) {995 if (!values[index]) {
881 i += 1;996 i += 1;