I just found this very cool GCD algorithm by Stehle & Zimmerman; it is built on a *** backwards *** division algorithm (!), more reminiscent of something you might want to do with 2-adic ('backwards') integers. Divide (ahi,alo) by b to get q, r. But instead of getting (ahi,alo) - q*b to get r<b, we want q', r', s.t. (ahi,alo) - q'*b = (r',0), where r'<b. (Here, '*' & '-' have their usual meaning.) I.e., instead of q killing off the _high order_ digit ahi, we us q' to kill off the _low order_ digit alo !! http://perso.ens-lyon.fr/damien.stehle/downloads/recbinary.pdf Stehle, Damien & Zimmerman, Paul. "A Binary Recursive GCD Algorithm". ~2004. So when calculating gcd(m,n) we can use a gcd on the _low order_ parts of m & n to kill off these low order bits, and then work on the high order part. Henry Baker At 04:58 PM 3/25/2013, Warren D Smith wrote:
From: Henry Baker <hbaker1@pipeline.com> I'm looking for the _simplest_ binary _recursive_ GCD algorithm that has "pretty good" performance.
I.e., a GCD algorithm that takes two length 2m integers and reduces this to a number of computations on length m integers.
--Let the two input numbers be X and Y with X<Y. X=(Xhi,Xlo) and Y=(Yhi,Ylo) for the top and bottom halves, same lengths for all four. We assume Xhi>0 otherwise you can immediately replace Y by Y mod X to reduce to half precision.
Algorithm idea: Use extended GCD to find a,b and g so that a*Xhi+b*Yhi=g=gcd(Xhi,Yhi). Replace Y by a*X+b*Y. Replace X by X mod Y. Replace Y by Y mod X. swap(X,Y). Now continue on.
Seems to me, each iteration will reduce Yhi to about its square root (or less) thus multiplying the bitlength by 0.75 (or less).
I've cheated a bit by using extended GCD as a subroutine inside plain GCD. However, you can add extra bells & whistles to the plain gcd algorithm to keep track of what happens to the coeffs, thus converting it to an extended GCD algorithm...