# Lehmer's algorithm for finding the greatest common divisor

I'm trying to understand Lehmer's GCD algorithm but I found two conflicting descriptions about when to end the reduction step.

One comes from Lehmer, and is what Knuth, Jebelean and Wikipedia describe, i.e. to perform the Euclidean algorithm "until the quotients differ".

The other is found in Modern Computer Arithmetic by Brent and Zimmermann: "HalfBezout takes as input two 2-word integers, performs Euclid’s algorithm until the smallest remainder fits in one word".

The problem is that these two methods give slightly different results. For example, for two 64-bits words 10499958131665514997 and 14799178230035213023, the first method computes the following coefficients:

[0, 1, 1, 0]
[1, 0, -1, 1]
[-1, 1, 3, -2]
[3, -2, -7, 5]
[-7, 5, 24, -17]
[24, -17, -31, 22]
[-31, 22, 148, -105]
[148, -105, -179, 127]
[-179, 127, 4623, -3280]
[4623, -3280, -4802, 3407]
[-4802, 3407, 9425, -6687]
[9425, -6687, -23652, 16781]
[-23652, 16781, 80381, -57030]
[80381, -57030, -104033, 73811]
[-104033, 73811, 496513, -352274]
[496513, -352274, -600546, 426085]
[-600546, 426085, 2298151, -1630529]
[2298151, -1630529, -2898697, 2056614]
[-2898697, 2056614, 10994242, -7800371]
[10994242, -7800371, -24887181, 17657356]
[-24887181, 17657356, 85655785, -60772439]
[85655785, -60772439, -110542966, 78429795]
[-110542966, 78429795, 196198751, -139202234]
[196198751, -139202234, -306741717, 217632029]
[-306741717, 217632029, 1116423902, -792098321]


while the second method computes one addtional step:

[1116423902, -792098321, -2539589521, 1801828671]


So my question is, which method is the "correct" one? Or perhaps more precisely, does the second method always produces correct reduction, i.e. the first method terminates a bit too early? Is the first method too "cautious"?

This is the Python code I used to generate the sequences with:

def first(x, y):
A, B, C, D = 1, 0, 0, 1
while True:
if y+C == 0 or y+D == 0:
break
q = (x+A)//(y+C)
if q != (x+B)//(y+D):
break
A, B, x, C, D, y = C, D, y, A-q*C, B-q*D, x-q*y
print(map(int, [A, B, C, D]))
return A, B, C, D

def second(x, y):
s = 0; old_s = 1
t = 1; old_t = 0
r = y; old_r = x
while r >= 2**32:
quotient = old_r//r
old_r, r = r, old_r - quotient*r
old_s, s = s, old_s - quotient*s
old_t, t = t, old_t - quotient*t
print(map(int, [old_s, old_t, s, t]))
return old_s, old_t, s, t

x, y = 10499958131665514997, 14799178230035213023
first(x, y)
second(x, y)

• I suspect that both algorithms are correct, meaning that both find solutions to $ax + by = \text{gcd}$ (that equation has many solutions). Have you checked that? If not there's a serious problem you should investigate, If so, use whichever algorithm you prefer. – Ethan Bolker Nov 11 '17 at 15:24
• @EthanBolker Yes, from what I tested both reductions compute correct GCD. But that's part of my question, whether I just got lucky or not... – Faaf Nov 11 '17 at 15:28
• The odds that luck played a part is miniscule given the fact that the algorithm is easy and well understood and the size of your numbers. I don' think you have anything to worry about. Run a few more test cases if you need more reassurance. – Ethan Bolker Nov 11 '17 at 15:32
• @MishaLavrov Please see the update. – Faaf Nov 11 '17 at 17:41
• @MishaLavrov Not sure what that implies, the idea is not to work with 32-bit, instead of 64-bit, the saving wouldn't probably be work it. It's about computing with top bits of numbers with much greater bit length. – Faaf Nov 11 '17 at 18:13

As far as validity of the corresponding algorithm goes, we have considerable freedom. At each step, we have two pairs of coefficients $(A,B)$ and $(C,D)$ that are:

• Computed from $x$ and $y$, the leading bits of long integer inputs $X$ and $Y$.
• Intended to be useful for $X$ and $Y$: we want $AX+BY$ and $CX+DY$ to be significantly smaller than $X$ and $Y$.
• Valid GCD simplifications for any inputs: we have $\gcd(Ax'+By',Cx'+Dy') = \gcd(x',y')$ for any $x'$ and $y'$, because we are multiplying by an integer-invertible matrix at each step, so we can express $x'$ and $y'$ as integer linear combinations of $Ax'+By'$ and $Cx'+Dy'$.

In both of the slightly-different implementations of this algorithm, we use $x$ and $y$ to compute pairs of coefficients until some stopping point. Then we replace $X$ and $Y$ with $AX+BY$ and $CX+DY$ and repeat the algorithm. If we stop sooner or later than necessary, we can't really go wrong; we will always compute the correct GCD in the end. So what happens if we change the stopping point?

If we put the stopping point too soon, then we won't have made as much progress, so $AX+BY$ and $CX+DY$ will be larger than necessary. We will have to do more iterations of the outer loop.

If we put the stopping point too late, then the progress we make on $\gcd(x,y)$ will have stopped reflecting progress on $\gcd(X,Y)$. The values $Ax+By$ and $Cx+Dy$ will keep shrinking, but $AX+BY$ and $CX+DY$ will behave essentially randomly: they may shrink or they may start growing again. We may have to do more iterations of the outer loop.

Lehmer's algorithm as outlined in Wikipedia puts the stopping point as exactly the point at which optimal steps for $(x,y)$ are no longer necessarily optimal for $(X,Y)$. We know that $\frac{AX+BY}{CX+DY}$ and $\frac{Ax+By}{Cx+Dy}$ are both always between $q_1 = \frac{A(x+1)+By}{C(x+1)+Dy}$ and $q_2 = \frac{Ax+B(y+1)}{Cx+D(y+1)}$, so if there is not an integer between $q_1$ and $q_2$, the calculation for $(X,Y)$ is unambiguously the same as the calculation for $(x,y)$. Once that stops being true, the algorithm stops.

The HalfBezout algorithm in Modern Computer Arithmetic doesn't bother calculating the upper and lower bounds. Instead, it knows that the stopping point will come approximately when the remainders are half the size of a word, so it just stops when that condition is satisfied.

This heuristic is sometimes too careful and sometimes too reckless. In particular, when $x$ or $y$ are smaller than expected, it will actually perform fewer steps than it could. In other cases, such as your example, it performs more steps than the quotient-comparison method.

However, it only needs to compute one quotient at each step rather than two, and it's a heuristic that guarantees a rate of progress (by approximately one half-word) in $\gcd(X,Y)$, so it's probably fine.