binary GCD algorithm

{{Short description|Algorithm for computing the greatest common divisor}}

File:binary_GCD_algorithm_visualisation.svg

The binary GCD algorithm, also known as Stein's algorithm or the binary Euclidean algorithm,{{r|brenta|brentb}} is an algorithm that computes the greatest common divisor (GCD) 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 physicist and programmer Josef Stein in 1967, it was known by the 2nd century BCE, in ancient China.{{r|Knuth98}}

Algorithm

The algorithm finds the GCD of two nonnegative numbers $u$ and $v$ by repeatedly applying these identities:

$\gcd(u, 0) = u$ : everything divides zero, and $u$ is the largest number that divides $u$ .
$\gcd(2u, 2v) = 2 \cdot \gcd(u, v)$ : $2$ is a common divisor.
$\gcd(u, 2v) = \gcd(u, v)$ if $u$ is odd: $2$ is then not a common divisor.
$\gcd(u, v) = \gcd(u, v - u)$ if $u, v$ odd and $u \leq v$ .

As GCD is commutative ( $\gcd(u, v) = \gcd(v, u)$ ), those identities still apply if the operands are swapped: $\gcd(0, v) = v$ , $\gcd(2u, v) = \gcd(u, v)$ if $v$ is odd, etc.

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 $2$ 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 iteratively rather than recursively: the resulting implementation can be laid out to avoid repeated work, invoking identity 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-free except for its exit condition ( $v = 0$ ): in the example below, the exchange of $u$ and $v$ (ensuring $u \leq v$ ) compiles down to conditional moves; hard-to-predict branches can have a large, negative impact on performance.

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

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();

}

Note: The implementation above accepts unsigned (non-negative) integers; given that $\gcd(u, v) = \gcd(\pm{}u, \pm{}v)$ , the signed case can be handled as follows:

/// 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())

}

Complexity

Asymptotically, the algorithm requires $O(n)$ steps, where $n$ is the number of bits in the larger of the two numbers, as every two steps reduce 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); when working with word-sized numbers, each arithmetic operation translates to a single machine operation, so the number of machine operations is on the order of $n$ , i.e. $\log_{2}(\max(u, v))$ .

