1//! Greatest common divisor (https://mathworld.wolfram.com/GreatestCommonDivisor.html)
 2const std = @import("std");
 3
 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`.
 6pub fn gcd(a: anytype, b: anytype) @TypeOf(a, b) {
 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,
11    };
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
18    std.debug.assert(a != 0 or b != 0);
19
20    if (a == 0) return b;
21    if (b == 0) return a;
22
23    var x: N = a;
24    var y: N = b;
25
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);
41    }
42    return y << @intCast(shift);
43}
44
45test gcd {
46    const expectEqual = std.testing.expectEqual;
47
48    try expectEqual(gcd(0, 5), 5);
49    try expectEqual(gcd(5, 0), 5);
50    try expectEqual(gcd(8, 12), 4);
51    try expectEqual(gcd(12, 8), 4);
52    try expectEqual(gcd(33, 77), 11);
53    try expectEqual(gcd(77, 33), 11);
54    try expectEqual(gcd(49865, 69811), 9973);
55    try expectEqual(gcd(300_000, 2_300_000), 100_000);
56    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);
59}