Lehmer's GCD algorithm: Difference between revisions

Revision as of 18:17, 16 August 2008

Revision as of 18:20, 16 August 2008

Line 1:

'''Lehmer's GCD algorithm''', named after [[Derrick Henry Lehmer]], is a rather fast [[greatest common divisor|GCD]] [[algorithm]], an improvement on the simpler but slower [[Euclidean algorithm]].

== Algorithm ==

[[Derrick Henry Lehmer|Lehmer]] noted that that most of the [[quotient]]s from each step of the division part of the standard algorithm are small. (For example, [[Donald Knuth|Knuth]] observed that the quotients 1, 2, and 3 comprise 67.7% of all quotients.)

Say we want to obtain the GCD of the two integers ''a'' and ''b''. Let <math>\textstyle ''a \''&nbsp;&ge ;&nbsp;''b</math>''.

* If ''b'' contains only one digit (in the chosen [[radix|base]], say <math>\beta=1000</math> or <math>\beta=2^{32}</math>), 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''.

Line 10:

** Decrease ''m'' by one. Let ''x'' be the leading (most significant) digit in ''a'', <math>x=a \text{ div } \beta^m</math> and ''y'' the leading digit in ''b'', <math>y=b \text{ div } \beta^m</math>.

** Initialize a 2 by 3 [[matrix (mathematics)|matrix]] <math>\textstyle \begin{bmatrix} A & B & x\\ C & D & y \end{bmatrix}</math> to an extended [[identity matrix]] <math>\textstyle \begin{bmatrix} 1 & 0 & x \\ 0 & 1 & y\end{bmatrix}</math>, and perform the euclidean algorithm simultaneously on the pairs <math>(x+A,y+C)</math> and <math>(x+B,y+D)</math>, until the quotients differ. That is, iterate:

*** Compute the quotients ''w₁'' of the long divisions of (''x'' + ''A'') by (''y'' + ''C'') and ''w''₂'' of (''x'' + ''B'') by (''y'' + ''D'') respectively. Also let ''w'' be the (not computed) quotient from the current long division in the chain of long divisions of the euclidean algorithm.

*** If ''w''₁'' &ne; ''w''₂'', then break out of the inner iteration. Else set ''w'' to ''w''₁'' (or ''w''₂'').

*** Set our current matrix <math>\textstyle \begin{bmatrix} A & B & x \\ C & D & y \end{bmatrix}</math> to <math>\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}</math>

*** If ''B'' = 0, we have reached a ''deadlock''; perform a normal step of the euclidean algorithm with ''a'' and ''b'', and recompute the matrix.