从gcd到exgcd
注:如果方程设成ax-by=0之类的,判断-b即可。详见本文。
Part 0 前言
gcd以及exgcd算法是oi数论的基础,不管是求逆、解不定方程还是扩展中国剩余定理等都需要用到这里的知识。
Part 1 引入——求解gcd
欧几里得算法
我们~小学二年级~就学过最大公约数,并了解过一种计算gcd的方法——辗转相除法。其定义如下:
| int gcd(int a,int b){
if(b==0) return a;
return gcd(b,a%b);
}
|
该算法用到了一个重要的性质:
gcd(a,b)=gcd(b,amodb)
笔者太菜,有关证明详见 OI-wiki-最大公约数
而exgcd,便是以这段代码为基础,转化而来。它可以求解形如ax+by=c的方程(整数域中)。我们先学习一下怎样判断其无解。
裴蜀定理
设 a,b 是不全为零的整数,则存在整数x,y,使得 ax+by=gcd(a,b)
这句话告诉我们什么道理呢?
- ax+by=c 这个方程有整数解,当且仅当 gcd(a,b)∣c
- 令f=ax+by,在整数域中,使得 f>0 且 f 尽可能的小,则fmin=gcd(a,b)
对于多元有推论:
- n元一次线性方程 i=1∑nAi×Xi=C有整数解,当且仅当 gcd(A1,A2,...,An)∣C
- 在整数域中,使得 C>0 且 C 尽可能的小,则C=gcd(A1,A2,...,An)
当你学会这个推论,那么便可以解决这道题了:
P4549 【模板】裴蜀定理
Part 2 推导
exgcd算法
别被标题吓到,其实这个挺简单的。 学会推式子,不会忘了咋打~
注意:标准的 exgcd 只能用于求解方程ax+by=gcd(a,b),若c是gcd(a,b)的倍数,我们在后文讨论。
如题,我们有ax+by=gcd(a,b)
考虑进行一次迭代,带入a=b,b=a%b,我们又有:bx′+(a%b)y′=gcd(b,a%b)
注意这里的x′,y′系数已经改变,新方程的解与原方程x,y不一样,故另作标号。
类似gcd,exgcd经过不断迭代递归,最后也会到达状态为exgcd(a,0)的形式,即y系数为0,x系数为a,方程为ax+0y=gcd(a,0)=a,故此时x=1,y=0
那么如何通过这个解求得上一个方程的x,y?
我们知道,gcd(a,b)=gcd(b,a%b) ,所以 ax+by=bx′+(a%b)y′
又a%b=a−b∗⌊a/b⌋,所以有:
ax+by
=bx′+(a−b∗⌊a/b⌋)y′
=bx′+ay′−b∗⌊a/b⌋y′
=ay′+b(x′−⌊a/b⌋y′)
待定系数,我们有x=y′,y=x′−⌊a/b⌋y′
由此,我们便可以一层一层向上推出原方程的一组解。
代码如下:
| LL exgcd(LL a, LL b, LL& x, LL& y) {
if (b == 0) {
x = 1, y = 0;
return a;
}
LL px, py; // 上一次的答案(x'和y')
LL ans = exgcd(b, a % b, px, py);
x = py, y = px - a / b * py;
return ans;
}
|
这里可以直接在exgcd中求得gcd(a,b)的值,但是若传入的参数a,b为负数,则gcd有可能是负数。
推导通解
那么已知一组特解(x0,y0),如何才能求出方程的通解呢?或者限定x或y的值为最小正整数,该如何求解?下面是通解的推导:
ax0+by0=gcd(a,b)
⇒ax0+k∗lcm(a,b)+by0−k∗lcm(a,b)=gcd(a,b)
我们知道lcm(a,b)=ab/gcd(a,b),故上式可化为:
ax0+a(k∗b/gcd)+by0−b(k∗a/gcd)=gcd(a,b)
⇒a(x0+k∗b/gcd)+b(y0−k∗a/gcd)=gcd(a,b)
待定系数得:{x=x0+k∗b/gcdy=y0−k∗a/gcd
由此,我们便有了求解x或y的最小正整数解的方法:
注意到式子形如x=a+kb,也即x≡x0(modb/gcd)(这对于y也同理)
设p=∣b/gcd∣,则x的最小整数解即为x=(x0∗gcd%p+p)%p
更一般的情况
上一节推导了ax+by=gcd(a,b)的通解,而一般往往求解的方程是ax+by=c=n∗gcd(a,b)
聪明的你会说,给算出来的x和y乘个n不就得了,即:
{x′=nx=n∗(x0+k∗b/gcd)y′=ny=n∗(y0−k∗a/gcd)
但是,如果你直接给上一节算出的最小整数解乘n,其实不一定是该方程的最小整数解:
{x′=nx0+kn∗b/gcd→nx0+q∗b/gcdy′=ny0−kn∗a/gcd→ny0−q∗a/gcd
在构造x′和y′时,我们只要求kn是整数,故k实际上可以取任意实数,可以一并视为一个倍数q。
从集合的角度看,其实下面这组解集是包含上面那组解集的。
我们不难得到最小整数解:
x′=(x0∗(c/gcd)%p+p)%p,其中p=∣b/gcd∣
那如果a,b有负数呢?我们不妨假设a>0,b>0:
−ax+by=n∗gcd(a,b)
我们有a(−x/n)+b(y/n)=gcd(a,b),说明exgcd求得的特解x0=−x/n,y0=y/n,移项解出x或y,再利用取模的公式即可(下文有图像)。
代码实现
1
2
3
4
5
6
7
8
9
10
11
12
13
14 | #include <cmath>
// 求解形如ax+by=0
// 若形如ax-by=0,系数变为a,-b
LL a, b, c, X0, Y0;
LL gcd = exgcd(abs(a), abs(b), X0, Y0); // gcd一定要是正值
if (c % gcd != 0) continue; // 先判方程无解
LL p = abs(b / gcd);
if (a < 0) X0 = -X0; // 系数为负 解也为负
// if (b < 0) Y0 = -Y0;
// 取决于要取哪个值,算一个带入求另一个即可
LL x = (X0 * c / gcd % p + p) % p;
LL y = (c - a * x) / b;
|
Part 3 几何意义与一些问题
裴蜀定理到底说了个啥?为啥取模求最小正整数的方法对于所有情况都适用?
我们不妨画出直线ax+by=c
Q:不符合裴蜀定理的情况——与整数格点无交点

Q:为什么模数偏偏是b/gcd,a/gcd?
eq2为原方程,eq1为exgcd所求的方程,可以发现,eq1取(1,0),乘倍数3后并不是x为最小整数解的点。
并且可以观察到,转化前后直线斜率相同,斜率三角形相同,故一旦确定一个格点(x0,y0),则每一个解的x,y坐标在直线上挪动的距离也相同,所以模数也相同。

Q:(x,y)是否都能取正整数?
我们可以先算出x的最小正整数解,验证对应的y是否大于0;如小于0,则再算y的最小正整数,求出对应的x是否大于0(一定要二次确认,会存在y>0而x不是最小正整数解的情况)

(下面这两幅图配的不太对) 凑合着看看,建议手画理解一下。

