[模板] 扩展欧几里得算法详解

算法简介

P.S.扩展欧几里得非常重要，本人在此纠结多时，将经验总结于此。

本算法由欧几里得辗转相除得出：

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

如此简洁又高效的代码可以快速地解决最大公约数的问题，那么我们能否对其增加一些应用呢？

算法应用

扩展欧几里得算法主要应用在下面三个方面：

①求解不定方程ax+by=c

②求解模的逆元ax≡1（mod m）

③求解同余方程ax≡b（mod m）

可以看出，第②项是第③项的特例，第③项可以化为第①项：ax+my=b

算法核心

要求ax+by=c的解，先求出ax+by=gcd(a,b)（为什么？看后面去）

下面使用扩展欧几里得解出这样的一组特解

long long Exgcd(long long a,long long b,long long &x,long long &y) {
if(b==0) {
x=1;
y=0;
return a;
}
long long ans=Exgcd(b,a%b,x,y),tmp=x;
x=y;
y=tmp-a/b*y;
return ans;
}

对于算法进行解释：

(i)若b=0，则gcd(a,b)=a，x=1，y=0（为什么y=0？因为不定方程有无数组解我们只需要求出一组特解，后面再进行解的放大缩小即可）

(ii)若b≠0，由欧几里得辗转相除可得到

gcd=gcd(a,b)=gcd(b,a%b) ①

因为我们要求gcd=ax+by，又由①得

gcd=bx’+(a%b)y’ ②

因为m%n可以表示为m-(m/n)*n （此处的’/’是整除符号）

故将②变形为：

gcd=bx’+(a-(a/b)*b)y’=ay’+b(x’-(a/b)*y’) ③

又因为gcd=ax+by，则有：x=y’，y=x’-(a/b)*y’

以上过程较为复杂，建议读者在草稿纸上推算一遍，并将算法记忆下来。

下面解决之前留下来的一个问题（蓝色部分）：

我们求出了ax+by=gcd(a,b)的一组解，如何解出ax+by=c的所有解呢？

仔细观察这两个方程，就会发现他们相差c/gcd倍

故将原方程乘上c/gcd即可得出答案：

x0=x*c/gcd

y0=y*c/gcd

其中若c%gcd!=0，则无解。

我们可以用这一组特解表示出不定方程的通解：

X=x0+(b/gcd)*t

Y=y0-(a/gcd)*t

下面是证明（伪证）：

ax1+by1=c ①

ax2+by2=c ②

①+②得：

ax1+by1=ax2+by2 ③

合并同类项：

a(x1-x2)=b(y2-y1) ④

同除gcd（因为求的就是a、b的最大公约数，所以必能除尽）：

a/gcd*(x1-x2)=b/gcd*(y2-y1) ⑤

移项得：

x1-x2=b/gcd*(y2-y1)*gcd/a ⑥

x1=x2+b/gcd*(y2-y1)*gcd/a ⑦

令t=(y2-y1)*gcd/a

故x1=x2+b/gcd*t ⑧

原题得证

不定方程有无数组解，因此有些题目要求我们求最小正整数解，步骤如下：

①用Exgcd求出ax+by=gcd(a,b)的解

②若c%gcd!=0则无解

③用上述通解表示最小正整数解：

xmin=((x*c/gcd)%(b/gcd)+(b/gcd))%(b/gcd)

读者也许不解为什么有这么多b/gcd？

因为(x*c/gcd)%(b/gcd)可能为负，需要再加一次再模（膜）一下得到正整数解

long long Cal(long long a,long long b,long long c) { //求出ax+by=c的x的最小正整数解
long long x,y;
long long gcd=Exgcd(a,b,x,y);
if(c%gcd!=0) { //无解
cout<<"Impossible"<<endl;
exit(0);
}
long long b0=b/gcd;
return ((x*c/gcd)%b0+b0)%b0;
}

算法扩展

扩展一：求解一元线性同余方程组

（非中国剩余定理）

对于方程组，我们可以进行方程组合并，求出合并后方程的解，不论方程有多少个，都可以两两解决。

考虑两个方程的情况：

有解的条件是：gcd(m1,m2)|(a1-a2)

若同余方程组为：

x≡r1(mod a1) ①

x≡r2(mod a2) ②

x≡r3(mod a3) ③

……

x≡rn(mod an)

将①②根据扩展欧几里得求出在lcm(a1,a2)内的最小非负整数解x=r1+a1y1，接着将两个方程合并为一个方程：

x’≡x(mod lcm(a1,a2))

则新的r1=r1+a1y1,a1=a1/d*a2

合并过程如下：

考虑①、②方程可得

x=r1+a1x1 ④

x=r2+a2x2 ⑤

由④、⑤

a1x1+a2x2=r2-r1 ⑥（符号归于系数中）

用Exgcd解出特解x0

将x0代入④得

x=a1x0+r1 ⑤

用通解形式表示x可得

x=a1(x0+i*a2/gcd(a1,a2))+r1 ⑥

=a1x0+r1+i*a2/gcd(a1,a2) ⑦

故合并为方程：x’≡(a1x0+r1)(mod lcm(a1,a2))

将该方程与后面方程一一联立即可求解

long long Solve() {
long long a1=a[1],r1=r[1];
for(int i=2; i<=n; i++) {
long long a2=a[i],r2=r[i],a=a1,b=a2,c=r2-r1,x0,y0;
long long gcd=Exgcd(a,b,x0,y0);
if(c%gcd)return -1; //无解
long long b1=b/gcd;

4000
x0=(x0*(c/gcd)%b1+b1)%b1;
r1+=a1*x0; //带入原方程组的第一个
a1*=a2/gcd;
}
return r1;
}

扩展二：求解多元同余方程

即求解：

a1x1+a2x2+a3x3+...+anxn+b≡0(modm)

若有解，则满足：gcd(a1,a2,a3…an,m)|b

解的个数为m^n-1*gcd(a1,a2,a3…an,m)

求解此方程一般不需要掌握，有兴趣的同学可以自行研究。

经典例题

[ZJOI2002] 青蛙的约会

[Vijos1009] 清帝之惑之康熙

[NOIP2012] 同余方程

[poj2115] C Looooops

[poj2891] Strange Way to Express Integers

[hdu1573] X问题

[hdu3579] Hello Kiki

[poj2142] 天平 The Balance

[poj1091] 跳蚤

后记

感谢大家的支持，若有误请指出本人将及时更改。