authorgravatar for 102751849+Fri3dNstuff@users.noreply.github.comFri3dNstuff <102751849+Fri3dNstuff@users.noreply.github.com> 2024-09-24 03:15:57+03:00
committergravatar for noreply@github.comGitHub <noreply@github.com> 2024-09-23 17:15:57-07:00
logb2c53eb0d717989caf66dc8dfc0a7b1afb2607b3
tree2115b8d85aa75ffb3c73b570b66f0ebc2e72171f
parenta08f8d44dace9d7e7833bce092b47faa5bc0796e
signaturebadge-check Signed by PGP key B5690EEEBB952194

std.math: change gcd's implementation to use Stein's algorithm instead of Euclid's (#21077)

1 files changed, 35 insertions(+), 24 deletions(-)

lib/std/math/gcd.zig+35-24
......@@ -1,41 +1,50 @@
11//! Greatest common divisor (https://mathworld.wolfram.com/GreatestCommonDivisor.html)
22const std = @import("std");
3const expectEqual = std.testing.expectEqual;
43
5/// Returns the greatest common divisor (GCD) of two unsigned integers (a and b) which are not both zero.
6/// For example, the GCD of 8 and 12 is 4, that is, gcd(8, 12) == 4.
4/// Returns the greatest common divisor (GCD) of two unsigned integers (`a` and `b`) which are not both zero.
5/// For example, the GCD of `8` and `12` is `4`, that is, `gcd(8, 12) == 4`.
76pub fn gcd(a: anytype, b: anytype) @TypeOf(a, b) {
8
9 // only unsigned integers are allowed and not both must be zero
10 comptime switch (@typeInfo(@TypeOf(a, b))) {
11 .int => |int| std.debug.assert(int.signedness == .unsigned),
12 .comptime_int => {
13 std.debug.assert(a >= 0);
14 std.debug.assert(b >= 0);
15 },
16 else => unreachable,
7 const N = switch (@TypeOf(a, b)) {
8 // convert comptime_int to some sized int type for @ctz
9 comptime_int => std.math.IntFittingRange(@min(a, b), @max(a, b)),
10 else => |T| T,
1711 };
12 if (@typeInfo(N) != .int or @typeInfo(N).int.signedness != .unsigned) {
13 @compileError("`a` and `b` must be usigned integers");
14 }
15
16 // using an optimised form of Stein's algorithm:
17 // https://en.wikipedia.org/wiki/Binary_GCD_algorithm
1818 std.debug.assert(a != 0 or b != 0);
1919
20 // if one of them is zero, the other is returned
2120 if (a == 0) return b;
2221 if (b == 0) return a;
2322
24 // init vars
25 var x: @TypeOf(a, b) = a;
26 var y: @TypeOf(a, b) = b;
27 var m: @TypeOf(a, b) = a;
23 var x: N = a;
24 var y: N = b;
2825
29 // using the Euclidean algorithm (https://mathworld.wolfram.com/EuclideanAlgorithm.html)
30 while (y != 0) {
31 m = x % y;
32 x = y;
33 y = m;
26 const xz = @ctz(x);
27 const yz = @ctz(y);
28 const shift = @min(xz, yz);
29 x >>= @intCast(xz);
30 y >>= @intCast(yz);
31
32 var diff = y -% x;
33 while (diff != 0) : (diff = y -% x) {
34 // ctz is invariant under negation, we
35 // put it here to ease data dependencies,
36 // makes the CPU happy.
37 const zeros = @ctz(diff);
38 if (x > y) diff = -%diff;
39 y = @min(x, y);
40 x = diff >> @intCast(zeros);
3441 }
35 return x;
42 return y << @intCast(shift);
3643}
3744
38test "gcd" {
45test gcd {
46 const expectEqual = std.testing.expectEqual;
47
3948 try expectEqual(gcd(0, 5), 5);
4049 try expectEqual(gcd(5, 0), 5);
4150 try expectEqual(gcd(8, 12), 4);
......@@ -45,4 +54,6 @@ test "gcd" {
4554 try expectEqual(gcd(49865, 69811), 9973);
4655 try expectEqual(gcd(300_000, 2_300_000), 100_000);
4756 try expectEqual(gcd(90000000_000_000_000_000_000, 2), 2);
57 try expectEqual(gcd(@as(u80, 90000000_000_000_000_000_000), 2), 2);
58 try expectEqual(gcd(300_000, @as(u32, 2_300_000)), 100_000);
4859}