Stein算法 - 求最大公约数

欧几里德算法（辗转相除法）是计算两个数最大公约数的传统算法，它无论从理论还是从效率上都是很好的。但是他有一个致命的缺陷，这个缺陷只有在大素数时才会显现出来。

Stein算法由 J. Stein 1961年提出，和欧几里德算法算法不同的是，Stein算法只有整数的移位和加减法。

为了说明Stein算法的正确性，首先必须注意到以下结论：

gcd(a,a) = a，也就是一个数和他自身的公约数是其自身
gcd(ka,kb) = k gcd(a,b)，也就是最大公约数运算和倍乘运算可以交换，特殊的，当k=2时，说明两个偶数的最大公约数必然能被2整除。

有了上述规律就可以给出Stein算法如下：

1.如果A=0，B是最大公约数，算法结束
2.如果B=0，A是最大公约数，算法结束
3.设置A1 = A、B1=B和C1 = 1
4.如果An和Bn都是偶数，则An+1 =An /2，Bn+1 =Bn /2，Cn+1 =Cn *2(注意，乘2只要把整数左移一位即可，除2只要把整数右移一位即可)
5.如果An是偶数，Bn不是偶数，则An+1 =An /2，Bn+1 =Bn ，Cn+1 =Cn (很显然啦，2不是奇数的约数)
6.如果Bn是偶数，An不是偶数，则Bn+1 =Bn /2，An+1 =An ，Cn+1 =Cn (很显然啦，2不是奇数的约数)
7.如果An和Bn都不是偶数，则An+1 =|An -Bn|，Bn+1 =min(An,Bn)，Cn+1 =Cn
8.n++，转4

C++实现：

unsigned int gcd(unsigned int a,unsigned int b)
{
    if(a<b) swap(a,b);
    if(b==0) return a;
    if(a%2==0&&b%2==0) return 2*gcd(a/2,b/2);
    if(a%2==0)    return gcd(a/2,b);
    if(b%2==0)    return gcd(a,b/2);
    return gcd((a+b)/2,(a-b)/2);
}

ANSI C:

来自http://blog.csdn.net/ztj111/archive/2007/11/27/1905015.aspx

unsigned int stein( unsigned int x, unsigned int y )
/* return the greatest common divisor of x and y */
{
    unsigned int factor = 0;
    unsigned int temp;
    if ( x < y ){
        temp = x;
        x = y;
        y = temp;
    }
    if ( 0 == y )    return 0;
    while ( x != y )
    {
        if ( x & 0x1 )
        {/* when x is odd */
            if ( y & 0x1 )
            {/* when x and y are both odd */
                y = ( x - y ) >> 1;
                x -= y;
            }
            else
            {/* when x is odd and y is even */
                y >>= 1;
            }
        }
        else
        {/* when x is even */
            if ( y & 0x1 )
            {/* when x is even and y is odd */
                x >>= 1;
                if ( x < y )
                {
                    temp = x;
                    x = y;
                    y = temp;
                }
            }
            else
            {/* when x and y are both even */
                x >>= 1;
                y >>= 1;
                ++factor;
            }
        }
    }
    return ( x << factor );
}