Lehmer's GCD algorithm: Difference between revisions

Revision as of 00:31, 23 December 2006

Latest revision as of 18:19, 11 January 2020

(32 intermediate revisions by 20 users not shown)

Line 1:

'''Lehmer's{{short GCDdescription|Fast algorithm''' is a rather fast [[greatest common divisor|GCD]] [[algorithm]], an improvement on the simpler but slower [[Euclidean algorithm]].}}

'''Lehmer's GCD algorithm''', named after [[Derrick Henry Lehmer]], is a fast [[greatest common divisor|GCD]] [[algorithm]], an improvement on the simpler but slower [[Euclidean algorithm]]. It is mainly used for big integers that have a representation as a string of digits relative to some chosen numeral system [[radix|base]], say ''β'' = 1000 or ''β'' = 2³².

== Algorithm ==

[[Derrick 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 name="knuth2">[[Donald Knuth|0}}Knuth]], ''[[The Art of Computer Programming]] vol 2 "Seminumerical algorithms"'', chapter 4.5.3 Theorem E.</ref>) 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 <math>\textstyle ''a'' \&ge; ''b</math>''.

* If ''b'' contains only one digit (in the chosen [[radix|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, perform a division so that ''a'' and ''b'' are equal in length, with length equal to ''m''.

* Let ''xOuter loop:'' beIterate theuntil leadingone (most significant) digit inof ''a'' andor ''yb'' theis leading digit in ''b''.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 square [[matrix (mathematics)|matrix]] <math>\textstyle \begin{Bmatrix} A & B \\ C & D \end{Bmatrix}</math> to the [[identity matrix]] <math>\textstyle \begin{Bmatrix} 1 & 0 \\ 0 & 1 \end{Bmatrix}</math>, and iterate:

** Initialize a 2 by 3 [[matrix (mathematics)|matrix]]

** Compute the quotients ''w₁'' of (''x'' + ''A'')(''y'' + ''C'') and ''w₂'' of (''x'' + ''B'')(''y'' + ''D'') respectively. Also let ''w'' be the quotient from ''a''/''b''.

*::<math>\textstyle

** If ''w₁'' = ''w₂'', set ''w'' to ''w₁'' (or ''w₂'').

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

** Set our current matrix <math>\textstyle \begin{Bmatrix} A & B \\ C & D \end{Bmatrix}</math> to <math>\textstyle \begin{Bmatrix} 0 & 1 \\ 1 & -w \end{Bmatrix} \bullet \begin{Bmatrix} A & B \\ C & D \end{Bmatrix} = \begin{Bmatrix} C & D \\ A - wC & B - wD \end{Bmatrix}</math>

</math> to an extended [[identity matrix]] <math>\textstyle

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

\begin{bmatrix} 1 & 0 & x \\ 0 & 1 & y\end{bmatrix},

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

</math>

*: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. 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''₂'').

*** Replace the current matrix

**::<math>\textstyle \begin{bmatrix} A & B & x \\ C & D & y \end{bmatrix}</math>

**: with the matrix product

**:: <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>

**: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 divisionstep with ''a'' and ''b'', and recomputeof the matrix.euclidean Otherwise,algorithm setwith ''a'' to ''aA'' + ''bB'' and ''b'', toand ''Ca''restart +the ''Db'' (again simultaneously), andouter iterateloop.

** 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''. If ''b'' ≠ 0 go to the start of the outer loop.

== References ==

<references/>

* {{note|0}} Knuth, ''[[The Art of Computer Programming]]'' volume 4, section 5.3, Theorem E

* Kapil Paranjape, [http://www.imsc.res.in/~kapil/crypto/notes/node11.html Lehmer's Algorithm]

* [[Kapil Hari Paranjape|Kapil Paranjape]], [http://www.imsc.res.in/~kapil/crypto/notes/node11.html Lehmer's Algorithm]

{{numtheory-stub}}

{{Number-theoretic algorithms}}

{{DEFAULTSORT:Lehmer's Gcd Algorithm}}

[[Category:Number theoretic algorithms]]