Lehmer's GCD algorithm: Difference between revisions

Revision as of 11:19, 7 October 2007

Revision as of 11:58, 7 May 2008

Line 2:

== 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.{{ref|0}})

Say we want to obtain the GCD of the two integers ''a'' and ''b''. Let <math>\textstyle a \ge 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''.

* Iterate until one of ''a'' or ''b'' is zero:

** 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 square2 by 3 [[matrix (mathematics)|matrix]] <math>\textstyle \begin{Bmatrixbmatrix} A & B & x\\ C & D & y \end{Bmatrixbmatrix}</math> to thean extended [[identity matrix]] <math>\textstyle \begin{Bmatrixbmatrix} 1 & 0 & x \\ 0 & 1 & y\end{Bmatrixbmatrix}</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 ''a''/''b''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{Bmatrixbmatrix} A & B & x \\ C & D & y \end{Bmatrixbmatrix}</math> to <math>\textstyle \begin{Bmatrixbmatrix} 0 & 1 \\ 1 & -w \end{Bmatrixbmatrix} \bulletcdot \begin{Bmatrixbmatrix} A & B & x \\ C & D & y \end{Bmatrixbmatrix} = \begin{Bmatrixbmatrix} C & D &y \\ A - wC & B - wD & x-wy \end{Bmatrixbmatrix}</math>

** Set ''x'' to ''y'' and ''y'' to ''x'' - ''wy'' (simultaneously).

*** If ''B'' = 0, we have reached a ''deadlock''; perform a normal divisionstep of the euclidean algorithm with ''a'' and ''b'', and recompute the matrix. Otherwise, set ''a'' to ''aA'' + ''bB'' and ''b'' to ''Ca'' + ''Db'' (again simultaneously), and iterate.

** Set ''a'' to ''aA'' + ''bB'' and ''b'' to ''Ca'' + ''Db'' (again simultaneously). 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''. Restart the outer iteration.

== References ==