Rust benchmark of Binary GCD algorithms

gcd.rs

//! Benchmarking GCD algorithms on primitive integers
//! 
//! Time consumption with my test (on a 64bit system):
//! u32 : gcd_bin < gcd_bin_bfree < gcd_bin_compact < gcd_euclid; gcd_euclid_ext < gcd_bin_ext;
//! u64 : gcd_bin_compact < gcd_bin_bfree < gcd_bin < gcd_euclid; gcd_euclid_ext < gcd_bin_ext;
//! u128: gcd_bin_compact < gcd_bin_bfree < gcd_bin < gcd_euclid; gcd_euclid_ext < gcd_bin_ext;

#[macro_use]
extern crate criterion;
use criterion::Criterion;
use rand::random;

#[allow(non_camel_case_types)]
type utarget = u64;
#[allow(non_camel_case_types)]
type itarget = i64;

fn gcd_euclid(mut a: utarget, mut b: utarget) -> utarget {
    let mut r = a % b;
    while r > 0 {
        a = b % r;
        b = r;
        r = a;
    }
    b
}

fn gcd_euclid_ext(a: utarget, b: utarget) -> (utarget, itarget, itarget) {
    let (mut last_r, mut r) = (a, b);
    let (mut last_s, mut s) = (1, 0);
    let (mut last_t, mut t) = (0, 1);

    while r > 0 {
        let quo = last_r / r;
        let new_r = last_r - quo*r;
        last_r = std::mem::replace(&mut r, new_r);
        let new_s = last_s - quo as itarget*s;
        last_s = std::mem::replace(&mut s, new_s);
        let new_t = last_t - quo as itarget*t;
        last_t = std::mem::replace(&mut t, new_t);
    }

    (last_r, last_s, last_t)
}

// Use Stein's algorithm
fn gcd_bin(mut m: utarget, mut n: utarget) -> utarget {
    if m == 0 || n == 0 {
        return m | n;
    }

    // find common factors of 2
    let shift = (m | n).trailing_zeros();
    m >>= m.trailing_zeros();
    n >>= n.trailing_zeros();

    while m != n {
        if m > n {
            m -= n;
            m >>= m.trailing_zeros();
        } else {
            n -= m;
            n >>= n.trailing_zeros();
        }
    }
    m << shift
}

// Use Stein's algorithm
fn gcd_bin_compact(mut m: utarget, mut n: utarget) -> utarget {
    if m == 0 || n == 0 {
        return m | n;
    }

    // find common factors of 2
    let shift = (m | n).trailing_zeros();

    // divide n and m by 2 until odd
    // m inside loop
    n >>= n.trailing_zeros();

    while m != 0 {
        m >>= m.trailing_zeros();
        if n > m {
            std::mem::swap(&mut n, &mut m)
        }
        m -= n;
    }

    n << shift  
}

// Branch free version
fn gcd_bin_bfree(m: utarget, mut n: utarget) -> utarget {
    if m == 0 || n == 0 {
        return m | n;
    }

    // find common factors of 2
    let shift = (m | n).trailing_zeros();

    // divide n and m by 2 until odd
    // m inside loop
    n >>= n.trailing_zeros();
    
    // In this loop, we represent the odd numbers u and v
    // without the redundant least significant one bit. This reduction
    // in size by one bit ensures that the high bit of t, below, is set
    // if and only if v > u.

    let mut u = (m >> 1) as itarget;
    let mut v = (n >> 1) as itarget;

    while u != v {
        let t = u - v;
        let vgtu = t >> (utarget::BITS - 1);

        /* v <-- min (u, v) */
        v += vgtu & t;

        /* u <-- |u - v| */
        u = (t ^ vgtu) - vgtu;
        u = u >> (t.trailing_zeros() + 1);
    }
    let g = ((u << 1) + 1) as utarget;
    g << shift
}

