浅谈 exgcd & 二元一次不定方程
linch
·
·
题解
Update on 2025/10/12：修改一处警示框中的笔误（不影响算法学习）。
本文使用严谨的数学计算，详细推导每一个结论。在阅读之前，请您先学习前置知识，否则可能在阅读时遇到困难。
:::info[前置知识：欧几里得定理[^1]]
证明见 [最大公约数 - OI Wiki](https://oi-wiki.org//math/number-theory/gcd/)。
:::
:::info[前置知识：裴蜀定理[^2]]
当 $a,b\in \mathbb Z$ 时，
- 对于 $\forall x,y\in \mathbb Z$，均有 $\gcd(a,b)\mid (ax+by)$。
- $\exists x,y \in \mathbb Z$，使得 $ax+by=\gcd(a,b)$。
证明见 [裴蜀定理 & 一次不定方程 - OI Wiki](https://oi-wiki.org/math/number-theory/bezouts/)，模板题 [P4549](https://www.luogu.com.cn/problem/P4549)。
:::
## 算法介绍
### 扩展欧几里德算法 exgcd
可以通过该算法求 $ax+by=\gcd(a,b)$ 的一组可行解。
对于方程 $ax+by=\gcd(a,b)$，由欧几里得定理得 $ax+by=\gcd(a,b)=\gcd(b,a\bmod b)$。
设 $x',y'$ 为 $bx+(a\bmod b)y=\gcd(b,a\bmod b)$ 一组解，则有：
$$
\begin{align}
ax+by &=\gcd(b,a\bmod b)\nonumber\\
&=bx'+(a\bmod b)y'\nonumber\\
&=bx'+(a-\left \lfloor \frac{a}{b} \right \rfloor \times b)y'\nonumber\\
&=bx'+ay'-\left \lfloor \frac{a}{b} \right \rfloor \times by'\nonumber\\
&=ay'+b(x'-\left \lfloor \frac{a}{b} \right \rfloor \times y')\nonumber
\end{align}
$$
于是有 $ax+by=ay'+b(x'-\left \lfloor \frac{a}{b} \right \rfloor \times y')$，那么可以得出
$$
\begin{cases}
x=y'\\y=x'-\left \lfloor \frac{a}{b} \right \rfloor \times y'
\end{cases}
$$
递归边界为 $b\mid a$，此时原方程化为 $ax+by=b$，则解为 $x=0,y=1$。每次将 $ax+by=\gcd(a,b)$ 递归至 $bx+(a\bmod b)y=\gcd(b,a\bmod b)$ 求解后回到上一层计算结果即可。
### 解方程
已经使用上述方法通过 exgcd 求出了 $ax+by=\gcd(a,b)$ 的解并获得了 $\gcd(a,b)$ 的值。
由裴蜀定理得 $\gcd(a,b)\mid (ax+by)$，故 $ax+by=c$ 有整数解与 $\gcd(a,b) \mid c$ 互为充要条件。当 $\gcd(a,b) \nmid c$ 时即为无整数解。
接下来考虑当 $\gcd(a,b) \mid c$ 时，对于 exgcd 计算出解的方程 $ax+by=\gcd(a,b)$，设一组解为 $x_0,y_0$，那么在其两边同乘 $\frac{c}{\gcd(a,b)}$ 有 $a\times \frac{c}{\gcd(a,b)}x_0+b\times \frac{c}{\gcd(a,b)}y_0=c$，容易得出方程 $ax+by=c$ 的一组特解为 $x=\frac{c}{\gcd(a,b)}x_0,y=\frac{c}{\gcd(a,b)}y_0$。
对于一个整数 $x_1$，满足要求的 $y$ 解为 $\frac{c-x_1\times a}{b}$，那么 $x$ 越大，$y$ 越小，反之同理。容易发现对 $ax$ 和 $by$ 增加、减少相同值后仍是一个解。假设分别需要增加（减少） $s$，则有 $(ax+s)+(by-s)=c$ 成立。此时将其转化为 $ax+by=c$ 形式的方程：$a(x+\frac{s}{a})+b(y+\frac{s}{b})=c$，则此时解为：
$$
\begin{cases}
x'=x+\frac{s}{a}\\
y'=y-\frac{s}{b}\\
\end{cases}
$$
要求解是整数，那么需要 $a\mid s$ 且 $b\mid s$。因此 $s$ 既是 $a$ 的倍数，又是 $b$ 的倍数，那么 $s$ 一定为 $\operatorname{lcm}(a,b)$ 的倍数。下面设 $d=\operatorname{lcm}(a,b)$，那么 $s$ 一定可以写成 $s=pd$ 的形式，此时有
$$
\begin{align}
x'&=x+\frac{pd}{a}\nonumber\\&=x+p\times\frac{d}{a}\nonumber\\&=x+p\times\frac{d}{a}\nonumber\\&=x+p\times \frac{ab}{\gcd(a,b)\times a}\nonumber\\&=x+p\times \frac{b}{\gcd(a,b)}\nonumber
\end{align}
$$
可以得出 $x'$ 可以由 $x$ 增加或减少若干个 $\frac{b}{\gcd(a,b)}$ 得来。同理，$y'$ 可以由 $y$ 增加或减少若干个 $\frac{a}{\gcd(a,b)}$ 得来。
不妨设 $d_x=\frac{b}{\gcd(a,b)}$，$d_y=\frac{a}{\gcd(a,b)}$。
可得出通解形式为
$$
\begin{cases}
x'=x+p\times d_x\\
y'=y-p\times d_y\\
\end{cases}
$$
根据该通解求出问题答案即可，求解过程会在【代码实现】板块详细说明。
## 正确性证明
上面已经详细推出了每一个结论，接下来证明一下 exgcd 时间复杂度是 $O(\log n)$ 的[^3]。
下面设数对 $(a,b)$ 为此时递归的 $a,b$ 值，$(a',b')$ 为操作一次后的值。分类讨论：
- $a<b$：显然此时经过一次操作后变为 $(b,a)$，即转为 $a>b$ 的情况。
- $b\mid a$：直接返回。
- $a>b 且 b\nmid a$：可以通过不超过 $2$ 次操作使其规模（即 $\max(a,b)$，由于 $a\bmod b<b$，因此一定为 $a$）减少一半。
:::info[证明]{open}
首先 $a>2b$ 时，对于操作后数对 $(b,a\bmod b)$，有 $2a'=2b<a$，$a'$ 的值一定较 $a$ 缩小至少一半，原命题得证。
再考虑 $b<a<2b$ 时，此时 对于操作后数对 $(b,a\bmod b)$ 即 $(b,a-\left \lfloor \frac{a}{b} \right \rfloor\times b)$ 变为 $(b,a-b)$，再操作一次 $a''=b'=a-b$ 由于 $b>\frac{a}{2}$，那么 $a-b<\frac{a}{2}$，即 $a''$ 至少较 $a$ 减小一半，命题得证。
:::
按照该规模减小速度，于是可以得出时间复杂度为 $O(\log n)$。
## 代码实现
### 核心实现
exgcd 部分，按照上述方法，对于 $(a,b)$，递归至 $(b,a\bmod b)$ 查询后，解出
$$
\begin{cases}
x=y'\\y=x'-\left \lfloor \frac{a}{b} \right \rfloor \times y'
\end{cases}
$$
即可。递归到 $b\mid a$ 时解得 $x=0,y=1$ 返回。实现上，$x,y$ 使用 `&` 标识符方便地获取到返回的解计算。
:::success[exgcd 部分代码]
```cpp line-numbers
int exgcd(int a,int b,long long &x,long long &y){
if(a%b==0){
x=0,y=1;
return b;
}
int k=exgcd(b,a%b,x,y);
int t=x;
x=y;
y=t-a/b*y;
return k;
}
```
:::
接下来计算一下要求输出的答案。
首先，当 $\gcd(a,b) \nmid c$ 时即为无整数解，输出 `-1` 即可。
然后当 $\gcd(a,b) \mid c$ 时，对于 exgcd 计算出解的方程 $ax+by=\gcd(a,b)$，设一组解为 $x_0,y_0$，根据【算法说明】部分的推导，得出方程 $ax+by=c$ 的一组特解为 $x=\frac{c}{\gcd(a,b)}x_0,y=\frac{c}{\gcd(a,b)}y_0$。
计算出一组特解后，根据 $d_x=\frac{b}{\gcd(a,b)}$，$d_y=\frac{a}{\gcd(a,b)}$，得出通解形式
$$
\begin{cases}
x'=x+p\times d_x\\
y'=y-p\times d_y\\
\end{cases}
$$
考虑正整数解，
我们限制 $x'>0$，则有
$$
\begin{align}
x+p\times d_x&>0\nonumber\\
p\times d_x&>-x\nonumber\\
p&>-\frac{x}{d_x}\nonumber
\end{align}
$$
由于 $p\in \mathbb Z$，故有 $p\ge \left \lceil \frac{1-x}{d_x} \right \rceil$。同理，限制 $y'>0$ 可推出 $p<\frac{y}{d_y}$，于是有 $p\le \left \lfloor \frac{y-1}{d_y} \right \rfloor$，那么正整数解的限制为 $\left \lceil \frac{1-x}{d_x} \right \rceil\le p \le \left \lfloor \frac{y-1}{d_y} \right \rfloor$。
:::warning[注意]{open}
注意不是 $\left \lceil -\frac{x}{d_x} \right \rceil\le p \le \left \lfloor \frac{y}{d_y} \right \rfloor$，因为这样在 $d_x \mid x$ 或 $d_y \mid y$ 时会出现 $p$ 取等，那么解就是 $0$，不是正整数。
:::
那么当 $\left \lceil \frac{1-x}{d_x} \right \rceil > \left \lfloor \frac{y-1}{d_y} \right \rfloor$ 时，上述不等式 $p$ 无解，也就是没有正整数解。显然 $p$ 越小，$x$ 越小。因此此时 $x$ 的最小正整数解即为 $p$ 取最小值即 $p=\left \lfloor -\frac{x}{d_x} \right \rfloor$ 时，得出 $x_{\min}=x+\left \lfloor -\frac{x}{d_x} \right \rfloor \times d_x$。同理，$p$ 越大，$y$ 越小。$y$ 的最小正整数解即为 $p$ 取最大值即 $p= \left\lfloor \frac{y-1}{d_y} \right \rfloor$ 时 $y_{\min}=y-\left\lfloor \frac{y-1}{d_y} \right \rfloor \times d_y$。
否则，便存在 $p$ 使得 $x,y$ 均为正整数。由于 $p$ 越大，$x$ 越大，$y$ 越小，那么 $x$ 的最大正整数解和 $y$ 的最小正整数解便取在 $p$ 的最大值，则有 $x_{\max}=x+\left\lfloor \frac{y-1}{d_y} \right \rfloor \times d_x,y_{\min}=y-\left\lfloor \frac{y-1}{d_y} \right \rfloor \times d_y$。反之，$p$ 越小，$x$ 越小，$y$ 越大，$x$ 的最小正整数解和 $y$ 的最大正整数解便取在 $p$ 的最小值，则有 $x_{\min}=x+\left \lceil \frac{1-x}{d_x} \right \rceil \times d_x,y_{\max}=y-\left \lceil \frac{1-x}{d_x} \right \rceil \times d_y$。每个 $p$ 都对应一组解，总正整数解数即为最大的 $p$ 和最小的 $p$ 之差。
根据上述公式求出结果即可。
### 常见错误
1. 取整。本题中存在负数除法取整，由于 C++ 默认向 $0$ 取整，而不是向下取整，因此需要使用 `ceil` 和 `floor` 函数，而不能使用正数的写法。
:::error[错误的写法]
- 向上取整：`(x-1)/y+1`；
- 向下取整：`x/y`。
:::
:::success[正确的写法]
- 向上取整：`ceil(1.0*x/y)`；
- 向下取整：`floor(1.0*x/y)`。
:::
2. 开 `long long`。
3. 各种不等号是否取等，特别是上文警告框中提到的 $p$ 范围。
4. 整数除法取整需要强制转换类型。
:::success[[AC](https://www.luogu.com.cn/record/231841609) 代码]
```cpp line-numbers
#include<bits/stdc++.h>
using namespace std;
long long a,b,c;
int exgcd(int a,int b,long long &x,long long &y){
if(a%b==0){
x=0,y=1;
return b;
}
int k=exgcd(b,a%b,x,y);
int t=x;
x=y;
y=t-a/b*y;
return k;
}
int main(){
ios::sync_with_stdio(false);
cin.tie(0);
int t;
cin>>t;
while(t--){
cin>>a>>b>>c;
long long x=0,y=0,z;
z=exgcd(a,b,x,y);
if(c%z!=0){
cout<<"-1\n";
continue;
}
long long dx=b/z,dy=a/z;
x*=c/z,y*=c/z;
long long m=ceil(1.0*(1-x)/dx),n=floor(1.0*(y-1)/dy);
if(m>n) cout<<x+m*dx<<" "<<y-n*dy<<"\n";
else{
long long xmn=(x%dx+dx)%dx;
xmn=xmn==0?dx:xmn;
long long ymn=(y%dy+dy)%dy;
ymn=ymn==0?dy:ymn;
long long xmx=(c-ymn*b)/a,ymx=(c-xmn*a)/b,cnt=n-m+1;
cout<<cnt<<" "<<xmn<<" "<<ymn<<" "<<xmx<<" "<<ymx<<"\n";
}
}
return 0;
}
```
:::
## 参考文献
[^1]:[最大公约数 - OI Wiki](https://oi-wiki.org//math/number-theory/gcd/)
[^2]:[裴蜀定理 & 一次不定方程 - OI Wiki](https://oi-wiki.org/math/number-theory/bezouts/)
[^3]:[辗转相除的时间复杂度 - Vegdie - 博客园](https://www.cnblogs.com/CYLSY/p/16816154.html)