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

欧几里得算法（辗转相除法）

求两个正整数的最大公约数，时间复杂度 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))。