From 1ee59c5c31efca394d7e6d9da3f64c289a996b99 Mon Sep 17 00:00:00 2001 From: joachimschmidt557 Date: Mon, 6 Apr 2020 23:22:03 +0200 Subject: [PATCH] move big.rational.gcd to big.int.gcd --- lib/std/math/big/int.zig | 187 +++++++++++++++++++++++++++++++++ lib/std/math/big/rational.zig | 189 +--------------------------------- 2 files changed, 188 insertions(+), 188 deletions(-) diff --git a/lib/std/math/big/int.zig b/lib/std/math/big/int.zig index 51e97abf61..54f4a48d35 100644 --- a/lib/std/math/big/int.zig +++ b/lib/std/math/big/int.zig @@ -10,6 +10,7 @@ const minInt = std.math.minInt; pub const Limb = usize; pub const DoubleLimb = std.meta.IntType(false, 2 * Limb.bit_count); +pub const SignedDoubleLimb = std.meta.IntType(true, DoubleLimb.bit_count); pub const Log2Limb = math.Log2Int(Limb); comptime { @@ -1359,8 +1360,129 @@ pub const Int = struct { r[i] = a[i]; } } + + pub fn gcd(rma: *Int, x: Int, y: Int) !void { + rma.assertWritable(); + var r = rma; + var aliased = rma.limbs.ptr == x.limbs.ptr or rma.limbs.ptr == y.limbs.ptr; + + var sr: Int = undefined; + if (aliased) { + sr = try Int.initCapacity(rma.allocator.?, math.max(x.len(), y.len())); + r = &sr; + aliased = true; + } + defer if (aliased) { + rma.swap(r); + r.deinit(); + }; + + try gcdLehmer(r, x, y); + } + + fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void { + var x = try xa.clone(); + x.abs(); + defer x.deinit(); + + var y = try ya.clone(); + y.abs(); + defer y.deinit(); + + if (x.cmp(y) == .lt) { + x.swap(&y); + } + + var T = try Int.init(r.allocator.?); + defer T.deinit(); + + while (y.len() > 1) { + debug.assert(x.isPositive() and y.isPositive()); + debug.assert(x.len() >= y.len()); + + var xh: SignedDoubleLimb = x.limbs[x.len() - 1]; + var yh: SignedDoubleLimb = if (x.len() > y.len()) 0 else y.limbs[x.len() - 1]; + + var A: SignedDoubleLimb = 1; + var B: SignedDoubleLimb = 0; + var C: SignedDoubleLimb = 0; + var D: SignedDoubleLimb = 1; + + while (yh + C != 0 and yh + D != 0) { + const q = @divFloor(xh + A, yh + C); + const qp = @divFloor(xh + B, yh + D); + if (q != qp) { + break; + } + + var t = A - q * C; + A = C; + C = t; + t = B - q * D; + B = D; + D = t; + + t = xh - q * yh; + xh = yh; + yh = t; + } + + if (B == 0) { + // T = x % y, r is unused + try Int.divTrunc(r, &T, x, y); + debug.assert(T.isPositive()); + + x.swap(&y); + y.swap(&T); + } else { + var storage: [8]Limb = undefined; + const Ap = FixedIntFromSignedDoubleLimb(A, storage[0..2]); + const Bp = FixedIntFromSignedDoubleLimb(B, storage[2..4]); + const Cp = FixedIntFromSignedDoubleLimb(C, storage[4..6]); + const Dp = FixedIntFromSignedDoubleLimb(D, storage[6..8]); + + // T = Ax + By + try r.mul(x, Ap); + try T.mul(y, Bp); + try T.add(r.*, T); + + // u = Cx + Dy, r as u + try x.mul(x, Cp); + try r.mul(y, Dp); + try r.add(x, r.*); + + x.swap(&T); + y.swap(r); + } + } + + // euclidean algorithm + debug.assert(x.cmp(y) != .lt); + + while (!y.eqZero()) { + try Int.divTrunc(&T, r, x, y); + x.swap(&y); + y.swap(r); + } + + r.swap(&x); + } + }; +// Storage must live for the lifetime of the returned value +fn FixedIntFromSignedDoubleLimb(A: SignedDoubleLimb, storage: []Limb) Int { + std.debug.assert(storage.len >= 2); + + var A_is_positive = A >= 0; + const Au = @intCast(DoubleLimb, if (A < 0) -A else A); + storage[0] = @truncate(Limb, Au); + storage[1] = @truncate(Limb, Au >> Limb.bit_count); + var Ap = Int.initFixed(storage[0..2]); + Ap.setSign(A_is_positive); + return Ap; +} + // NOTE: All the following tests assume the max machine-word will be 64-bit. // // They will still run on larger than this and should pass, but the multi-limb code-paths @@ -2738,3 +2860,68 @@ test "big.int var args" { defer d.deinit(); testing.expect(a.cmp(d) != .gt); } + +test "big.int gcd non-one small" { + var a = try Int.initSet(testing.allocator, 17); + defer a.deinit(); + var b = try Int.initSet(testing.allocator, 97); + defer b.deinit(); + var r = try Int.init(testing.allocator); + defer r.deinit(); + + try r.gcd(a, b); + + testing.expect((try r.to(u32)) == 1); +} + +test "big.int gcd non-one small" { + var a = try Int.initSet(testing.allocator, 4864); + defer a.deinit(); + var b = try Int.initSet(testing.allocator, 3458); + defer b.deinit(); + var r = try Int.init(testing.allocator); + defer r.deinit(); + + try r.gcd(a, b); + + testing.expect((try r.to(u32)) == 38); +} + +test "big.int gcd non-one large" { + var a = try Int.initSet(testing.allocator, 0xffffffffffffffff); + defer a.deinit(); + var b = try Int.initSet(testing.allocator, 0xffffffffffffffff7777); + defer b.deinit(); + var r = try Int.init(testing.allocator); + defer r.deinit(); + + try r.gcd(a, b); + + testing.expect((try r.to(u32)) == 4369); +} + +test "big.int gcd large multi-limb result" { + var a = try Int.initSet(testing.allocator, 0x12345678123456781234567812345678123456781234567812345678); + defer a.deinit(); + var b = try Int.initSet(testing.allocator, 0x12345671234567123456712345671234567123456712345671234567); + defer b.deinit(); + var r = try Int.init(testing.allocator); + defer r.deinit(); + + try r.gcd(a, b); + + testing.expect((try r.to(u256)) == 0xf000000ff00000fff0000ffff000fffff00ffffff1); +} + +test "big.int gcd one large" { + var a = try Int.initSet(testing.allocator, 1897056385327307); + defer a.deinit(); + var b = try Int.initSet(testing.allocator, 2251799813685248); + defer b.deinit(); + var r = try Int.init(testing.allocator); + defer r.deinit(); + + try r.gcd(a, b); + + testing.expect((try r.to(u64)) == 1); +} diff --git a/lib/std/math/big/rational.zig b/lib/std/math/big/rational.zig index 8d2932a920..21e89ff8c4 100644 --- a/lib/std/math/big/rational.zig +++ b/lib/std/math/big/rational.zig @@ -447,7 +447,7 @@ pub const Rational = struct { const sign = r.p.isPositive(); r.p.abs(); - try gcd(&a, r.p, r.q); + try a.gcd(r.p, r.q); r.p.setSign(sign); const one = Int.initFixed(([_]Limb{1})[0..]); @@ -463,193 +463,6 @@ pub const Rational = struct { } }; -const SignedDoubleLimb = std.meta.IntType(true, DoubleLimb.bit_count); - -fn gcd(rma: *Int, x: Int, y: Int) !void { - rma.assertWritable(); - var r = rma; - var aliased = rma.limbs.ptr == x.limbs.ptr or rma.limbs.ptr == y.limbs.ptr; - - var sr: Int = undefined; - if (aliased) { - sr = try Int.initCapacity(rma.allocator.?, math.max(x.len(), y.len())); - r = &sr; - aliased = true; - } - defer if (aliased) { - rma.swap(r); - r.deinit(); - }; - - try gcdLehmer(r, x, y); -} - -// Storage must live for the lifetime of the returned value -fn FixedIntFromSignedDoubleLimb(A: SignedDoubleLimb, storage: []Limb) Int { - std.debug.assert(storage.len >= 2); - - var A_is_positive = A >= 0; - const Au = @intCast(DoubleLimb, if (A < 0) -A else A); - storage[0] = @truncate(Limb, Au); - storage[1] = @truncate(Limb, Au >> Limb.bit_count); - var Ap = Int.initFixed(storage[0..2]); - Ap.setSign(A_is_positive); - return Ap; -} - -fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void { - var x = try xa.clone(); - x.abs(); - defer x.deinit(); - - var y = try ya.clone(); - y.abs(); - defer y.deinit(); - - if (x.cmp(y) == .lt) { - x.swap(&y); - } - - var T = try Int.init(r.allocator.?); - defer T.deinit(); - - while (y.len() > 1) { - debug.assert(x.isPositive() and y.isPositive()); - debug.assert(x.len() >= y.len()); - - var xh: SignedDoubleLimb = x.limbs[x.len() - 1]; - var yh: SignedDoubleLimb = if (x.len() > y.len()) 0 else y.limbs[x.len() - 1]; - - var A: SignedDoubleLimb = 1; - var B: SignedDoubleLimb = 0; - var C: SignedDoubleLimb = 0; - var D: SignedDoubleLimb = 1; - - while (yh + C != 0 and yh + D != 0) { - const q = @divFloor(xh + A, yh + C); - const qp = @divFloor(xh + B, yh + D); - if (q != qp) { - break; - } - - var t = A - q * C; - A = C; - C = t; - t = B - q * D; - B = D; - D = t; - - t = xh - q * yh; - xh = yh; - yh = t; - } - - if (B == 0) { - // T = x % y, r is unused - try Int.divTrunc(r, &T, x, y); - debug.assert(T.isPositive()); - - x.swap(&y); - y.swap(&T); - } else { - var storage: [8]Limb = undefined; - const Ap = FixedIntFromSignedDoubleLimb(A, storage[0..2]); - const Bp = FixedIntFromSignedDoubleLimb(B, storage[2..4]); - const Cp = FixedIntFromSignedDoubleLimb(C, storage[4..6]); - const Dp = FixedIntFromSignedDoubleLimb(D, storage[6..8]); - - // T = Ax + By - try r.mul(x, Ap); - try T.mul(y, Bp); - try T.add(r.*, T); - - // u = Cx + Dy, r as u - try x.mul(x, Cp); - try r.mul(y, Dp); - try r.add(x, r.*); - - x.swap(&T); - y.swap(r); - } - } - - // euclidean algorithm - debug.assert(x.cmp(y) != .lt); - - while (!y.eqZero()) { - try Int.divTrunc(&T, r, x, y); - x.swap(&y); - y.swap(r); - } - - r.swap(&x); -} - -test "big.rational gcd non-one small" { - var a = try Int.initSet(testing.allocator, 17); - defer a.deinit(); - var b = try Int.initSet(testing.allocator, 97); - defer b.deinit(); - var r = try Int.init(testing.allocator); - defer r.deinit(); - - try gcd(&r, a, b); - - testing.expect((try r.to(u32)) == 1); -} - -test "big.rational gcd non-one small" { - var a = try Int.initSet(testing.allocator, 4864); - defer a.deinit(); - var b = try Int.initSet(testing.allocator, 3458); - defer b.deinit(); - var r = try Int.init(testing.allocator); - defer r.deinit(); - - try gcd(&r, a, b); - - testing.expect((try r.to(u32)) == 38); -} - -test "big.rational gcd non-one large" { - var a = try Int.initSet(testing.allocator, 0xffffffffffffffff); - defer a.deinit(); - var b = try Int.initSet(testing.allocator, 0xffffffffffffffff7777); - defer b.deinit(); - var r = try Int.init(testing.allocator); - defer r.deinit(); - - try gcd(&r, a, b); - - testing.expect((try r.to(u32)) == 4369); -} - -test "big.rational gcd large multi-limb result" { - var a = try Int.initSet(testing.allocator, 0x12345678123456781234567812345678123456781234567812345678); - defer a.deinit(); - var b = try Int.initSet(testing.allocator, 0x12345671234567123456712345671234567123456712345671234567); - defer b.deinit(); - var r = try Int.init(testing.allocator); - defer r.deinit(); - - try gcd(&r, a, b); - - testing.expect((try r.to(u256)) == 0xf000000ff00000fff0000ffff000fffff00ffffff1); -} - -test "big.rational gcd one large" { - var a = try Int.initSet(testing.allocator, 1897056385327307); - defer a.deinit(); - var b = try Int.initSet(testing.allocator, 2251799813685248); - defer b.deinit(); - var r = try Int.init(testing.allocator); - defer r.deinit(); - - try gcd(&r, a, b); - - testing.expect((try r.to(u64)) == 1); -} - fn extractLowBits(a: Int, comptime T: type) T { testing.expect(@typeInfo(T) == .Int);