For arbitrarily large numbers, the asymptotic complexity of this algorithm is $O(n^2)$ , as each arithmetic operation (subtract and shift) involves a linear number of machine operations (one per word in the numbers' binary representation).

If the numbers can be represented in the machine's memory, i.e. each number's size can be represented by a single machine word, this bound is reduced to:

$O\left(\frac{n^2}{\log_2 n}\right)$

This is the same as for the Euclidean algorithm, though a more precise analysis by Akhavi and Vallée proved that binary GCD uses about 60% fewer bit operations.

Extensions

The binary GCD algorithm can be extended in several ways, either to output additional information, deal with arbitrarily large integers more efficiently, or to compute GCDs in domains other than the integers.

The extended binary GCD algorithm, analogous to the extended Euclidean algorithm, fits in the first kind of extension, as it provides the Bézout coefficients in addition to the GCD: integers $a$ and $b$ such that $a\cdot{}u + b\cdot{}v = \gcd(u, v)$ .

In the case of large integers, the best asymptotic complexity is $O(M(n) \log n)$ , with $M(n)$ the cost of $n$ -bit multiplication; this is near-linear and vastly smaller than the binary GCD algorithm's $O(n^2)$ , though concrete implementations only outperform older algorithms for numbers larger than about 64 kilobits (i.e. greater than 8×10¹⁹²⁶⁵). This is achieved by extending the binary GCD algorithm using ideas from the Schönhage–Strassen algorithm for fast integer multiplication.

The binary GCD algorithm has also been extended to domains other than natural numbers, such as Gaussian integers, Eisenstein integers, quadratic rings, and integer rings of number fields.

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:

{{quote| title= Fangtian – Land surveying| source=The Nine Chapters on the Mathematical Art |text=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.}}

The phrase "if possible halve it" is ambiguous,

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.

References

{{Citation |first=J. |last=Stein |title=Computational problems associated with Racah algebra |journal=Journal of Computational Physics |volume= 1 |issue=3 |pages=397–405 |date=February 1967 |issn=0021-9991 |doi=10.1016/0021-9991(67)90047-2|bibcode=1967JCoPh...1..397S }}

{{Citation |first=Donald |last=Knuth |author-link=Donald Knuth |series=The Art of Computer Programming |volume=2 |title=Seminumerical Algorithms |edition=3rd |publisher=Addison-Wesley |isbn=978-0-201-89684-8 |year=1998 }}

{{Cite web | title=Avoiding the Cost of Branch Misprediction | first=Rajiv | last=Kapoor | date=21 February 2009 | website=Intel Developer Zone | url=https://software.intel.com/content/www/us/en/develop/articles/avoiding-the-cost-of-branch-misprediction.html }}

{{Cite web | title=Mispredicted branches can multiply your running times | first=Daniel | last=Lemire | date=15 October 2019 | url=https://lemire.me/blog/2019/10/15/mispredicted-branches-can-multiply-your-running-times/ }}

{{Cite web | title=Compiler Explorer | first=Matt | last=Godbolt | access-date=4 February 2024 | url=https://rust.godbolt.org/z/56jva3KPn }}

{{Harvnb|Knuth|1998|p=646}}, answer to exercise 39 of section 4.5.2

{{Cite book |chapter=§14.4 Greatest Common Divisor Algorithms |chapter-url=http://cacr.uwaterloo.ca/hac/about/chap14.pdf#page=17 |title=Handbook of Applied Cryptography |url=http://cacr.uwaterloo.ca/hac/ |pages=606–610 |date=October 1996 |publisher=CRC Press |first1=Alfred J. |last1=Menezes |first2=Paul C. |last2=van Oorschot |first3=Scott A. |last3=Vanstone |isbn=0-8493-8523-7 |access-date=9 September 2017}}

{{cite book |last=Cohen |first=Henri |author-link= Henri Cohen (number theorist) |year=1993 |title=A Course In Computational Algebraic Number Theory

|pages=17–18 |chapter=Chapter 1 : Fundamental Number-Theoretic Algorithms

|isbn= 0-387-55640-0|publisher=Springer-Verlag|series=Graduate Texts in Mathematics|volume=138

|url= https://books.google.com/books?id=hXGr-9l1DXcC}}

{{Cite web | url=http://gmplib.org/manual/Binary-GCD.html | title=GNU MP 6.1.2: Binary GCD}}

{{citation

| last1 = Stehlé | first1 = Damien

| last2 = Zimmermann | first2 = Paul | author-link2 = Paul Zimmermann (mathematician)

| contribution = A binary recursive gcd algorithm

| doi = 10.1007/978-3-540-24847-7_31

| mr = 2138011

| pages = 411–425

| publisher = Springer, Berlin

| series = Lecture Notes in Comput. Sci.

| title = Algorithmic number theory

| contribution-url = http://hal.archives-ouvertes.fr/docs/00/07/15/33/PDF/RR-5050.pdf

| volume = 3076

| year = 2004 | isbn = 978-3-540-22156-2

| citeseerx = 10.1.1.107.8612

| s2cid = 3119374

| id = INRIA Research Report RR-5050

}}.

{{Citation | first1 = Ali | last1 = Akhavi | first2 = Brigitte | last2 = Vallée | title = Average Bit-Complexity of Euclidean Algorithms | year = 2000 | journal = Proceedings ICALP'00, Lecture Notes Computer Science 1853 | pages = 373–387 | url = https://vallee.users.greyc.fr/Publications/icalp8-2000.ps | citeseerx = 10.1.1.42.7616}}

{{cite conference |first=Richard P. |last=Brent |author-link=Richard P. Brent |url=http://wwwmaths.anu.edu.au/~brent/pub/pub183.html |title=Twenty years' analysis of the Binary Euclidean Algorithm |conference=1999 Oxford-Microsoft Symposium in honour of Professor Sir Antony Hoare |date=13–15 September 1999 |location=Oxford}}

{{cite tech report |first=Richard P. |last=Brent |author-link=Richard P. Brent |title=Further analysis of the Binary Euclidean algorithm |publisher=Oxford University Computing Laboratory |id=PRG TR-7-99 |date=November 1999 |arxiv=1303.2772 |url=https://www.cs.ox.ac.uk/techreports/oucl/tr-7-99.html}}

{{cite journal |first=André |last=Weilert |title=(1+i)-ary GCD Computation in Z[i] as an Analogue to the Binary GCD Algorithm

|date=July 2000 |journal=Journal of Symbolic Computation |volume=30 |issue=5 |pages=605–617 |doi=10.1006/jsco.2000.0422|doi-access=free }}

{{cite conference |last1=Damgård |first1=Ivan Bjerre |last2=Frandsen |first2=Gudmund Skovbjerg

|title=Efficient Algorithms for GCD and Cubic Residuosity in the Ring of Eisenstein Integers |doi=10.1007/978-3-540-45077-1_11

|doi-access= |pages=109–117 |conference=14th International Symposium on the Fundamentals of Computation Theory |location=Malmö, Sweden |date= 12–15 August 2003}}

{{cite conference |last1=Agarwal |first1=Saurabh |last2=Frandsen |first2=Gudmund Skovbjerg

|title=Binary GCD Like Algorithms for Some Complex Quadratic Rings |doi=10.1007/978-3-540-24847-7_4 |pages=57–71

|conference=Algorithmic Number Theory Symposium |date= 13–18 June 2004 |location=Burlington, VT, USA}}

{{cite conference |last1=Agarwal |first1=Saurabh |last2=Frandsen |first2=Gudmund Skovbjerg

|title=A New GCD Algorithm for Quadratic Number Rings with Unique Factorization |doi=10.1007/11682462_8 |pages=30–42

|conference=7th Latin American Symposium on Theoretical Informatics |date= 20–24 March 2006 |location=Valdivia, Chile}}

{{cite conference |last=Wikström |first=Douglas

|title=On the l-Ary GCD-Algorithm in Rings of Integers |doi=10.1007/11523468_96 |pages=1189–1201 |date= 11–15 July 2005

|conference=Automata, Languages and Programming, 32nd International Colloquium |location=Lisbon, Portugal}}

External links

[https://xlinux.nist.gov/dads/HTML/binaryGCD.html NIST Dictionary of Algorithms and Data Structures: binary GCD algorithm]
[http://www.cut-the-knot.org/blue/binary.shtml Cut-the-Knot: Binary Euclid's Algorithm] at cut-the-knot
[http://wwwmaths.anu.edu.au/~brent/pub/pub037.html Analysis of the Binary Euclidean Algorithm] (1976), a paper by Richard P. Brent, including a variant using left shifts

Category:Number theoretic algorithms

Category:Articles with example Rust code