跳到主要内容
原文884 词数4 分钟

题解:P5656 【模板】二元一次不定方程 (exgcd)

概述

求解二元一次不定方程 ax+by=c:无整数解输出 -1;有正整数解时输出正整数解的个数以及 x、y 的最小值与最大值;有整数解但无正整数解时输出 x、y 的最小正值。先用扩展欧几里得算法求出 gcd(a,b) 和一组整数解,由裴蜀定理判无解,再把解按通解公式平移,用取模求出 x、y 的最小正值,进而得到最大值与解的个数。文中一并推导欧几里得算法、裴蜀定理与扩展欧几里得算法。

参考资料​

题意简述​

给定 1≤a,b,c≤1091\le a,b,c\le 10^9,求解不定方程:

ax+by=cax+by=c
  • 若无整数解,输出 -1;
  • 若有正整数解,输出正整数解中,解的数量、xx 的最小值、yy 的最小值、xx 的最大值、yy 的最大值;
  • 若无正整数解,输出 xx 的最小正数值和 yy 的最小正数值。

基础知识​

一次不定方程​

一次不定方程(Linear Diophantine Equation)是形如以下形式的方程:

a1x1+a2x2+⋯+anxn=ba_1x_1+a_2x_2+\dots+a_nx_n=b

其中 aia_i 和 bb 都是整数,我们会研究 xix_i 的整数解。

特别地,当 n=2n=2 时为二元一次不定方程:

a1x1+a2x2=ba_1x_1+a_2x_2=b

我们可以调整一下变量名:

ax+by=cax+by=c

裴蜀定理​

