Lehmer's GCD algorithm

Lehmer's GCD algorithm, named after D. H. Lehmer, is a fast GCD algorithm for multiple-precision arithmetic, which improves on the simpler Euclidean algorithm by doing most operations using only the leading digits of the values. "Digit" here is not necessarily a decimal digit; the algorithm is often used with integers represented using a base β such as β = 1000 or β = 2³².

Algorithm

Lehmer noted that most of the quotients from each step of the division part of the standard algorithm are small. (For example, Knuth observed that the quotients 1, 2, and 3 comprise 67.7% of all quotients.) Those small quotients can be identified from only a few leading digits. Thus the algorithm starts by splitting off those leading digits and computing the sequence of quotients as long as it is correct.

Say we want to obtain the GCD of the two integers a and b. Let a ≥ b.

• If b contains only one digit (in the chosen base, say β = 1000 or β = 2³²), use some other method, such as the Euclidean algorithm, to obtain the result.
• If a and b differ in the length of digits, reduce a modulo b and exchange them as in the standard Euclidean algorithm. Repeat until the two are of the same length m.
• Outer loop: Iterate until one of a or b is zero:
  • Decrease m by one. Let x be the leading (most significant) digit in a, x = a div β^m and y the leading digit in b, y = b div β^m.
  • Initialize a 2 by 3 matrix

  $\textstyle {\begin{bmatrix}A&B&x\\C&D&y\end{bmatrix}}$ to an extended identity matrix $\textstyle {\begin{bmatrix}1&0&x\\0&1&y\end{bmatrix}},$

  and perform the Euclidean algorithm simultaneously on the pairs (x + A, y + C) and (x + B, y + D), until the quotients differ. That is, iterate as an inner loop:

  • Compute the quotients w₁ of the long divisions of (x + A) by (y + C) and w₂ of (x + B) by (y + D) respectively. If these are equal, they also equal w, the (not computed) quotient from the corresponding step of the ordinary Euclidean algorithm.

  • • If w₁ ≠ w₂, then break out of the inner iteration. Else set w to w₁ (or w₂).
    • Replace the current matrix

    $\textstyle {\begin{bmatrix}A&B&x\\C&D&y\end{bmatrix}}$

    with the matrix product

    $\textstyle {\begin{bmatrix}0&1\\1&-w\end{bmatrix}}\cdot {\begin{bmatrix}A&B&x\\C&D&y\end{bmatrix}}={\begin{bmatrix}C&D&y\\A-wC&B-wD&x-wy\end{bmatrix}}$

    according to the matrix formulation of the extended Euclidean algorithm.

    • If B ≠ 0, go to the start of the inner loop.
  • If B = 0, we have reached a deadlock; perform a normal step of the Euclidean algorithm with a and b, and restart the outer loop.
  • Set a to aA + bB and b to Ca + Db (again simultaneously). This applies the steps of the Euclidean algorithm which were performed on the leading digits in compressed form to the long integers a and b. If b ≠ 0 go to the start of the outer loop.

References

• Kapil Paranjape, Lehmer's Algorithm