先介绍什么叫做欧几里德算法
有两个数 a b,现在,我们要求 a b 的最大公约数,怎么求?枚举他们的因子?不现实,当 a b 很大的时候,枚举显得那么的naïve ,那怎么做?
欧几里德有个十分又用的定理: gcd(a, b) = gcd(b , a%b) ,这样,我们就可以在几乎是 log 的时间复杂度里求解出来 a 和 b 的最大公约数了,这就是欧几里德算法,用 C++ 语言描述如下:
//欧几里得算法
int gcd(int a, int b) {
return b? gcd(b, a%b) : a;
}
那么什么是扩展欧几里德呢?
a 和 b 的最大公约数是 gcd ,那么,我们一定能够找到这样的 x 和 y ,使得: a*x + b*y = gcd
现在,我们知道了一定存在 x 和 y 使得 : a*x + b*y = gcd , 那么,怎么求出这个特解 x 和 y 呢?只需要在欧几里德算法的基础上加点改动就行了。
我们观察到:欧几里德算法停止的状态是: a= gcd , b = 0 ,那么,这是否能给我们求解 x y 提供一种思路呢?因为,这时候,只要 a = gcd 的系数是 1 ,那么只要 b 的系数是 0 或者其他值(无所谓是多少,反正任何数乘以 0 都等于 0 但是a 的系数一定要是 1),这时,我们就会有: a*1 + b*0 = gcd
当然这是最终状态,但是我们是否可以从最终状态反推到最初的状态呢?
假设当前我们要处理的是求出 a 和 b的最大公约数,并求出 x 和 y 使得 a*x + b*y= gcd ,而我们已经求出了下一个状态:b 和 a%b 的最大公约数,并且求出了一组x1 和y1 使得: b*x1 + (a%b)*y1 = gcd , 那么这两个相邻的状态之间是否存在一种关系呢?
我们知道: a%b = a - (a/b)*b(这里的 “/” 指的是整除,例如 5/2=2 , 1/3=0),那么,我们可以进一步得到:
gcd = b*x1 + (a-(a/b)*b)*y1
= b*x1 + a*y1 – (a/b)*b*y1
= a*y1 + b*(x1 – a/b*y1)
对比之前我们的状态:求一组 x 和 y 使得:a*x + b*y = gcd ,是否发现了什么?
这里:
x = y1
y = x1 – a/b*y1
以上就是扩展欧几里德算法的全部过程,依然用递归写:
//扩展欧几里得
int Extend_gcd(int a, int b, int &x, int &y) {
if(b == 0) {
x = 1;
y = 0;
return a;
}
int ans = Extend_gcd(b, a%b, x, y);
int temp = x;
x = y;
y = temp - a/b*y;
return ans;
}
欧几里得可以求两个数的最大公约数,扩展欧几里德算法的应用主要有以下三方面:
(1)求解不定方程;
(2)求解模线性方程(线性同余方程);
(3)求解模的逆元;
下面的证明过程转载自jumping_frog,解释的特别详细,稍微整理综合了一下,我感觉扩展欧几里得解决的这三个问题其实是同一问题,只是一开始的表达形式不一样,经过转化以后其实本质都是一样的,就是求解不定方程
(1)求解不定方程
对于不定整数方程ax+by=c,若 c mod Gcd(a, b)=0,则该方程存在整数解,否则不存在整数解。
设 d= gcd(a,b),假如整数 x 和 y,满足 d= ax+ by(用扩展欧几里德得出),
如果 d | c(意思是c可以整除d),则方程 a* x0+ n* y0= d(最基本的扩展欧几里得公式), 方程两边乘以 c / d,(因为 d | c,所以能够整除,正是因为可以整出才能这样乘,不然如果c不能被d整除,那么a* x0+ n* y0= d的右边的d乘以c / d就没办法约成c了,在实际的数学中可以,而在计算机中/除法是整除,结果是整数,也就是说不一定是准确的除,就是说如果结果为1.5,但是实际的值为1,因为只保留整数部分(前提是进行除法运算的两个数是全是整数),所以如果像这样用的话,必须保证d | c),得到 a* x0* c/ d+ n* y0* c/ d= c, 这样就转换成了一开始的这个式子的形式ax+by=c,对应位置相对应就行,x= x0* c/ d,y= y0* c/ d 。所以 x= x0* c/d,y= y0* c/ d 为 ax+by=c 的一个解,其余的解为:
//不定整数方程pa+qb=c,求出p和q
int linear_equation(int a, int b, int c, int &x, int &y) {
int gcd = Extend_gcd(a, b, x, y);
if(c % gcd) //只有满足该条件该不定方程才有解
return 0;
x *= c / gcd; //求的只是其中的一组解
y *= c / gcd;
//int x1 = x + b/gcd * t; t为任意整数,x1,y1就是所有的整数解
//int y1 = y - a/gdc * t;
return 1;
}
(2)用扩展欧几里德算法求解模线性方程的方法:
同余方程 ax≡b (mod n)对于未知数 x 有解,当且仅当 gcd(a,n) | b。且方程有解时,方程有 gcd(a,n) 个解。
求解方程 ax≡b (mod n) 相当于求解方程 ax+ ny= b, (x, y为整数,x, y可能为负数)
设 d= gcd(a,n),假如整数 x 和 y,满足 d= ax+ ny(用扩展欧几里德得出)。如果 d| b,则方程
a* x0+ n* y0= d, 方程两边乘以 b/ d,(因为 d|b,所以能够整除),得到 a* x0* b/ d+ n* y0* b/ d= b。
所以 x= x0* b/ d,y= y0* b/ d 为 ax+ ny= b 的一个解,所以 x= x0* b/ d 为 ax+ ny= b的解。所以,
ax≡b (mod n)的一个解为 x0= x* (b/ d ) mod n(注意这里模了一个n,不能忘了),且方程的 d=gcd(a, n) 个解分别为 xi= (x0+ i* (n/ d ))mod n {i= 0... d-1}。
设ans=x*(b/d),s=n/d;
方程ax≡b (mod n)的最小整数解为:(ans%s+s)%s;
相关证明:
证明方程有一解是: x0 = x'(b/d) mod n;
由 a*x0 = a*x'(b/d) (mod n)
a*x0 = d (b/d) (mod n) (由于 ax' = d (mod n))
= b (mod n)
证明方程有d个解: xi = x0 + i*(n/d) (mod n);
由 a*xi (mod n) = a * (x0 + i*(n/d)) (mod n)
= (a*x0+a*i*(n/d)) (mod n)
= a * x0 (mod n) (由于 d | a)
= b
首先看一个简单的例子:
5x=4(mod3)
解得x = 2,5,8,11,14.......
由此可以发现一个规律,就是解的间隔是3.
那么这个解的间隔是怎么决定的呢?
如果可以设法找到第一个解,并且求出解之间的间隔,那么就可以求出模的线性方程的解集了.
我们设解之间的间隔为dx.
那么有
a*x = b(mod n);
a*(x+dx) = b(mod n);
两式相减,得到:
a*dx(mod n)= 0;
也就是说a*dx就是a的倍数,同时也是n的倍数,即a*dx是a 和 n的公倍数.为了求出dx,我们应该求出a 和 n的最小公倍数,此时对应的dx是最小的.
设a 和 n的最大公约数为d,那么a 和 n 的最小公倍数为(a*n)/d.
即a*dx = a*n/d;
所以dx = n/d.
因此解之间的间隔就求出来了.
代码如下:
//求同余方程 ax≡b (mod n), x的解的个数
int modular_linear_equation(int a, int b, int n) {
int x, y, x0, i;
int gcd = Extend_gcd(a, n, x, y);
if(b % gcd) //同余方程没有解
return 0;
x0 = x * b / gcd % n; //特解
for(i=1; i<gcd; ++i)
printf("%d ",(x0 + i * n / gcd) % n);
return 1;
}
(3)用欧几里德算法求模的逆元:
同余方程ax≡b (mod n),如果 gcd(a,n)== 1,则方程只有唯一解。
在这种情况下,如果 b== 1,同余方程就是 ax=1 (mod n ),gcd(a,n)= 1。
这时称求出的 x 为 a 的对模 n 乘法的逆元。
对于同余方程 ax= 1(mod n ), gcd(a,n)= 1 的求解就是求解方程
ax+ ny= 1,x, y 为整数。这个可用扩展欧几里德算法求出,原同余方程的唯一解就是用扩展欧几里德算法得出的 x 。
求乘法逆远的代码: