Visualisation of using the binary GCD algorithm to find the greatest common divisor (GCD) of 36 and 24. Thus, the GCD is 2² × 3 = 12.

The binary GCD algorithm, also known as Stein's algorithm, is an algorithm that computes the greatest common divisor of two nonnegative integers. Stein's algorithm uses simpler arithmetic operations than the conventional Euclidean algorithm; it replaces division with arithmetic shifts, comparisons, and subtraction.

Although the algorithm in its contemporary form was first published by the Israeli physicist and programmer Josef Stein in 1967,^[1] it may have been known by the 2nd century BCE, in ancient China.^[2][3]

Algorithm

The algorithm reduces the problem of finding the GCD of two nonnegative numbers v and u by repeatedly applying these identities:

• gcd(0, v) = v, because everything divides zero, and v is the largest number that divides v. Similarly, gcd(u, 0) = u.
• gcd(2u, 2v) = 2·gcd(u, v)
• gcd(2u, v) = gcd(u, v), if v is odd (2 is not a common divisor). Similarly, gcd(u, 2v) = gcd(u, v) if u is odd.
• gcd(u, v) = gcd(|u − v|, min(u, v)), if u and v are both odd.

An extended binary GCD, analogous to the extended Euclidean algorithm, can provide the Bézout coefficients in addition to the GCD, i.e. integers a and b such that a·u + b·v = gcd(u, v).^[4][5]

Implementation

Recursive version in C

Following is a recursive implementation of the algorithm in C. The implementation is similar to the description of the algorithm given above, and optimised for readability rather than speed, though all but one of the recursive calls are tail recursive.

unsigned int gcd(unsigned int u, unsigned int v)
{
    // Base cases
    //  gcd(n, n) = n
    if (u == v)
        return u;
    
    //  Identity 1: gcd(0, n) = gcd(n, 0) = n
    if (u == 0)
        return v;
    if (v == 0)
        return u;

    if (u % 2 == 0) { // u is even
        if (v % 2 == 1) // v is odd
            return gcd(u/2, v); // Identity 3
        else // both u and v are even
            return 2 * gcd(u/2, v/2); // Identity 2

    } else { // u is odd
        if (v % 2 == 0) // v is even
            return gcd(u, v/2); // Identity 3

        // Identities 4 and 3 (u and v are odd, so u-v and v-u are known to be even)
        if (u > v)
            return gcd((u - v)/2, v);
        else
            return gcd((v - u)/2, u);
    }
}

Iterative version in Rust

Following is an implementation of the algorithm in Rust, adapted from uutils. It removes all common factors of 2, using identity 2, then computes the GCD of the remaining numbers using identities 3 and 4, combining these to form the final answer.

pub fn gcd(mut u: u64, mut v: u64) -> u64 {
    use std::cmp::min;
    use std::mem::swap;

    // 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 u and 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 = {} is even", u);
        debug_assert!(v % 2 == 1, "v = {} is even", v);

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

        // Using identity 4 (gcd(u, v) = gcd(|v-u|, min(u, v))
        v -= u;

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

        // Identity 3: gcd(u, 2ʲ v) = gcd(u, v) (u is known to be odd)
        v >>= v.trailing_zeros();
    }
}

This implementation showcases several performance optimisations:

• Trial division by 2 is eschewed in favour of a single bitshift and the count trailing zeros primitive u64::trailing_zeros.
• The loop is laid out so as to avoid repeated work; for instance, eliminating 2 as a factor of v was moved to the back of the loop, and after the exit condition, as v is known to be odd upon entering the loop.
• The body of the loop is branch-free except for its exit condition (v == 0), as the exchange of u and v (ensuring u ≤ v) compiles down to conditional moves.^[6] Hard-to-predict branches can have a large, negative impact on performance.^[7][8]

Efficiency

The algorithm requires O(n) steps, where n is the number of bits in the larger of the two numbers, as every step reduces at least one of the operands by at least a factor of 2. Each step involves only a few arithmetic operations (O(1) with a small constant), so the number of machine operations is on the order of log(max(u, v)) when working with word-sized numbers.

However, the asymptotic complexity of this algorithm is O(n²),^[9] as those arithmetic operations (subtract and shift) each take linear time for arbitrarily-sized numbers (one machine operation per word of the representation). This is the same as for the Euclidean algorithm, though neither is the fastest for arbitrary-precision arithmetic; instead, recursive methods that combine ideas from the binary GCD algorithm with the Schönhage–Strassen algorithm for fast integer multiplication can find GCDs in near-linear time, but only outperform older algorithms for numbers larger than about 64 kilobits (i.e. greater than 8×10¹⁹²⁶⁵).^[10]

A more precise analysis by Akhavi and Vallée proved that binary GCD uses about 60% fewer bit operations than the Euclidean algorithm.^[11]

Historical description

An algorithm for computing the GCD of two numbers was known in ancient China, under the Han dynasty, as a method to reduce fractions:

If possible halve it; otherwise, take the denominator and the numerator, subtract the lesser from the greater, and do that alternately to make them the same. Reduce by the same number.

— Fangtian - Land surveying, The Nine Chapters on the Mathematical Art

The phrase "if possible halve it" is ambiguous,^[2][3]

• if this applies when either of the numbers become even, the algorithm is the binary GCD algorithm;
• if this only applies when both numbers are even, the algorithm is similar to the Euclidean algorithm.

See also

• Euclidean algorithm
• Extended Euclidean algorithm
• Least common multiple

References

Further reading

• Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein. Introduction to Algorithms, Second Edition. MIT Press and McGraw-Hill, 2001. ISBN 0-262-03293-7. Problem 31-1, pg.902.
• Richard P. Brent, Twenty years' analysis of the Binary Euclidean Algorithm, 2000

External links

• NIST Dictionary of Algorithms and Data Structures: binary GCD algorithm
• Cut-the-Knot: Binary Euclid's Algorithm at cut-the-knot
• Analysis of the Binary Euclidean Algorithm (1976), a paper by Richard P. Brent, including a variant using left shifts
• "Dynamics of the Binary Euclidean Algorithm: Functional Analysis and Operators" (1998), a paper by Brigitte Vallée.