From c577d1e61e7243a3f3fdc8f34b7b38fcad32db34 Mon Sep 17 00:00:00 2001 From: Hila Friedman Date: Wed, 14 Aug 2024 20:45:10 +0300 Subject: [PATCH 1/3] Change gcd implementation to Stein's algorithm --- lib/std/math/gcd.zig | 65 ++++++++++++++++++++++++++++---------------- 1 file changed, 41 insertions(+), 24 deletions(-) diff --git a/lib/std/math/gcd.zig b/lib/std/math/gcd.zig index 298bf5fc2b70..0c732a5c9746 100644 --- a/lib/std/math/gcd.zig +++ b/lib/std/math/gcd.zig @@ -1,41 +1,53 @@ //! Greatest common divisor (https://mathworld.wolfram.com/GreatestCommonDivisor.html) const std = @import("std"); -const expectEqual = std.testing.expectEqual; -/// Returns the greatest common divisor (GCD) of two unsigned integers (a and b) which are not both zero. -/// For example, the GCD of 8 and 12 is 4, that is, gcd(8, 12) == 4. +/// Returns the greatest common divisor (GCD) of two unsigned integers (`a` and `b`) which are not both zero. +/// For example, the GCD of `8` and `12` is `4`, that is, `gcd(8, 12) == 4`. pub fn gcd(a: anytype, b: anytype) @TypeOf(a, b) { - - // only unsigned integers are allowed and not both must be zero - comptime switch (@typeInfo(@TypeOf(a, b))) { - .Int => |int| std.debug.assert(int.signedness == .unsigned), - .ComptimeInt => { - std.debug.assert(a >= 0); - std.debug.assert(b >= 0); - }, - else => unreachable, + const N = switch (@TypeOf(a, b)) { + // convert comptime_int to some sized int so we can @ctz on it. + // type coercion takes care of the conversion back to comptime_int + // at function's return + comptime_int => std.math.IntFittingRange(@min(a, b), @max(a, b)), + else => |T| T, }; + // integers are unsigned, at least one is nonzero + comptime std.debug.assert(@typeInfo(N).Int.signedness == .unsigned); std.debug.assert(a != 0 or b != 0); - // if one of them is zero, the other is returned + // using Stein's algorithm (https://en.wikipedia.org/wiki/Binary_GCD_algorithm) if (a == 0) return b; if (b == 0) return a; - // init vars - var x: @TypeOf(a, b) = a; - var y: @TypeOf(a, b) = b; - var m: @TypeOf(a, b) = a; + var x: N = a; + var y: N = b; + + const i = @ctz(x); + const j = @ctz(y); + // x, y are nonzero, @intCast(@ctz(self)) does not overflow + x >>= @intCast(i); + y >>= @intCast(j); + + // invariants: x, y are odd + while (true) { + // ensure x ≥ y + if (y > x) std.mem.swap(N, &x, &y); + + // odd - odd is even + x -= y; - // using the Euclidean algorithm (https://mathworld.wolfram.com/EuclideanAlgorithm.html) - while (y != 0) { - m = x % y; - x = y; - y = m; + // gcd(0, y) == y, remultiply by the common power of 2 + if (x == 0) return y << @intCast(@min(i, j)); + + // x is nonzero, @intCast(@ctz(self)) does not overflow + // x is even, its value decreases + x >>= @intCast(@ctz(x)); } - return x; } -test "gcd" { +test gcd { + const expectEqual = std.testing.expectEqual; + try expectEqual(gcd(0, 5), 5); try expectEqual(gcd(5, 0), 5); try expectEqual(gcd(8, 12), 4); @@ -45,4 +57,9 @@ test "gcd" { try expectEqual(gcd(49865, 69811), 9973); try expectEqual(gcd(300_000, 2_300_000), 100_000); try expectEqual(gcd(90000000_000_000_000_000_000, 2), 2); + try expectEqual(gcd(@as(u80, 90000000_000_000_000_000_000), 2), 2); + + // an important test case! @ctz's returned type is different from the + // type of >>'s rhs - an @intCast is required! + try expectEqual(gcd(300_000, @as(u32, 2_300_000)), 100_000); } From fc4b98fea48ba418ecb6c173d9f5f46c615b281d Mon Sep 17 00:00:00 2001 From: Hila Friedman Date: Wed, 28 Aug 2024 10:30:32 +0300 Subject: [PATCH 2/3] Improve implementation performance --- lib/std/math/gcd.zig | 48 ++++++++++++++++++-------------------------- 1 file changed, 20 insertions(+), 28 deletions(-) diff --git a/lib/std/math/gcd.zig b/lib/std/math/gcd.zig index 0c732a5c9746..f35d81125f75 100644 --- a/lib/std/math/gcd.zig +++ b/lib/std/math/gcd.zig @@ -5,44 +5,39 @@ const std = @import("std"); /// For example, the GCD of `8` and `12` is `4`, that is, `gcd(8, 12) == 4`. pub fn gcd(a: anytype, b: anytype) @TypeOf(a, b) { const N = switch (@TypeOf(a, b)) { - // convert comptime_int to some sized int so we can @ctz on it. - // type coercion takes care of the conversion back to comptime_int - // at function's return + // convert comptime_int to some sized int type for @ctz comptime_int => std.math.IntFittingRange(@min(a, b), @max(a, b)), else => |T| T, }; - // integers are unsigned, at least one is nonzero comptime std.debug.assert(@typeInfo(N).Int.signedness == .unsigned); + + // using an optimised form of Stein's algorithm: + // https://en.wikipedia.org/wiki/Binary_GCD_algorithm std.debug.assert(a != 0 or b != 0); - // using Stein's algorithm (https://en.wikipedia.org/wiki/Binary_GCD_algorithm) if (a == 0) return b; if (b == 0) return a; var x: N = a; var y: N = b; - const i = @ctz(x); - const j = @ctz(y); - // x, y are nonzero, @intCast(@ctz(self)) does not overflow - x >>= @intCast(i); - y >>= @intCast(j); - - // invariants: x, y are odd - while (true) { - // ensure x ≥ y - if (y > x) std.mem.swap(N, &x, &y); - - // odd - odd is even - x -= y; - - // gcd(0, y) == y, remultiply by the common power of 2 - if (x == 0) return y << @intCast(@min(i, j)); - - // x is nonzero, @intCast(@ctz(self)) does not overflow - // x is even, its value decreases - x >>= @intCast(@ctz(x)); + const xz = @ctz(x); + const yz = @ctz(y); + const shift = @min(xz, yz); + x >>= @intCast(xz); + y >>= @intCast(yz); + + var diff = y -% x; + while (diff != 0) : (diff = y -% x) { + // ctz is invariant under negation, we + // put it here to ease data dependencies, + // makes the CPU happy. + const zeros = @ctz(diff); + if (x > y) diff = -%diff; + y = @min(x, y); + x = diff >> @intCast(zeros); } + return y << @intCast(shift); } test gcd { @@ -58,8 +53,5 @@ test gcd { try expectEqual(gcd(300_000, 2_300_000), 100_000); try expectEqual(gcd(90000000_000_000_000_000_000, 2), 2); try expectEqual(gcd(@as(u80, 90000000_000_000_000_000_000), 2), 2); - - // an important test case! @ctz's returned type is different from the - // type of >>'s rhs - an @intCast is required! try expectEqual(gcd(300_000, @as(u32, 2_300_000)), 100_000); } From 7d0b2b433fcd29f9311035eba2449d8b2f312d09 Mon Sep 17 00:00:00 2001 From: Hila Friedman Date: Fri, 20 Sep 2024 10:02:23 +0300 Subject: [PATCH 3/3] Add compile error for bad types --- lib/std/math/gcd.zig | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lib/std/math/gcd.zig b/lib/std/math/gcd.zig index f35d81125f75..16ca7846f19a 100644 --- a/lib/std/math/gcd.zig +++ b/lib/std/math/gcd.zig @@ -9,7 +9,9 @@ pub fn gcd(a: anytype, b: anytype) @TypeOf(a, b) { comptime_int => std.math.IntFittingRange(@min(a, b), @max(a, b)), else => |T| T, }; - comptime std.debug.assert(@typeInfo(N).Int.signedness == .unsigned); + if (@typeInfo(N) != .int or @typeInfo(N).int.signedness != .unsigned) { + @compileError("`a` and `b` must be usigned integers"); + } // using an optimised form of Stein's algorithm: // https://en.wikipedia.org/wiki/Binary_GCD_algorithm