裴蜀定理(Bézout's Lemma)给出了二元一次不定方程 有整数解 的充要条件。

对于整数 a,ba,b(不全为 00),存在整数 x,yx,y 使得 ax+by=gcd⁡(a,b)ax+by=\gcd(a,b)。

换句话说,ax+by=cax+by=c 有整数解当且仅当 gcd⁡(a,b)∣c\gcd(a,b)\mid c。

证明详见 Bézout's identity - Wikipedia。

欧几里得算法​

欧几里得算法(辗转相除法)可以求两个整数的 最大公约数(Greatest Common Divisor,GCD)。

gcd⁡(a,b)=gcd⁡(b,a mod b)\gcd(a,b)=\gcd(b,a\bmod b)

证明:设 a=kb+ca=kb+c。若 d∣a,bd\mid a,b,则 d∣cd\mid c;若 d∣b,cd\mid b,c,则 d∣ad\mid a。故 a,ba,b 与 b,cb,c 的公约数集合相同,因此 gcd⁡(a,b)=gcd⁡(b,c)=gcd⁡(b,a mod b)\gcd(a,b)=\gcd(b,c)=\gcd(b,a\bmod b)。

显然 a mod b<ba\bmod b<b,因此可以递归求解。特别地,gcd⁡(a,0)=a\gcd(a,0)=a。

gcd⁡(a,b)={ab=0gcd⁡(b,a mod b)otherwise\gcd(a,b)= \begin{cases} a & b=0 \\ \gcd(b,a\bmod b) & \text{otherwise} \end{cases}
46 Bcpp
ll Gcd(ll a,ll b)
{
return !b?a:Gcd(b,a%b);
}

扩展欧几里得算法​

扩展欧几里得算法(Extended Euclidean Algorithm,EXGCD)可以求 ax+by=gcd⁡(a,b)ax+by=\gcd(a,b) 的 一组整数解。

ax+by=gcd⁡(a,b)  ⟺  ay+b(x−⌊ab⌋y)=gcd⁡(a,b)ax+by=\gcd(a,b)\iff ay+b\left(x-\left\lfloor\frac{a}{b}\right\rfloor y\right)=\gcd(a,b)

设一组整数解为 (x0,y0)(x_0,y_0):

ax0+by0=gcd⁡(a,b)ax_0+by_0=\gcd(a,b) bx1+(a mod b)y1=gcd⁡(b,a mod b)bx_1+(a\bmod b)y_1=\gcd(b,a\bmod b)

根据欧几里得算法:

gcd⁡(a,b)=gcd⁡(b,a mod b)\gcd(a,b)=\gcd(b,a\bmod b)

所以:

ax0+by0=bx1+(a mod b)y1ax_0+by_0=bx_1+(a\bmod b)y_1

又因为:

a mod b=a−⌊ab⌋×ba\bmod b=a-\left\lfloor\frac{a}{b}\right\rfloor\times b

所以:

ax0+by0=bx1+(a−⌊ab⌋×b)y1ax_0+by_0=bx_1+\left(a-\left\lfloor\frac{a}{b}\right\rfloor\times b\right)y_1 ax0+by0=ay1+bx1−⌊ab⌋×by1=ay1+b(x1−⌊ab⌋y1)ax_0+by_0=ay_1+bx_1-\left\lfloor\frac{a}{b}\right\rfloor\times by_1=ay_1+b\left(x_1-\left\lfloor\frac{a}{b}\right\rfloor y_1\right)

所以:

x0=y1,y0=x1−⌊ab⌋y1x_0=y_1,y_0=x_1-\left\lfloor\frac{a}{b}\right\rfloor y_1

将 x1,y1x_1,y_1 不断代入递归求解直至 gcd⁡\gcd 为 00 递归 x=1,y=0x=1,y=0 回去求解。

110 Bcpp
tuple<ll,ll,ll> exgcd(ll a,ll b)
{
if(!b)return {a,1,0};
auto [g,x,y]=exgcd(b,a%b);
return {g,y,x-a/b*y};
}

解题思路​

我们先用 扩展欧几里得算法 求出最大公约数 g=gcd⁡(a,b)g=\gcd(a,b) 和一组整数解 (x0,y0)(x_0,y_0)。

ax0+by0=gax_0+by_0=g

根据裴蜀定理,当 cc 不是 gg 的倍数时无解,输出 -1。

将方程两边同时乘以 cg\frac{c}{g}:

cg⋅ax0+cg⋅by0=cg⋅g=c\frac{c}{g}\cdot ax_0+\frac{c}{g}\cdot by_0=\frac{c}{g}\cdot g=c

此时我们令 x0←cg⋅x0,y0←cg⋅y0x_0\gets\frac{c}{g}\cdot x_0,y_0\gets\frac{c}{g}\cdot y_0,就得到了原问题的一组解 (x0,y0)(x_0,y_0)。

每两组解的 xx 和 yy 分别相差了:

dx=bg,dy=agd_x=\frac{b}{g},d_y=\frac{a}{g}

而所有正整数解分别为:

(xmin⁡,ymax⁡),(xmin⁡+dx,ymax⁡−dy),…,(xmax⁡,ymin⁡)(x_{\min},y_{\max}),(x_{\min}+d_x,y_{\max}-d_y),\dots,(x_{\max},y_{\min})

显然 xmin⁡∈[1,dx],ymin⁡∈[1,dy]x_{\min}\in[1,d_x],y_{\min}\in[1,d_y],可以通过取模得到:

xmin⁡=(x0−1) mod dx+1,ymin⁡=(y0−1) mod dy+1x_{\min}=(x_0-1)\bmod d_x+1,y_{\min}=(y_0-1)\bmod d_y+1

当 xx 最小时 yy 最大,yy 最小时 xx 最大,因此:

xmax⁡=c−bymin⁡a,ymax⁡=c−axmin⁡bx_{\max}=\frac{c-by_{\min}}{a},y_{\max}=\frac{c-ax_{\min}}{b}

因此正整数解的数量为:

n=xmax⁡−xmin⁡dx+1=ymax⁡−ymin⁡dy+1n=\frac{x_{\max}-x_{\min}}{d_x}+1=\frac{y_{\max}-y_{\min}}{d_y}+1

参考代码​

638 Bcpp
#include <bits/stdc++.h>
using namespace std;

using ll=long long;
tuple<ll,ll,ll> exgcd(ll a,ll b)
{
if(!b)return {a,1,0};
auto [g,x,y]=exgcd(b,a%b);
return {g,y,x-a/b*y};
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr);
int T;
cin>>T;
while(T--)
{
ll a,b,c;
cin>>a>>b>>c;
auto [g,x,y]=exgcd(a,b);
if(c%g){cout<<-1<<'\n';continue;}
x*=c/g;y*=c/g;
ll dx=b/g,dy=a/g;
ll x_min=(x%dx+dx-1)%dx+1,y_min=(y%dy+dy-1)%dy+1,x_max=(c-b*y_min)/a,y_max=(c-a*x_min)/b;
if(y_max>0)cout<<(x_max-x_min)/dx+1<<' '<<x_min<<' '<<y_min<<' '<<x_max<<' '<<y_max<<'\n';
else cout<<x_min<<' '<<y_min<<'\n';
}
return 0;
}