求最大公约数算法

一、欧几里德算法（辗转相除法）

定理：对于任意整数a,b,有gcd(a,b) = gcd(b,a mod b)

证明：

a = kb + r --> r = a % b

假设a、b的公约数为x，即 x|a、x|b  ("|"表示a、b可被x整除)，用同余符号表示为：a≡b(mod x)。

x|b --> x|kb

由取模运算性质

可得：x|(a - b) --> x|(a - kb) --> x|(kb+r - kb) --> x|(r)

即r可以被x整除，所以x是a，b，r的公约数，也就是说a，b，r有相同的公约数，那么最大公约数自然也相等。

定理证毕。

 

C代码实现：

1.递归实现

int gcd(int a, int b)
{
    if (b == 0)
        return a;

    return gcd(b, a%b);
}

 

2.迭代实现

int gcd(int a, int b)
{
    int tmp;

    while (b != 0)
    {
        tmp = a;
        a = b;
        b = tmp % a;
    }

    return a;
}

 

二、拓展欧几里德算法

扩展欧几里得算法是欧几里得算法（又叫辗转相除法）的扩展。除了计算a、b两个整数的最大公约数，此算法还能找到整数x、y（其中一个很可能是负数），使它们满足贝祖等式

 注意x、y的解不唯一。

原理：

假设a' = b,b' = a%b.

内层中，假设存在x,y,使得a'x + b'y = gcd(a',b') --> bx + (a%b)*y = gcd(a',b') --> bx + (a - a/b * b)*y = gcd(a',b')

在欧几里德算法中已经证明gcd(a,b) = gcd(b,a%b) = gcd(a',b')

所以：bx + (a - a/b * b)*y = gcd(a,b) --> ay + b(x - y*a/b) = gcd(a,b).

显然：只要求出内层a',b'对应的x,y,就能够求出外层a,b对应的x,y.

 

C代码实现：

int gcd(int a, int b, int *x, int *y)
{
    int ret;
    int tmp;

    if (b == 0)
    {
        *x = 1;
        *y = 0;
        return a;
    }
    
    ret = gcd(b, a%b, x, y);

    tmp = *x;
    *x = *y;
    *y = tmp - (*y)*(a/b);

    return ret;
}

 

从代码中可以看出，递归极限（b=0）的情况下，a*x + b*y = gcd(a,0) = a,这样就很容找出一组x，y满足此方程。

我们选择x=1，y=0。

 

三、Stein算法

当求两个非常大的数的最大公约数且其中有素数时，欧几里德算法在效率上有缺陷，Stein算法可以弥补此缺陷。

算法规则如下：

x=0,则y是最大公约数。

y=0,则x是最大公约数。

x偶y偶：gcd(x,y) = 2gcd(x/2,y/2)

x偶y奇：gcd(x,y) = gcd(x/2,y)

x奇y偶：gcd(x,y) = gcd(x,y/2)

x奇y奇：gcd(x,y) = gcd((x+y)/2,|x-y|/2)

基本思想就是把要求的两个数不断化小，直到其中一个为0.

 

为了证明上述等式成立，需要先了解几个预备知识：

1.一个奇数的所有约数都是奇数

2.gcd(nx,ny) = n * gcd(x,y)

3.如果x%a=0,y%a=0,则(x+y)%a=0,(x-y)%a=0　　（这个很容易理解）

4.奇数 * 偶数 = 偶数　　奇数 + 奇数 = 偶数　　奇数 - 奇数 = 偶数

 

下面开始证明：

1.x偶y偶：gcd(x,y) = 2gcd(x/2,y/2)

这个很简单，根据gcd(nx,ny) = n * gcd(x,y)，令n=2即可的上述等式。

 

2.x偶y奇：gcd(x,y) = gcd(x/2,y)

假设a=gcd(x,y),根据预备知识1可知，a是奇数。x%a=0,设x=na,根据预备知识4可知，n必定为偶数，

所以x/2=(n/2)*a,也就是x/2能够被a整除，(x/2)%a=0.

x/2,y有公约数a --> gcd(x,y) <= gcd(x/2,y)　　(这里应该很容易理解吧)

gcd(x,y) >= gcd(x/2,y)　　(显而易见)

所以，gcd(x,y) = gcd(x/2,y)，证毕。

 

3.x奇y偶：gcd(x,y) = gcd(x,y/2)

和上面的2相同，不重复。

 

4.x奇y奇：gcd(x,y) = gcd((x+y)/2,|x-y|/2)

当x>y时，

根据预备知识4可知，(x+y),(x-y)都是偶数，所以gcd((x+y),(x-y)) = 2*gcd((x+y)/2,(x-y)/2).

假设(x+y)/2=m,(x-y)/2=n,那么m+n=x,m-n=y.

如果a=gcd(m,n) --> m%a=0,n%a=0，根据预备知识3，(m+n)%a=0,(m-n)%a=0 --> x%a=0,y%a=0.

所以a是x,y的公约数 --> gcd(m,n) <= gcd(x,y)

现在来证gcd(m,n) >= gcd(x,y)：

设b=gcd(x,y).因为x,y都是奇数，所以b必定为奇数。

x%b=0,y%b=0 --> (x+y)%b=0,(x-y)%b=0.根据预备知识4可知，(x+y),(x-y)都是偶数，那么((x+y)/2)%b=0,((x-y)/2)%b=0.(这一点在证明2时已经证明过了)

即m%b=0,n%b=0 --> b是m,n的公约数 --> gcd(m,n) >= gcd(x,y)

所以，gcd(x,y) = gcd(m,n) = gcd((x+y)/2,(x-y)/2)

当y>x时，gcd(x,y) = gcd((x+y)/2,(y-x)/2)

这就是为什么要加绝对值符号的原因。

 

int stein(int x, int y)
{
    if (x == 0)
        return y;
    if (y == 0)
        return x;

    /* 右移1位代替除以2 */
    if ((x%2 ==0) && (y%2 == 0))
        return 2 * stein(x>>1, y>>1);

    else if ((x%2 ==0) && (y%2 != 0))
        return stein(x>>1, y);

    else if ((x%2 !=0) && (y%2 == 0))
        return stein(x, y>>1);

    else
        return stein((x+y)>>1, abs(x-y)>>1);
}

 

转载于:https://www.cnblogs.com/oyld/p/3380794.html