authorgravatar for marc@tiehu.isMarc Tiehuis <marc@tiehu.is> 2019-04-03 17:20:23+13:00
committergravatar for marc@tiehu.isMarc Tiehuis <marc@tiehu.is> 2019-04-11 19:36:35+12:00
log87d8ecda462688c597c726b1da5dbd5f8478e0fc
treecdce715471f9868cda4bf0c19b2645fb81471284
parent30788a98b1bb58886d409464baf7eb01e428d939

Fix math.big.Int divN/gcdLehmer and fuzz-test failures

2 files changed, 96 insertions(+), 88 deletions(-)

std/math/big/int.zig+77-54
......@@ -67,7 +67,7 @@ pub const Int = struct {
6767 .len = limbs.len,
6868 };
6969
70 self.normN(limbs.len);
70 self.normalize(limbs.len);
7171 return self;
7272 }
7373
......@@ -508,28 +508,12 @@ pub const Int = struct {
508508 return cmp(a, b) == 0;
509509 }
510510
511 // Normalize for a possible single carry digit.
512 //
513 // [1, 2, 3, 4, 0] -> [1, 2, 3, 4]
514 // [1, 2, 3, 4, 5] -> [1, 2, 3, 4, 5]
515 // [0] -> [0]
516 fn norm1(r: *Int, length: usize) void {
517 debug.assert(length > 0);
518 debug.assert(length <= r.limbs.len);
519
520 if (r.limbs[length - 1] == 0) {
521 r.len = if (length > 1) length - 1 else 1;
522 } else {
523 r.len = length;
524 }
525 }
526
527511 // Normalize a possible sequence of leading zeros.
528512 //
529513 // [1, 2, 3, 4, 0] -> [1, 2, 3, 4]
530514 // [1, 2, 0, 0, 0] -> [1, 2]
531515 // [0, 0, 0, 0, 0] -> [0]
532 fn normN(r: *Int, length: usize) void {
516 fn normalize(r: *Int, length: usize) void {
533517 debug.assert(length > 0);
534518 debug.assert(length <= r.limbs.len);
535519
......@@ -577,11 +561,11 @@ pub const Int = struct {
577561 if (a.len >= b.len) {
578562 try r.ensureCapacity(a.len + 1);
579563 lladd(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);
580 r.norm1(a.len + 1);
564 r.normalize(a.len + 1);
581565 } else {
582566 try r.ensureCapacity(b.len + 1);
583567 lladd(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);
584 r.norm1(b.len + 1);
568 r.normalize(b.len + 1);
585569 }
586570
587571 r.positive = a.positive;
......@@ -630,12 +614,12 @@ pub const Int = struct {
630614 if (a.cmp(b) >= 0) {
631615 try r.ensureCapacity(a.len + 1);
632616 llsub(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);
633 r.normN(a.len);
617 r.normalize(a.len);
634618 r.positive = true;
635619 } else {
636620 try r.ensureCapacity(b.len + 1);
637621 llsub(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);
638 r.normN(b.len);
622 r.normalize(b.len);
639623 r.positive = false;
640624 }
641625 } else {
......@@ -643,12 +627,12 @@ pub const Int = struct {
643627 if (a.cmp(b) < 0) {
644628 try r.ensureCapacity(a.len + 1);
645629 llsub(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);
646 r.normN(a.len);
630 r.normalize(a.len);
647631 r.positive = false;
648632 } else {
649633 try r.ensureCapacity(b.len + 1);
650634 llsub(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);
651 r.normN(b.len);
635 r.normalize(b.len);
652636 r.positive = true;
653637 }
654638 }
......@@ -708,7 +692,7 @@ pub const Int = struct {
708692 }
709693
710694 r.positive = a.positive == b.positive;
711 r.normN(a.len + b.len);
695 r.normalize(a.len + b.len);
712696 }
713697
714698 // a + b * c + *carry, sets carry to the overflow bits
......@@ -817,7 +801,7 @@ pub const Int = struct {
817801 try quo.ensureCapacity(a.len);
818802
819803 lldiv1(quo.limbs[0..], &rem.limbs[0], a.limbs[ab_zero_limb_count..a.len], b.limbs[b.len - 1]);
820 quo.norm1(a.len - ab_zero_limb_count);
804 quo.normalize(a.len - ab_zero_limb_count);
821805 quo.positive = a.positive == b.positive;
822806
823807 rem.len = 1;
......@@ -846,13 +830,10 @@ pub const Int = struct {
846830 try divN(quo.allocator.?, quo, rem, &x, &y);
847831
848832 quo.positive = a.positive == b.positive;
833 }
849834
850 // If dividend had trailing zeros beyond divisor, add extra trailing limbs.
851 // Single-limb division never has multi-limb remainder so nothing to add.
852 if (a_zero_limb_count > b_zero_limb_count) {
853 const shift = a_zero_limb_count - b_zero_limb_count;
854 try rem.shiftLeft(rem.*, shift * Limb.bit_count);
855 }
835 if (ab_zero_limb_count != 0) {
836 try rem.shiftLeft(rem.*, ab_zero_limb_count * Limb.bit_count);
856837 }
857838 }
858839
......@@ -933,7 +914,7 @@ pub const Int = struct {
933914 tmp.limbs[0] = if (i >= 2) x.limbs[i - 2] else 0;
934915 tmp.limbs[1] = if (i >= 1) x.limbs[i - 1] else 0;
935916 tmp.limbs[2] = x.limbs[i];
936 tmp.normN(3);
917 tmp.normalize(3);
937918
938919 while (true) {
939920 // 2x1 limb multiplication unrolled against single-limb q[i-t-1]
......@@ -941,7 +922,7 @@ pub const Int = struct {
941922 r.limbs[0] = addMulLimbWithCarry(0, if (t >= 1) y.limbs[t - 1] else 0, q.limbs[i - t - 1], &carry);
942923 r.limbs[1] = addMulLimbWithCarry(0, y.limbs[t], q.limbs[i - t - 1], &carry);
943924 r.limbs[2] = carry;
944 r.normN(3);
925 r.normalize(3);
945926
946927 if (r.cmpAbs(tmp) <= 0) {
947928 break;
......@@ -964,10 +945,10 @@ pub const Int = struct {
964945 }
965946
966947 // Denormalize
967 q.normN(q.len);
948 q.normalize(q.len);
968949
969950 try r.shiftRight(x.*, norm_shift);
970 r.normN(r.len);
951 r.normalize(r.len);
971952 }
972953
973954 // r = a << shift, in other words, r = a * 2^shift
......@@ -976,7 +957,7 @@ pub const Int = struct {
976957
977958 try r.ensureCapacity(a.len + (shift / Limb.bit_count) + 1);
978959 llshl(r.limbs[0..], a.limbs[0..a.len], shift);
979 r.norm1(a.len + (shift / Limb.bit_count) + 1);
960 r.normalize(a.len + (shift / Limb.bit_count) + 1);
980961 r.positive = a.positive;
981962 }
982963
......@@ -1076,11 +1057,11 @@ pub const Int = struct {
10761057 if (a.len > b.len) {
10771058 try r.ensureCapacity(b.len);
10781059 lland(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);
1079 r.normN(b.len);
1060 r.normalize(b.len);
10801061 } else {
10811062 try r.ensureCapacity(a.len);
10821063 lland(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);
1083 r.normN(a.len);
1064 r.normalize(a.len);
10841065 }
10851066 }
10861067
......@@ -1102,11 +1083,11 @@ pub const Int = struct {
11021083 if (a.len > b.len) {
11031084 try r.ensureCapacity(a.len);
11041085 llxor(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);
1105 r.normN(a.len);
1086 r.normalize(a.len);
11061087 } else {
11071088 try r.ensureCapacity(b.len);
11081089 llxor(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);
1109 r.normN(b.len);
1090 r.normalize(b.len);
11101091 }
11111092 }
11121093
......@@ -1130,7 +1111,9 @@ pub const Int = struct {
11301111// They will still run on larger than this and should pass, but the multi-limb code-paths
11311112// may be untested in some cases.
11321113
1133const al = debug.global_allocator;
1114var buffer: [64 * 8192]u8 = undefined;
1115var fixed = std.heap.FixedBufferAllocator.init(buffer[0..]);
1116const al = &fixed.allocator;
11341117
11351118test "big.int comptime_int set" {
11361119 comptime var s = 0xefffffff00000001eeeeeeefaaaaaaab;
......@@ -1179,7 +1162,7 @@ test "big.int to target too small error" {
11791162 testing.expectError(error.TargetTooSmall, a.to(u8));
11801163}
11811164
1182test "big.int norm1" {
1165test "big.int normalize" {
11831166 var a = try Int.init(al);
11841167 try a.ensureCapacity(8);
11851168
......@@ -1187,26 +1170,26 @@ test "big.int norm1" {
11871170 a.limbs[1] = 2;
11881171 a.limbs[2] = 3;
11891172 a.limbs[3] = 0;
1190 a.norm1(4);
1173 a.normalize(4);
11911174 testing.expect(a.len == 3);
11921175
11931176 a.limbs[0] = 1;
11941177 a.limbs[1] = 2;
11951178 a.limbs[2] = 3;
1196 a.norm1(3);
1179 a.normalize(3);
11971180 testing.expect(a.len == 3);
11981181
11991182 a.limbs[0] = 0;
12001183 a.limbs[1] = 0;
1201 a.norm1(2);
1184 a.normalize(2);
12021185 testing.expect(a.len == 1);
12031186
12041187 a.limbs[0] = 0;
1205 a.norm1(1);
1188 a.normalize(1);
12061189 testing.expect(a.len == 1);
12071190}
12081191
1209test "big.int normN" {
1192test "big.int normalize multi" {
12101193 var a = try Int.init(al);
12111194 try a.ensureCapacity(8);
12121195
......@@ -1214,24 +1197,24 @@ test "big.int normN" {
12141197 a.limbs[1] = 2;
12151198 a.limbs[2] = 0;
12161199 a.limbs[3] = 0;
1217 a.normN(4);
1200 a.normalize(4);
12181201 testing.expect(a.len == 2);
12191202
12201203 a.limbs[0] = 1;
12211204 a.limbs[1] = 2;
12221205 a.limbs[2] = 3;
1223 a.normN(3);
1206 a.normalize(3);
12241207 testing.expect(a.len == 3);
12251208
12261209 a.limbs[0] = 0;
12271210 a.limbs[1] = 0;
12281211 a.limbs[2] = 0;
12291212 a.limbs[3] = 0;
1230 a.normN(4);
1213 a.normalize(4);
12311214 testing.expect(a.len == 1);
12321215
12331216 a.limbs[0] = 0;
1234 a.normN(1);
1217 a.normalize(1);
12351218 testing.expect(a.len == 1);
12361219}
12371220
......@@ -2065,7 +2048,9 @@ test "big.int div multi-multi zero-limb trailing (with rem)" {
20652048 try Int.divTrunc(&q, &r, a, b);
20662049
20672050 testing.expect((try q.to(u128)) == 0x10000000000000000);
2068 testing.expect((try r.to(u128)) == 0x44444443444444431111111111111111);
2051
2052 const rs = try r.toString(al, 16);
2053 testing.expect(std.mem.eql(u8, rs, "4444444344444443111111111111111100000000000000000000000000000000"));
20692054}
20702055
20712056test "big.int div multi-multi zero-limb trailing (with rem) and dividend zero-limb count > divisor zero-limb count" {
......@@ -2077,7 +2062,9 @@ test "big.int div multi-multi zero-limb trailing (with rem) and dividend zero-li
20772062 try Int.divTrunc(&q, &r, a, b);
20782063
20792064 testing.expect((try q.to(u128)) == 0x1);
2080 testing.expect((try r.to(u128)) == 0x44444443444444431111111111111111);
2065
2066 const rs = try r.toString(al, 16);
2067 testing.expect(std.mem.eql(u8, rs, "444444434444444311111111111111110000000000000000"));
20812068}
20822069
20832070test "big.int div multi-multi zero-limb trailing (with rem) and dividend zero-limb count < divisor zero-limb count" {
......@@ -2095,6 +2082,42 @@ test "big.int div multi-multi zero-limb trailing (with rem) and dividend zero-li
20952082 testing.expect(std.mem.eql(u8, rs, "4e11f2baa5896a321d463b543d0104e30000000000000000"));
20962083}
20972084
2085test "big.int div multi-multi fuzz case #1" {
2086 var a = try Int.init(al);
2087 var b = try Int.init(al);
2088
2089 try a.setString(16, "ffffffffffffffffffffffffffffc00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff80000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffc00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");
2090 try b.setString(16, "3ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe0000000000000000000000000000000000001ffffffffffffffffffffffffffffffffffffffffffffffffffc000000000000000000000000000000007fffffffffff");
2091
2092 var q = try Int.init(al);
2093 var r = try Int.init(al);
2094 try Int.divTrunc(&q, &r, a, b);
2095
2096 const qs = try q.toString(al, 16);
2097 testing.expect(std.mem.eql(u8, qs, "3ffffffffffffffffffffffffffff0000000000000000000000000000000000001ffffffffffffffffffffffffffff7fffffffe000000000000000000000000000180000000000000000000003fffffbfffffffdfffffffffffffeffff800000100101000000100000000020003fffffdfbfffffe3ffffffffffffeffff7fffc00800a100000017ffe000002000400007efbfff7fe9f00000037ffff3fff7fffa004006100000009ffe00000190038200bf7d2ff7fefe80400060000f7d7f8fbf9401fe38e0403ffc0bdffffa51102c300d7be5ef9df4e5060007b0127ad3fa69f97d0f820b6605ff617ddf7f32ad7a05c0d03f2e7bc78a6000e087a8bbcdc59e07a5a079128a7861f553ddebed7e8e56701756f9ead39b48cd1b0831889ea6ec1fddf643d0565b075ff07e6caea4e2854ec9227fd635ed60a2f5eef2893052ffd54718fa08604acbf6a15e78a467c4a3c53c0278af06c4416573f925491b195e8fd79302cb1aaf7caf4ecfc9aec1254cc969786363ac729f914c6ddcc26738d6b0facd54eba026580aba2eb6482a088b0d224a8852420b91ec1"));
2098
2099 const rs = try r.toString(al, 16);
2100 testing.expect(std.mem.eql(u8, rs, "310d1d4c414426b4836c2635bad1df3a424e50cbdd167ffccb4dfff57d36b4aae0d6ca0910698220171a0f3373c1060a046c2812f0027e321f72979daa5e7973214170d49e885de0c0ecc167837d44502430674a82522e5df6a0759548052420b91ec1"));
2101}
2102
2103test "big.int div multi-multi fuzz case #2" {
2104 var a = try Int.init(al);
2105 var b = try Int.init(al);
2106
2107 try a.setString(16, "3ffffffffe00000000000000000000000000fffffffffffffffffffffffffffffffffffffffffffffffffffffffffe000000000000000000000000000000000000000000000000000000000000001fffffffffffffffff800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000ffffffffffffffffffffc000000000000000000000000000000000000000000000000000000000000000");
2108 try b.setString(16, "ffc0000000000000000000000000000000000000000000000000");
2109
2110 var q = try Int.init(al);
2111 var r = try Int.init(al);
2112 try Int.divTrunc(&q, &r, a, b);
2113
2114 const qs = try q.toString(al, 16);
2115 testing.expect(std.mem.eql(u8, qs, "40100400fe3f8fe3f8fe3f8fe3f8fe3f8fe4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f91e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4992649926499264991e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4792e4b92e4b92e4b92e4b92a4a92a4a92a4"));
2116
2117 const rs = try r.toString(al, 16);
2118 testing.expect(std.mem.eql(u8, rs, "a900000000000000000000000000000000000000000000000000"));
2119}
2120
20982121test "big.int shift-right single" {
20992122 var a = try Int.initSet(al, 0xffff0000);
21002123 try a.shiftRight(a, 16);
std/math/big/rational.zig+19-34
......@@ -222,6 +222,7 @@ pub const Rational = struct {
222222 exp += 1;
223223 }
224224 if (mantissa >> msize1 != 1) {
225 // NOTE: This can be hit if the limb size is small (u8/16).
225226 @panic("unexpected bits in result");
226227 }
227228
......@@ -421,8 +422,6 @@ pub const Rational = struct {
421422 }
422423};
423424
424var al = debug.global_allocator;
425
426425const SignedDoubleLimb = @IntType(true, DoubleLimb.bit_count);
427426
428427fn gcd(rma: *Int, x: Int, y: Int) !void {
......@@ -441,11 +440,7 @@ fn gcd(rma: *Int, x: Int, y: Int) !void {
441440 r.deinit();
442441 };
443442
444 if (x.cmp(y) > 0) {
445 try gcdLehmer(r, x, y);
446 } else {
447 try gcdLehmer(r, y, x);
448 }
443 try gcdLehmer(r, x, y);
449444}
450445
451446// Storage must live for the lifetime of the returned value
......@@ -461,9 +456,6 @@ fn FixedIntFromSignedDoubleLimb(A: SignedDoubleLimb, storage: []Limb) Int {
461456 return Ap;
462457}
463458
464// Handbook of Applied Cryptography, 14.57
465//
466// r = gcd(x, y) where x, y > 0
467459fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void {
468460 var x = try xa.clone();
469461 x.abs();
......@@ -484,18 +476,8 @@ fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void {
484476 debug.assert(x.positive and y.positive);
485477 debug.assert(x.len >= y.len);
486478
487 // chop the leading zeros of the limbs and normalize
488 const offset = @clz(x.limbs[x.len - 1]);
489
490 var xh: SignedDoubleLimb = math.shl(Limb, x.limbs[x.len - 1], offset) |
491 math.shr(Limb, x.limbs[x.len - 2], Limb.bit_count - offset);
492
493 var yh: SignedDoubleLimb = if (y.len == x.len)
494 math.shl(Limb, y.limbs[y.len - 1], offset) | math.shr(Limb, y.limbs[y.len - 2], Limb.bit_count - offset)
495 else if (y.len == x.len - 1)
496 math.shr(Limb, y.limbs[y.len - 2], Limb.bit_count - offset)
497 else
498 0;
479 var xh: SignedDoubleLimb = x.limbs[x.len - 1];
480 var yh: SignedDoubleLimb = if (x.len > y.len) 0 else y.limbs[x.len - 1];
499481
500482 var A: SignedDoubleLimb = 1;
501483 var B: SignedDoubleLimb = 0;
......@@ -546,13 +528,7 @@ fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void {
546528 try r.add(x, r.*);
547529
548530 x.swap(&T);
549 x.abs();
550531 y.swap(r);
551 y.abs();
552
553 if (x.cmp(y) < 0) {
554 x.swap(&y);
555 }
556532 }
557533 }
558534
......@@ -568,6 +544,10 @@ fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void {
568544 r.swap(&x);
569545}
570546
547var buffer: [64 * 8192]u8 = undefined;
548var fixed = std.heap.FixedBufferAllocator.init(buffer[0..]);
549var al = &fixed.allocator;
550
571551test "big.rational gcd non-one small" {
572552 var a = try Int.initSet(al, 17);
573553 var b = try Int.initSet(al, 97);
......@@ -608,6 +588,16 @@ test "big.rational gcd large multi-limb result" {
608588 testing.expect((try r.to(u256)) == 0xf000000ff00000fff0000ffff000fffff00ffffff1);
609589}
610590
591test "big.rational gcd one large" {
592 var a = try Int.initSet(al, 1897056385327307);
593 var b = try Int.initSet(al, 2251799813685248);
594 var r = try Int.init(al);
595
596 try gcd(&r, a, b);
597
598 testing.expect((try r.to(u64)) == 1);
599}
600
611601fn extractLowBits(a: Int, comptime T: type) T {
612602 testing.expect(@typeId(T) == builtin.TypeId.Int);
613603
......@@ -721,12 +711,7 @@ test "big.rational toFloat" {
721711}
722712
723713test "big.rational set/to Float round-trip" {
724 // toFloat allocates memory in a loop so we need to free it
725 var buf: [512 * 1024]u8 = undefined;
726 var fixed = std.heap.FixedBufferAllocator.init(buf[0..]);
727
728 var a = try Rational.init(&fixed.allocator);
729
714 var a = try Rational.init(al);
730715 var prng = std.rand.DefaultPrng.init(0x5EED);
731716 var i: usize = 0;
732717 while (i < 512) : (i += 1) {