title: Lehmer's extended GCD algorithm implementation
body: <p>While doing my own BigInteger implementation, I got stuck with the extended GCD algorithm, which is fundamental for finding modular multiplicative inverse. As the well-known Euclidean approach performs too slow, with hybrid and binary algorithms only 5­-10 times faster, the choice was for the Lehmer's modification to the classic algorithm. But the difficulty is that, when it comes to describing the Lehmer's, all books that I found (Knuth, Handbook of Applied Cryptography, Internets, etc) have the same shortcomings:</p>
  
  <ol>
  <li>Explanation is based on several tricks:
  <ul>
  <li>the input numbers are always of the same length;</li>
  <li>the abstract CPU has <em>signed</em> registers, which can hold both the digit and the sign;</li>
  <li>the abstract CPU has <em>semi-unlimited</em> registers, i. e. it never overflows.</li>
  </ul></li>
  <li>Only the basic GCD algorithm is provided, without focusing on the inverse cofactors.</li>
  </ol>
  
  <p>As for the first problem, I was initially surprised by being unable to find any real-world implementation (don't point me to the GNU MP library — it's not a source to learn from), but finally took inspiration by decompiling the Microsoft's implementation from .Net 4.0, which is obviously based on the ideas from the paper “<a href="http://www.csie.nuk.edu.tw/~cychen/gcd/A%20double-digit%20Lehmer-Euclid%20algorithm%20for%20finding%20the%20GCD%20of%20long%20integers.pdf" rel="nofollow">A double-digit Lehmer-Euclid algorithm for finding the GCD of long integers</a>” by Jebelean. The resulting function is large, it looks scary, but works just great.</p>
  
  <p>But Microsoft's library provides the basic functionality only, no cofactors are computed. Well, to be precise, some cofactors <em>are</em> computed during the shorthand step, and during the very first step those cofactors simply <em>are</em> the initial ones, but after the longhand step is performed then they do not match anymore. My current solution is to update the “real” cofactors in parallel with the “substitute” ones (except the very first step), but it makes the performance to drop below zero: the function now completes only 25­-50 % faster than the binary method in basic mode. So, the problem is that, <strong>while the input numbers are fully updated during longhand steps only, the cofactors are updated on each shorthand step's iteration as well</strong>, thus destroying almost any benefit from Lehmer's approach.</p>
  
  <p>To speed up things a little, I implemented a “fused multiply-add” function, because a “fused multiply-multiply-subtract” really does help updating the input numbers, — but this time the impact was negligible. Another improvement is based on the fact that usually only one cofactor is necessary, so the other one can be just not computed at all. This should halve the overhead (or even more so, since the second number is usually significantly smaller than the first one), however in practice the overhead reduces only by 25 to 50 % <em>of expected</em>.</p>
  
  <p>Consequently, my questions come down to this:</p>
  
  <ol>
  <li>Is there any full-scale explanation of Lehmer's algorithm, tied to practical implementation on real-world hardware (with <em>unsigned</em> words of <em>limited</em> size)?</li>
  <li>Same as above, but regarding the <em>extended</em> GCD computation.</li>
  </ol>
  
  <p>So, as much as I'm happy with the performance of basic algorithm, the opposite applies to the extended mode of operation, which is the primary in my case.</p>