欧几里得算法和扩展欧几里得算法

欧几里得算法(辗转相除法)

求两个正整数的最大公约数,时间复杂度 O(logn)。

欧几里得算法基于下面这个原理:
设a,b均为正整数,则gcd(a,b)=gcd(b,a%b)。

证明:设 a = kb + r,其中 k 和 r 分别为 a 除以 b 得到的商和余数。
则有 r = a - kb 成立。
设 d 为 a 和 b 的一个公约数,
那么由 r = a - kb,得 d 也是 r 的一个约数。
因此 d 是 b 和 r 的一个公约数。
又由 r = a % b,得 d 为 b 和 a % b 的一个公约数。
因此 d 既是 a 和 b 的公约数,也是 b 和 a % b的公约数。
由 d 的任意性,得 a 和 b 的公约数都是 b 和 a % b 的公约数。
由 a = kb + r,同理可证 b 和 a % b的公约数都是 a 和 b 的公约数。
因此 a 和 b 的公约数与 b 和 a % b的公约数都是 a 和 b 的公约数。
即有 gcd(a,b) = gcd(b,a % b)。
证毕。

递归的两个关键:
① 递归式:gcd(a,b) = gcd(b,a % b) 。
② 递归边界:gcd(a,0) = a 。

于是可得到下面的求解最大公约数的代码:

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

更简洁的写法:

int gcd(int a, int b){
    return !b ? a : gcd (b, a % b);
}

对于最小公倍数的实现可以通过求解最大公约数来间接的实现,即a,b的最小公倍数lcm(a,b)=a/gcd(a,b)*b

扩展欧几里得算法

裴蜀定理:若 a,b 是整数,且 gcd(a,b)=d,那么对于任意的整数 x,y, ax+by 都一定是 d 的倍数,特别地,一定存在整数 x,y,使 ax+by=d 成立。

证明:设 a>b。
显然当 b=0,gcd(a,b)=a。此时 x=1,y=0;
当ab!=0 时
设 ax1+by1=gcd(a,b);
bx2+(a mod b)y2=gcd(b,a mod b);
根据朴素的欧几里德原理有 gcd(a,b)=gcd(b,a mod b);
则:ax1+by1=bx2+(a mod b)y2;
即:ax1+by1=bx2+(a-(a/b)b)y2=ay2+bx2-(a/b)by2;
根据恒等定理得:x1=y2; y1=x2-(a/b)*y2;
这样我们就得到了求解 x1,y1 的方法:x1,y1 的值基于 x2,y2.

扩展欧几里得算法可以在 O(logn) 的时间复杂度内求出系数 x, y。

int exgcd(int a, int b, int &x, int &y){
    if (b == 0){
        x = 1;
        y = 0;
        return a;
    }
    int g = exgcd(b, a % b, x, y);
    int temp = x;    //存放x的值
    x = y;    //更新x = y(old)
    y = temp - a / b * y;    //更新y = x(old) - a / b * y(old)
    return g;
}

应用

1. 方程ax+by=c的求解

已知:a * x+b * y=gcd

解:
a * x / gcd + b * y / gcd = 1

a * c * x / gcd + b * c * y / gcd = c

利用当前的x,y构造出新的X,Y
令X = x * c / gcd Y = y * c / gcd
则a * X + b * Y = c
因为此时上面的表达式是新构造出来的
为了使该表达式有解,必须满足 c%gcd==0

又因为 X = x0 + k * b / gcd

令B=b/gcd
则x0=X%B,y0与其类似

此时解出来的X,Y即,一对特解
而x0,y0为最小正整数解。

int cal(int a,int b,int c){
	int x,y;
	int gcd=exgcd(a,b,x,y);
	if(c%gcd) return -1;
	x*=c/gcd;
	b/=gcd;
	if(b<0)  b=-b;
	int ans=x%b;
	while(ans<=0) ans+=b;
	return ans;
}

重点

x*=c/gcd;
b/=gcd;


2. 同余式ax ≡ c(mod m)的求解

  先解释什么是同余式。对于整数a、b、m来说,如果m整除a-b(即(a-b)%m=0),那么就说a与b模m同余,对应的同余式为ax ≡ b(mod m),m称为同余式的模。例如10与13模3同余,10也与1模3同余,它们分别记为10 ≡ 13(mod 3)、10 ≡ 1(mod 3)。显然,每一个整数都各自与[0,m)中唯一的整数同余。
  根据同余式的定义,有(ax-c)%m=0成立,因此存在整数y,使得ax-c=my成立。移项并令y=-y后即得ax+my=c。
模  板上同。

2. 逆元的求解以及(b/a)%m的计算

  先解释下什么是逆元(此处特指乘法逆元)。假设a、b、m是整数,m>1,且有ab ≡ 1(mod m)c成立,那么就说a和b互为模m的逆元,一般也记作a≡(frac{1}{b})(mod m)或b≡(frac{1}{a})(mod m)。通俗的说,如果两个整数的乘积模m后等于1,就称为它们互为m的逆元。
  那么逆元有什么用处呢?当要计算 (b/a)%m 时,通过找到 a 模 m 的逆元 x,就有 (b/a)%m = (b* x)%m = ((b%m)* (x%m))%m 成立。
  由定义知,求 a 模 m 的逆元,就是求解同余式 ax ≡ 1(mod m),并且在实际使用中,一般把 x 的最小正整数解称为 a 模 m 的逆元

  • 如果 gcd(a,m) ≠ 1,那么同余式 ax ≡ 1(mod m) 无解,a 不存在模 m 的逆元
  • 如果 gcd(a,m) ≠ 1,那么同余式 ax ≡ 1(mod m) 在 (0,m) 上有唯一解,可以通过求解ax+my=1得到。注意:由于gcd(a,m)=1,因此ax+my=1=gcd(a,m),直接使用拓展欧几里得算法解出 x 之后就可以用 (x%m + m)%m 得到 (0,m) 范围内的解,也就是所需要的逆元。
int inverse(int a,int b){  //求解a模m的逆元
    int x,y;
    int g=exgcd(a,m,x,y);  //求解ax+my=1
    return (x%m+m)%m;  
}

  另外,如果 m 是素数,且 a 不是 m 的倍数,则还可以直接使用费马小定理来得到逆元。

费马小定理:设 m 是素数,a 是任意整数且 (a ≢ 0(mod m)),则 (a^{m-1} ≡ 1(mod m))

原文地址:https://www.cnblogs.com/transmigration-zhou/p/12302087.html