bench_gcd.zig

const std = @import("std");

// Old implementation (current master)
fn gcd_old(a: anytype, b: anytype) @TypeOf(a, b) {
    const N = switch (@TypeOf(a, b)) {
        comptime_int => std.math.IntFittingRange(@min(a, b), @max(a, b)),
        else => |T| T,
    };
    if (@typeInfo(N) != .int or @typeInfo(N).int.signedness != .unsigned) {
        @compileError("`a` and `b` must be unsigned integers");
    }

    std.debug.assert(a != 0 or b != 0);

    if (a == 0) return b;
    if (b == 0) return a;

    var x: N = a;
    var y: N = b;

    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) {
        const zeros = @ctz(diff);
        if (x > y) diff = -%diff;
        y = @min(x, y);
        x = diff >> @intCast(zeros);
    }
    return y << @intCast(shift);
}

// New implementation (PR #25533)
fn gcd_new(a: anytype, b: anytype) @TypeOf(a, b) {
    const N = switch (@TypeOf(a, b)) {
        comptime_int => std.math.IntFittingRange(@min(a, b), @max(a, b)),
        else => |T| T,
    };

    if (@typeInfo(N) != .int or @typeInfo(N).int.signedness != .unsigned) {
        @compileError("`a` and `b` must be unsigned integers");
    }

    std.debug.assert(a != 0 or b != 0);

    if (a == 0) return b;
    if (b == 0) return a;

    var x: N = a;
    var y: N = b;

    const xz = @ctz(x);
    const yz = @ctz(y);
    const shift = @min(xz, yz);
    x = @shrExact(x, @intCast(xz));
    y = @shrExact(y, @intCast(yz));

    var y_minus_x = y -% x;
    while (y_minus_x != 0) : (y_minus_x = y -% x) {
        const copy_x = x;
        const zeros = @ctz(y_minus_x);
        const carry = x < y;
        x -%= y;
        if (carry) {
            x = y_minus_x;
            y = copy_x;
        }
        x = @shrExact(x, @intCast(zeros));
    }

    return @shlExact(y, @intCast(shift));
}

const TestCase = struct {
    name: []const u8,
    pairs: []const [2]u64,
};

fn generateFibonacciPairs(allocator: std.mem.Allocator, count: usize) ![][2]u64 {
    const pairs = try allocator.alloc([2]u64, count);
    var a: u64 = 1;
    var b: u64 = 1;

    for (pairs) |*pair| {
        const next = a +% b;
        pair.* = .{ b, a };
        a = b;
        b = next;
    }

    return pairs;
}

fn generateRandomPairs(allocator: std.mem.Allocator, count: usize, seed: u64) ![][2]u64 {
    var prng = std.Random.DefaultPrng.init(seed);
    const random = prng.random();

    const pairs = try allocator.alloc([2]u64, count);
    for (pairs) |*pair| {
        pair.* = .{
            random.int(u64) | 1, // Ensure non-zero
            random.int(u64) | 1,
        };
    }

    return pairs;
}

fn benchmarkGcd(
    comptime gcd_fn: anytype,
    pairs: []const [2]u64,
    iterations: usize,
) !u64 {
    var timer = try std.time.Timer.start();

    const start = timer.read();
    for (0..iterations) |_| {
        for (pairs) |pair| {
            var result = gcd_fn(pair[0], pair[1]);
            std.mem.doNotOptimizeAway(&result);
        }
    }
    const end = timer.read();

    return end - start;
}

