扩展欧几里得算法详解



对于不完全为 0的非负整数a，b，gcd(a,b)表示a，b的最大公约数，必然存在整数x和y，使得gcd(a,b)=ax+by。

解析：

设a、b不全为0，令a>b，

当b=0时，gcd(a,b)=a,解的情况为x=1，y=0

当ab！=0，令a*x1+b*y1 = gcd(a,b)，所以b*x2+(a%b)*y2 =gcd(b,a%b)

又gcd(a,b) = gcd(b,a%b)（欧几里得算法）

故有a*x1+b*y1 = b*x2+(a%b)*y2

 = b*x2+( a-(a/b)*b )*y2
              = b*x2+a*y2-(a/b)*b*y2
              = a*y2+b*(x2-(a/b)*y2)
即有x1=y2, y1=x2-(a/b)*y2
因此x1与y1的值可由x2、y2推知，拓展欧几里得算法的求解过程就是不断地的将b放小，直至b等于0，最后反推求x和y。

//Date:2015.05.05
//扩展欧几里得算法的递归与非递归实现

#include <ctime>
#include <random>
#include <cstdio>
using namespace std;

int recursive_gcd(int a,int b,int &x,int &y){
    if(b==0){
        x=1;
        y=0;
        return a;
    }
    int r=recursive_gcd(b,a%b,x,y);
    int t=x;
    x=y;
    y=t-a/b*y;
    return r;
}

int iterative_gcd(int m,int n,int &x,int &y){
    int x1,y1,x0,y0;
    x0=1;x1=0;x=0;
    y0=0;y1=1;y=1;

    int r=m%n;
    int q=m/n;
    while(r){
        x=x0-q*x1;x0=x1;x1=x;
        y=y0-q*y1;y0=y1;y1=y;
        m=n; n=r; r=m%n;
        q=m/n;
    }
    return n;
}

int main(){
    int a,b,x,y,r;
    mt19937 engine(time(NULL));
    uniform_int_distribution<int> distribution(0,100000);
    a=distribution(engine);b=distribution(engine);

    printf("a*x+b*y=gcd(a,b)\n");
    r = recursive_gcd(a,b,x,y);
    printf("递归 : (%d) * (%d) + (%d) * (%d) = %d\n",a,x,b,y,r);
    r = iterative_gcd(a,b,x,y);
    printf("迭代 : (%d) * (%d) + (%d) * (%d) = %d\n",a,x,b,y,r);
    return 0;
}

应用：

1.求解不定方程ax+by=c :

对于不定整数方程ax+by=c ，若cmod gcd(a,b)=0,则该方程存在整数解，否则不存在整数解。 

故当c mod gcd(a,b)=0时，先用扩展欧几里得算法求出ax+by=gcd(a,b)的一组解x1,y1。

令k= c/gcd(a,b)；则x2=x1*k，y2=y1*k为ax+by=c的一组解。

x=x2+b/gcd(a,b)*t，y=y2-a/gcd(a,b)*t (t为整数)，即为ax+by=c的所有解。

 

2.求解乘法逆元

在密码学中，有限域GF(p)是一个很重要的域，其中p为素数。简单来说，GF(p)就是mod p，因为一个数模p后，结果在[0, p-1]之间。GF(p)里面的乘法即为一般的乘法运算。为什么p一定要是一个素数呢？这是因为当p为素数时，才能保证集合中的所有的元素都有乘法逆元(0除外)。假设p等于10，对于元素2，找不到一个数x，使得2*x mod10等于1。如果p是素数，那么对于域中的任一个元素a，总能在域中找到另外一个元素x，使得a*xmod p等于1，其中a与互为乘法逆元。也即若gcd(a,p)=1；则ax+py=1；

 

但是我们希望0到255这256个数字也能组成一个域，因为很多领域需要用到。mod 256的余数范围是0到255，不过256不是素数。为此我们定义一个素多项式m(x)=x^8+ x^4 + x^3 +x +1（类似素数，素多项式不能表示为其他两个多项式的乘积），于是0到255可以通过mod m(x)这样的方式组成一个域GF(2^8)，域中的每一个元素都存在乘法逆元(0的乘法逆元定义为0)。

注：有限域GF(2^8)中的加法运算为异或运算，减法与加法等价。

//Date:2015.05.05
//扩展欧几里得算法求解有限域GF(2^8)内的乘法逆元

#include<cstdio>
using namespace std;
const int Bit_Num=sizeof(int)*8;

//求解非零最高位
int index_of_max(int value){
    int index=0;
    for(int i=0;i<Bit_Num;++i)
        if(value & (1<<i))
            index=i;
    return index;
}

//有限域GF(2^8)内的除法运算
int divide(int m,int b,int &r){
    int m_MSB=index_of_max(m);
    int b_MSB=index_of_max(b);
    if(m_MSB<b_MSB){
        r=m;
        return 0;
    }
    int d=m_MSB-b_MSB;
    int temp=b;
    temp=temp<<d;
    m=m^temp;
    return (1<<d)|divide(m,b,r);
}

//迭代求解x，y
int iterate(int d0,int q,int d1){
    int value=0;
    for(int i=0;i<Bit_Num;++i){
          if(q & (1<<i))
             value=value^((d1<<i));
    }
    return d0^(value);
}

int exgcd(int m,int a,int &x,int &y){
    int x0,x1,y0,y1,q,r=0;
    x0=0;x1=1;
    y0=1;y1=0;

    while(1){
        if(a==0)
            return m;
        if(a==1)
            return a;
        q=divide(m,a,r);
        x=iterate(x0,q,x1);x0=x1;x1=x;
        y=iterate(y0,q,y1);y0=y1;y1=y;
        m=a;a=r;
    }
}

int main(void){
    int m=283,a,x,y;    //m(x) = x^8 + x^4 + x^3 + x + 1;

    while(1){
        printf("请输入:");
        scanf("%X",&a);
        if(a==0)    break;
        exgcd(m,a,x,y);
        printf("%X * %X + %X * %X = 1\n",a,x,m,y);
        printf("%X的乘法逆元是%X\n\n",a,x);
    }
    return 0;
}