1. 前言
扩展欧几里得(exgcd),是在欧几里得算法基础上求解任意形如 a x + b y = c ax+by=c ax+by=c 的二元一次方程的一组特解的一种算法。
在往下看之前,您只需要知道如何使用欧几里得算法求 gcd ( a , b ) \gcd(a,b) gcd(a,b)。
不知道也没关系,式子在这里:
gcd ( a , b ) = gcd ( b , a m o d b ) \gcd(a,b)=\gcd(b, a \bmod b) gcd(a,b)=gcd(b,amodb)
证明网络上面有很多,可以自行查找。
那么 exgcd \text{exgcd} exgcd 就是其扩展算法。
2. 详解
例题 1:给出 a , b ∈ N + , d = gcd ( a , b ) a,b \in N_+,d=\gcd(a,b) a,b∈N+,d=gcd(a,b),求方程 a x + b y = d ax+by=d ax+by=d 的一组整数解。
这就是 exgcd \text{exgcd} exgcd 要解决的问题。
接下来给出证明 a x + b y = d ax+by=d ax+by=d 一定有整数解的证明过程,该过程同时也给出了 exgcd \text{exgcd} exgcd 的求解过程。
考虑普通 gcd \gcd gcd 的求解过程: gcd ( a , b ) = gcd ( b , a m o d b ) \gcd(a,b)=\gcd(b,a \bmod b) gcd(a,b)=gcd(b,amodb)。
假设我们现在已经知道了一组解 x ′ , y ′ x',y' x′,y′ 是方程 b x + ( a m o d b ) y = d bx + (a \bmod b)y=d bx+(amodb)y=d 的一组特解,那么如何求出 a x + b y = d ax+by=d ax+by=d 的解呢?
首先我们知道这个式子成立: a m o d b = a − ⌊ a b ⌋ × b a \bmod b=a - \left\lfloor\dfrac{a}{b}\right\rfloor \times b amodb=a−⌊ba⌋×b
那么把上面这个式子带入到 b x + ( a m o d b ) y = d bx+(a \bmod b)y=d bx+(amodb)y=d 中,就有:
b x ′ + ( a − ⌊ a b ⌋ × b ) y ′ = d bx'+(a-\left\lfloor\dfrac{a}{b}\right\rfloor \times b)y'=d bx′+(a−⌊ba⌋×b)y′=d
拆掉括号:
b x ′ + a y ′ − ⌊ a b ⌋ × b × y ′ = d bx' + ay' -\left\lfloor\dfrac{a}{b}\right\rfloor \times b \times y' = d bx′+ay′−⌊ba⌋×b×y′=d
左边转化一下:
a y ′ + b ( x ′ − ⌊ a b ⌋ y ′ ) = d ay'+b(x'-\left\lfloor\dfrac{a}{b}\right\rfloor y')=d ay′+b(x′−⌊ba⌋y′)=d
上述方程跟下述方程是等价的:
a x + b y = d ax+by=d ax+by=d
两者对比,得到:
{ x = y ′ y = x ′ − ⌊ a b ⌋ y ′ \begin{cases}x=y'\\y=x'-\left\lfloor\dfrac{a}{b}\right\rfloor y'\end{cases} {x=y′y=x′−⌊ba⌋y′
于是我们成功的通过 b x ′ + ( a m o d b ) y ′ = d bx'+(a \bmod b)y'=d bx′+(amodb)y′=d 推出了 a x + b y = d ax+by=d ax+by=d 的解。
那么为什么一定会有解呢?
考虑欧几里得算法求 gcd \gcd gcd 的步骤: gcd ( a , b ) = gcd ( b , a m o d b ) \gcd(a,b)=\gcd(b,a \bmod b) gcd(a,b)=gcd(b,amodb)。
仿照上述过程求 exgcd \text{exgcd} exgcd,那么初始状态是 a ′ x + b ′ y = d a'x+b'y=d a′x+b′y=d,其中 a = d , b = 0 a=d,b=0 a=d,b=0。
因此 a = 1 , b = 0 a=1,b=0 a=1,b=0。
然后根据上述分析不断还原,最后一定能够还原出 a x + b y = d ax+by=d ax+by=d 的解。
证毕。
上述过程给出了 a x + b y = d ax+by=d ax+by=d 一定有解的证明,同时阐述了为什么有解。
同时根据上述过程应该明白为什么右边是 gcd ( a , b ) = d \gcd(a,b)=d gcd(a,b)=d 了吧~这跟 gcd \gcd gcd 的求法有很大关系。
代码:
void exgcd(int a, int b, LL &x, LL &y)//注意 x 跟 y 是引用
{
if (b == 0) {x = 1; y = 0; return ;}
exgcd(b, a % b, x, y);
LL p = x; x = y; y = p - ((LL)a / b) * y;//特别注意存一下 x 的初始值
}
例题 2:给出 a , b ∈ N + , d = gcd ( a , b ) a,b \in N_+,d = \gcd(a,b) a,b∈N+,d=gcd(a,b),求解 a x + b y = d ax+by=d ax+by=d 中 x x x 可以取到的最小正整数值。
首先显然还是需要利用 exgcd \text{exgcd} exgcd 求出特解 x 0 , y 0 x_0,y_0 x0,y0,但是怎样求出 x x x 的最小正整数值呢?
首先根据数学知识,我们知道通解如下:
{ x = x 0 − t b y = y 0 + t a , t ∈ Z \begin{cases}x=x_0-tb\\y=y_0+ta\end{cases},t \in Z {x=x0−tby=y0+ta,t∈Z
于是如果解出来的 x 0 < 0 x_0 <0 x0<0,那么不断加上 t b tb tb,否则不断减去 t b tb tb 就可以了。
而上述过程可以巧妙的使用下面这一句代码来解决:
x = (x % b + b) % b;
于是只需要在做完 exgcd \text{exgcd} exgcd 之后做一遍上述工作即可,就一行代码。
例题 3:求例题 2 方程的正整数解个数。
根据数学知识,当 x x x 取到最小值时 y y y 取到最大值,而这中间的变化量可以计算得到。
同理, y y y 取到最小值时 x x x 取到最大值。
设 x x x 最大为 x max ∈ N + x_{\max} \in N_+ xmax∈N+,最小为 x min ∈ N + x_{\min} \in N_+ xmin∈N+,那么正整数解的个数就是:
x max − x min b + 1 \dfrac{x_{\max}-x_{\min}}{b}+1 bxmax−xmin+1
很好理解吧~
例题 4:给出
a
,
b
,
c
∈
N
+
a,b,c \in N_+
a,b,c∈N+,求
a
x
+
b
y
=
c
ax+by=c
ax+by=c 的一组特解
x
0
,
y
0
∈
Z
x_0,y_0 \in Z
x0,y0∈Z,若无解输出 -1
。
首先需要判无解,而根据裴蜀定理可以知道: gcd ( a , b ) ∣ c \gcd(a,b) \mid c gcd(a,b)∣c 时有解,否则无解。
关于裴蜀定理参考这篇博文:妈妈我竟然证出了裴蜀定理!
判无解之后考虑如何解决这道题。
首先求出这个方程的特解: a x + b y = gcd ( a , b ) ax+by=\gcd(a,b) ax+by=gcd(a,b) 的解 x ′ , y ′ x',y' x′,y′。
因为 c = c gcd ( a , b ) × gcd ( a , b ) c=\dfrac{c}{\gcd(a,b)} \times \gcd(a,b) c=gcd(a,b)c×gcd(a,b),因此原方程的解就是 x = x ′ c gcd ( a , b ) , y = y ′ c gcd ( a , b ) x=\dfrac{x'c}{\gcd(a,b)},y=\dfrac{y'c}{\gcd(a,b)} x=gcd(a,b)x′c,y=gcd(a,b)y′c。
做完了。
3. 总结
其实 exgcd \text{exgcd} exgcd 还是很好写的,公式也好推,会推就可以了。