// Extended Binary GCD, From GMP
fn gcd_bin_ext(mut u: utarget, mut v: utarget) -> (utarget, itarget, itarget) {
    if u == 0 {
        return (v, 0, 1);
    }
    if v == 0 {
        return (u, 1, 0);
    }
    let (mut s0, mut s1) = (1 as itarget, 0);
    let (mut t0, mut t1) = (0, 1 as itarget);

    // fivd common factors of 2
    let zeros = (u | v).trailing_zeros();

    // divide v and u by 2 until odd
    u >>= zeros;
    v >>= zeros;
    let mut shift;
    if (u & 1) == 0 {
        shift = u.trailing_zeros();
        u >>= shift;
        t1 <<= shift;
    }
    else if (v & 1) == 0 {
        shift = v.trailing_zeros();
        v >>= shift;
        s0 <<= shift;
    }
    else {
        shift = 0
    };

    // main loop
    while u != v {
        let count;
        if u > v {
            u -= v;
            count = u.trailing_zeros();
            u >>= count;
            t0 += t1; t1 <<= count;
            s0 += s1; s1 <<= count;
        } else {
            v -= u;
            count = v.trailing_zeros();
            v >>= count;
            t1 += t0; t0 <<= count;
            s1 += s0; s0 <<= count;
        }
        shift += count;
    }

    // now u = v = g = gcd (u,v). Compute U/g and V/g
    let ug = t0 + t1;
    let vg = s0 + s1;

    // now 2^{shift} g = s0 U - t0 V. Get rid of the power of two, using
    // s0 U - t0 V = (s0 + V/g) U - (t0 + U/g) V.
    for _ in 0..shift {
        if (s0 | t0) & 1 > 0 {            
            s0 += vg;
            t0 += ug;
        }
        s0 /= 2;
        t0 /= 2;
    }
    if s0 > vg - s0 {
        s0 -= vg;
        t0 -= ug;
    }

    (u << zeros, s0, -t0)
}

fn bench_gcd(c: &mut Criterion) {
    let mut group = c.benchmark_group("gcd");

    const N: usize = 200;
    let mut lhs: [utarget; N] = [0; N];
    let mut rhs: [utarget; N] = [0; N];
    for i in 0..N {
        lhs[i] = random();
        rhs[i] = random();
    }

    assert_eq!(gcd_euclid(5, 2), 1);
    assert_eq!(gcd_euclid(10, 4), 2);
    group.bench_function("euclid", |b| {
        b.iter(|| {
            lhs.iter()
                .zip(rhs.iter())
                .map(|(&a, &b)| gcd_euclid(a, b))
                .reduce(|a, b| a.wrapping_add(b))
        })
    });

    assert_eq!(gcd_bin_compact(5, 2), 1);
    assert_eq!(gcd_bin_compact(10, 4), 2);
    group.bench_function("binary (compact)", |b| {
        b.iter(|| {
            lhs.iter()
                .zip(rhs.iter())
                .map(|(&a, &b)| gcd_bin_compact(a, b))
                .reduce(|a, b| a.wrapping_add(b))
        })
    });
    
    assert_eq!(gcd_bin(5, 2), 1);
    assert_eq!(gcd_bin(10, 4), 2);
    group.bench_function("binary", |b| {
        b.iter(|| {
            lhs.iter()
                .zip(rhs.iter())
                .map(|(&a, &b)| gcd_bin(a, b))
                .reduce(|a, b| a.wrapping_add(b))
        })
    });
    
    assert_eq!(gcd_bin_bfree(5, 2), 1);
    assert_eq!(gcd_bin_bfree(10, 4), 2);
    group.bench_function("binary (branch free)", |b| {
        b.iter(|| {
            lhs.iter()
                .zip(rhs.iter())
                .map(|(&a, &b)| gcd_bin_bfree(a, b))
                .reduce(|a, b| a.wrapping_add(b))
        })
    });
    group.finish();

    let mut group = c.benchmark_group("extended gcd");

    assert_eq!(gcd_euclid_ext(5, 2), (1, 1, -2));
    assert_eq!(gcd_euclid_ext(10, 4), (2, 1, -2));
    group.bench_function("euclid", |b| {
        b.iter(|| {
            lhs.iter()
                .zip(rhs.iter())
                .map(|(&a, &b)| gcd_euclid_ext(a, b))
                .reduce(|a, b| (a.0.wrapping_add(b.0), a.1.wrapping_add(b.1), a.2.wrapping_add(b.2)))
        })
    });
    
    assert_eq!(gcd_bin_ext(5, 2), (1, 1, -2));
    assert_eq!(gcd_bin_ext(10, 4), (2, 1, -2));
    group.bench_function("binary", |b| {
        b.iter(|| {
            lhs.iter()
                .zip(rhs.iter())
                .map(|(&a, &b)| gcd_bin_ext(a, b))
                .reduce(|a, b| (a.0.wrapping_add(b.0), a.1.wrapping_add(b.1), a.2.wrapping_add(b.2)))
        })
    });
}

criterion_group!(benches, bench_gcd);
criterion_main!(benches);