pub fn main() !void {
    var gpa = std.heap.GeneralPurposeAllocator(.{}){};
    defer _ = gpa.deinit();
    const allocator = gpa.allocator();

    var stdout_buffer: [8192]u8 = undefined;
    var stdout_writer = std.fs.File.stdout().writer(&stdout_buffer);
    const stdout = &stdout_writer.interface;

    try stdout.print("GCD Performance Benchmark\n", .{});
    try stdout.print("==========================\n\n", .{});

    const test_cases = [_]TestCase{
        .{
            .name = "Small Numbers",
            .pairs = &[_][2]u64{
                .{ 8, 12 },
                .{ 12, 8 },
                .{ 33, 77 },
                .{ 100, 50 },
                .{ 48, 18 },
                .{ 1071, 462 },
            },
        },
        .{
            .name = "Powers of 2",
            .pairs = &[_][2]u64{
                .{ 1024, 2048 },
                .{ 4096, 8192 },
                .{ 16384, 32768 },
                .{ 65536, 131072 },
                .{ 1 << 20, 1 << 21 },
            },
        },
        .{
            .name = "High Trailing Zeros",
            .pairs = &[_][2]u64{
                .{ 300_000, 2_300_000 },
                .{ 1_000_000, 5_000_000 },
                .{ 12_000_000, 18_000_000 },
            },
        },
        .{
            .name = "Coprime Numbers",
            .pairs = &[_][2]u64{
                .{ 97, 101 },
                .{ 1009, 1013 },
                .{ 10007, 10009 },
                .{ 100003, 100019 },
            },
        },
        .{
            .name = "Large Numbers",
            .pairs = &[_][2]u64{
                .{ 49865, 69811 },
                .{ 1299827, 1988891 },
                .{ 7919 * 7927, 7919 * 7933 },
                .{ 1_000_000_007, 1_000_000_009 },
            },
        },
    };

    const warmup_iterations: usize = 10_000;
    const bench_iterations: usize = 100_000;

    for (test_cases) |test_case| {
        try stdout.print("Test Case: {s}\n", .{test_case.name});
        try stdout.print("{s}\n", .{"-" ** 50});

        // Warmup
        _ = try benchmarkGcd(gcd_old, test_case.pairs, warmup_iterations);
        _ = try benchmarkGcd(gcd_new, test_case.pairs, warmup_iterations);

        // Benchmark multiple runs for statistical significance
        const num_runs = 10;
        var old_times: [num_runs]u64 = undefined;
        var new_times: [num_runs]u64 = undefined;

        for (0..num_runs) |i| {
            old_times[i] = try benchmarkGcd(gcd_old, test_case.pairs, bench_iterations);
            new_times[i] = try benchmarkGcd(gcd_new, test_case.pairs, bench_iterations);
        }

        // Calculate statistics
        var old_sum: u64 = 0;
        var new_sum: u64 = 0;
        var old_min: u64 = std.math.maxInt(u64);
        var old_max: u64 = 0;
        var new_min: u64 = std.math.maxInt(u64);
        var new_max: u64 = 0;

        for (old_times, new_times) |old_time, new_time| {
            old_sum += old_time;
            new_sum += new_time;
            old_min = @min(old_min, old_time);
            old_max = @max(old_max, old_time);
            new_min = @min(new_min, new_time);
            new_max = @max(new_max, new_time);
        }

        const old_avg = old_sum / num_runs;
        const new_avg = new_sum / num_runs;

        try stdout.print("  Old implementation:\n", .{});
        try stdout.print("    Min: {d:.3} ms\n", .{@as(f64, @floatFromInt(old_min)) / 1_000_000.0});
        try stdout.print("    Avg: {d:.3} ms\n", .{@as(f64, @floatFromInt(old_avg)) / 1_000_000.0});
        try stdout.print("    Max: {d:.3} ms\n", .{@as(f64, @floatFromInt(old_max)) / 1_000_000.0});

        try stdout.print("  New implementation:\n", .{});
        try stdout.print("    Min: {d:.3} ms\n", .{@as(f64, @floatFromInt(new_min)) / 1_000_000.0});
        try stdout.print("    Avg: {d:.3} ms\n", .{@as(f64, @floatFromInt(new_avg)) / 1_000_000.0});
        try stdout.print("    Max: {d:.3} ms\n", .{@as(f64, @floatFromInt(new_max)) / 1_000_000.0});

        const diff_pct = if (old_avg > 0)
            (@as(f64, @floatFromInt(new_avg)) - @as(f64, @floatFromInt(old_avg))) / @as(f64, @floatFromInt(old_avg)) * 100.0
        else
            0.0;

        if (diff_pct < -0.5) {
            try stdout.print("  Result: NEW IS FASTER by {d:.2}%\n", .{-diff_pct});
        } else if (diff_pct > 0.5) {
            try stdout.print("  Result: OLD IS FASTER by {d:.2}%\n", .{diff_pct});
        } else {
            try stdout.print("  Result: ROUGHLY EQUAL (diff: {d:.2}%)\n", .{diff_pct});
        }

        try stdout.print("\n", .{});
    }

    // Fibonacci numbers (worst case for Euclidean algorithm)
    try stdout.print("Test Case: Fibonacci Numbers (Worst Case)\n", .{});
    try stdout.print("{s}\n", .{"-" ** 50});

    const fib_pairs = try generateFibonacciPairs(allocator, 20);
    defer allocator.free(fib_pairs);

    _ = try benchmarkGcd(gcd_old, fib_pairs, warmup_iterations);
    _ = try benchmarkGcd(gcd_new, fib_pairs, warmup_iterations);

    const old_fib = try benchmarkGcd(gcd_old, fib_pairs, bench_iterations);
    const new_fib = try benchmarkGcd(gcd_new, fib_pairs, bench_iterations);

    try stdout.print("  Old: {d:.3} ms\n", .{@as(f64, @floatFromInt(old_fib)) / 1_000_000.0});
    try stdout.print("  New: {d:.3} ms\n", .{@as(f64, @floatFromInt(new_fib)) / 1_000_000.0});

    const fib_diff = (@as(f64, @floatFromInt(new_fib)) - @as(f64, @floatFromInt(old_fib))) / @as(f64, @floatFromInt(old_fib)) * 100.0;
    try stdout.print("  Difference: {d:.2}%\n", .{fib_diff});
    try stdout.print("\n", .{});

    // Random numbers
    try stdout.print("Test Case: Random Numbers\n", .{});
    try stdout.print("{s}\n", .{"-" ** 50});

    const rand_pairs = try generateRandomPairs(allocator, 100, 12345);
    defer allocator.free(rand_pairs);

    _ = try benchmarkGcd(gcd_old, rand_pairs, warmup_iterations / 10);
    _ = try benchmarkGcd(gcd_new, rand_pairs, warmup_iterations / 10);

    const old_rand = try benchmarkGcd(gcd_old, rand_pairs, bench_iterations / 10);
    const new_rand = try benchmarkGcd(gcd_new, rand_pairs, bench_iterations / 10);

    try stdout.print("  Old: {d:.3} ms\n", .{@as(f64, @floatFromInt(old_rand)) / 1_000_000.0});
    try stdout.print("  New: {d:.3} ms\n", .{@as(f64, @floatFromInt(new_rand)) / 1_000_000.0});

    const rand_diff = (@as(f64, @floatFromInt(new_rand)) - @as(f64, @floatFromInt(old_rand))) / @as(f64, @floatFromInt(old_rand)) * 100.0;
    try stdout.print("  Difference: {d:.2}%\n", .{rand_diff});

    try stdout.flush();
}