Editing Binary GCD algorithm (section)

==Implementation==

While the above description of the algorithm is mathematically correct, performant software implementations typically differ from it in a few notable ways:
* eschewing [[trial division]] by <math>2</math> in favour of a single bitshift and the [[count trailing zeros]] primitive; this is functionally equivalent to repeatedly applying identity 3, but much faster;
* expressing the algorithm [[Iteration#Computing|iteratively]] rather than [[Recursion (computer science)|recursively]]: the resulting implementation can be laid out to avoid repeated work, invoking identity&nbsp;2 at the start and maintaining as invariant that both numbers are odd upon entering the loop, which only needs to implement identities 3 and 4;
* making the loop's body [[Branch_(computer_science)#Branch-free_code|branch-free]] except for its exit condition (<math>v = 0</math>): in the example below, the exchange of <math>u</math> and <math>v</math> (ensuring <math>u \leq v</math>) compiles down to [[conditional moves]];<ref name="rust disassembly"/> hard-to-predict branches can have a large, negative impact on performance.<ref name="intel"/><ref name="lemire"/>

The following is an implementation of the algorithm in [[Rust (programming language)|Rust]] exemplifying those differences, adapted from [https://github.com/uutils/coreutils/blob/1eabda91cf35ec45c78cb95c77d5463607daed65/src/uu/factor/src/numeric/gcd.rs ''uutils'']:

<!-- PLEASE CHECK ANY CHANGES TO THIS ALGORITHM WITH THE TESTS AND OUTPUT HERE:
     https://play.rust-lang.org/?version=stable&mode=debug&edition=2021&gist=ea08036e202de879eb462dd5fc9747b3 -->
<syntaxhighlight lang="rust">
use std::cmp::min;
use std::mem::swap;

pub fn gcd(mut u: u64, mut v: u64) -> u64 {
    // Base cases: gcd(n, 0) = gcd(0, n) = n
    if u == 0 {
        return v;
    } else if v == 0 {
        return u;
    }

    // Using identities 2 and 3:
    // gcd(2ⁱ u, 2ʲ v) = 2ᵏ gcd(u, v) with u, v odd and k = min(i, j)
    // 2ᵏ is the greatest power of two that divides both 2ⁱ u and 2ʲ v
    let i = u.trailing_zeros();  u >>= i;
    let j = v.trailing_zeros();  v >>= j;
    let k = min(i, j);

    loop {
        // u and v are odd at the start of the loop
        debug_assert!(u % 2 == 1, "u = {} should be odd", u);
        debug_assert!(v % 2 == 1, "v = {} should be odd", v);

        // Swap if necessary so u ≤ v
        if u > v {
            swap(&mut u, &mut v);
        }

        // Identity 4: gcd(u, v) = gcd(u, v-u) as u ≤ v and u, v are both odd 
        v -= u;
        // v is now even

        if v == 0 {
            // Identity 1: gcd(u, 0) = u
            // The shift by k is necessary to add back the 2ᵏ factor that was removed before the loop
            return u << k;
        }

        // Identity 3: gcd(u, 2ʲ v) = gcd(u, v) as u is odd
        v >>= v.trailing_zeros();
    }
}
</syntaxhighlight>

'''Note''': The implementation above accepts ''unsigned'' (non-negative) integers; given that <math>\gcd(u, v) = \gcd(\pm{}u, \pm{}v)</math>, the signed case can be handled as follows:
<syntaxhighlight lang="rust">
/// Computes the GCD of two signed 64-bit integers
/// The result is unsigned and not always representable as i64: gcd(i64::MIN, i64::MIN) == 1 << 63
pub fn signed_gcd(u: i64, v: i64) -> u64 {
    gcd(u.unsigned_abs(), v.unsigned_abs())
}
</syntaxhighlight>