Stein算法求最大公约数 ( ANSI C )

首先引进一个符号：g_c_d是greatest common divisor（最大公约数）的缩写，g_c_d( x,y ) 表示x和y的最大公约数。然后有一个事实需要了解：一个奇数的所有约数都是奇数。这个很容易，下面我们要用到。
    来研究一下最大公约数的性质，我们发现有 g_c_d( k*x,k*y ) = k*g_c_d( x,y ) 这么一个非常好的性质（证明我就省去了）。说他好是因为他非常符合我们化小的思想。我们试取 k=2 ，则有 g_c_d( 2x,2y ) = 2 * g_c_d( x,y )。这使我们很快联想到将两个偶数化小的方法。那么一奇一个偶以及两个奇数的情况我们如何化小呢？

    先来看看一奇一偶的情况： 设有2x和y两个数，其中y为奇数。因为y的所有约数都是奇数，所以 a = g_c_d( 2x,y ) 是奇数。根据2x是个偶数不难联想到，a应该是x的约数。我们来证明一下：(2x)%a=0，设2x=n*a，因为a是奇数，2x是偶数，则必有n是偶数。又因为 x=(n/2)*a，所以 x%a=0，即a是x的约数。因为a也是y的约数，所以a是x和y的公约数，有 g_c_d( 2x,y ) <= g_c_d( x,y )。因为g_c_d( x,y )明显是2x和y的公约数，又有g_c_d( x,y ) <= g_c_d( 2x,y )，所以 g_c_d( 2x,y ) = g_c_d( x,y )。至此，我们得出了一奇一偶时化小的方法。

    再来看看两个奇数的情况：设有两个奇数x和y，似乎x和y直接向小转化没有什么太好的办法，我们可以绕个道，把x和y向偶数靠拢去化小。不妨设x>y，我们注意到x+y和x-y是两个偶数，则有 g_c_d( x+y,x-y ) = 2 * g_c_d( (x+y)/2,(x-y)/2 )，那么 g_c_d( x,y ) 与 g_c_d( x+y,x-y ) 以及 g_c_d( (x+y)/2,(x-y)/2 ) 之间是不是有某种联系呢？为了方便我设 m=(x+y)/2 ，n=(x-y)/2 ，容易发现 m+n=x ，m-n=y 。设 a = g_c_d( m,n )，则 m%a=0,n%a=0 ，所以 (m+n)%a=0，(m-n)%a=0 ，即 x%a=0 ，y%a=0 ，所以a是x和y的公约数，有 g_c_d( m,n )<= g_c_d(x,y)。再设 b = g_c_d( x,y )肯定为奇数，则 x%b=0,y%b=0 ，所以 (x+y)%b=0 ，(x-y)%b=0 ，又因为x+y和x-y都是偶数，跟前面一奇一偶时证明a是x的约数的方法相同，有 ((x+y)/2)%b=0,((x-y)/2)%b=0 ，即 m%b=0 ，n%b=0 ，所以b是m和n的公约数，有 g_c_d( x,y ) <= g_c_d( m,n )。所以 g_c_d( x,y ) = g_c_d( m,n ) = g_c_d( (x+y)/2,(x-y)/2 )。

我们来整理一下，对两个正整数 x>y ：
1.均为偶数 g_c_d( x,y ) =2g_c_d( x/2,y/2 )；
2.均为奇数 g_c_d( x,y ) = g_c_d( (x+y)/2,(x-y)/2 )；
2.x奇y偶   g_c_d( x,y ) = g_c_d( x,y/2 )；
3.x偶y奇   g_c_d( x,y ) = g_c_d( x/2,y )  或 g_c_d( x,y )=g_c_d( y,x/2 )；

现在我们已经有了递归式，还需要再找出一个退化情况。注意到 g_c_d( x,x ) = x ，我们就用这个。

int Stein( unsigned int x, unsigned int y )
/**//* return the greatest common divisor of x and y */
...{
        int factor = 0;
        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 );
}

 

Stein算法的优势是不言而喻的，上面的代码只有位运算和减法，不仅速度加快，而且结局了欧几里德算法求两个大数最大公约数的不便。不仅如此，Stein算法还可以很容易的写成汇编语言实现。