Lehmer's GCD algorithm

Algorithm

• If b contains only one digit, use some other method, such as the Euclidean algorithm, to obtain the result.
• If a and b differ in the length of digits, perform a division so that a and b are equal in length, with length equal to m.
• Outer loop: Iterate until one of a or b is zero:
• * Decrease m by one. Let x be the leading digit in a, x = a div β ^m and y the leading digit in b, y = b div β ^m.
• * Initialize a 2 by 3 matrix
• :: to an extended identity matrix
• :and perform the euclidean algorithm simultaneously on the pairs and, until the quotients differ. That is, iterate as an inner loop:
• :* Compute the quotients w₁ of the long divisions of by and w₂ of by respectively. Also let w be the quotient from the current long division in the chain of long divisions of the euclidean algorithm.
• ** If w₁ ≠ w₂, then break out of the inner iteration. Else set w to w₁.
• ** Replace the current matrix
• *::
• *: with the matrix product
• *::
• *: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. This applies the steps of the euclidean algorithm that 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.