题解 P5656 【【模板】二元一次不定方程(exgcd)】
linponess
2019-12-10 19:33:57
### **2020/1/25 UPD:更正了证明过程中符号的问题。（By [weapon_zipper](https://www.luogu.com.cn/user/88207)）**
### **2019/12/11 UPD:更正了在题解页面Latex代码漏出的问题。**
------------
## 该题解主要由以下内容组成：
- **扩展欧几里得算法（exgcd）**
- **求二元一次方程/线性同余方程特解**
- **求二元一次方程/线性同余方程通解**
- **暴力解决本题**
- **对暴力进行优化**
------------
------------
------------
### **扩展欧几里得算法（exgcd）**
扩展欧几里得算法：用于求解形如ax+by=gcd(a,b)的不定方程特解。
**求解方法&证明：** 当b=0时，可以看出gcd(a,b)=a，而
$$\left\{\begin{aligned}x & = 1\\y & = 0\\\end{aligned}\right.$$
是方程的一组特解。
当b$\neq$0时，递归求解exgcd(b,a%b,x,y)，设
$$\left\{\begin{aligned}a' & = b\\b' & = a\%b\\\end{aligned}\right.$$
可以求得a'x'+b'y'=gcd(a,b)的一组特解，即x',y'。
因为a'x'+b'y'等价于bx'+(a%b)y'，再根据模运算的定义可得bx'+(a-a/b*b)y'
bx'+ay'-(a/b*b)y'
ay'+b(x'-a/b*y')
此时就得到了
$$\left\{\begin{aligned}x & = y'\\y & = x'-a/b*y'\\\end{aligned}\right.$$
**代码：**
```cpp
void exgcd(int a,int b,int &x,int &y)
{
if(!b)
{
x=1;
y=0;
return;
}
int p;
exgcd(b,a%b,x,y);
p=x;
x=y;
y=p-(a/b)*y;
return;
}
```
------------
### **求二元一次方程/线性同余方程特解**
**二元一次方程：**
对于形如ax+by=c的二元一次方程，根据扩展欧几里得算法的定义可得，当且仅当gcd(a,b)|c(|为整除)时，存在整数解。
设g=gcd(a,b),a'=a/g,b'=b/g,c'=c/g,则ax+by=c$\Leftrightarrow$a'x+b'y=c'，
此时gcd(a',b')=1,可以利用exgcd求出a'x'+b'y'=1的一组特解，继而得出
$$\left\{\begin{aligned}x_0 & = c'x'\\y_0 & = c'y'\\\end{aligned}\right.$$
我们便求得了ax+by=c的一组特解。
**线性同余方程：**
对于形如ax$\equiv$c(mod b)的线性同余方程，根据模运算的定义，在方程左侧添加一个by不会对结果造成影响，其实质就等价于ax+by=c的不定方程，利用上述方法求解便可。
------------
### **求二元一次方程/线性同余方程通解**
**二元一次方程：**
对于形如ax+by=c的二元一次方程，整数解通解为
$$\left\{\begin{aligned}x_t & = x_0+b't\\y_t & = y_0-a't\\\end{aligned}\right.$$
(t$\in\mathbb{Z}$)
对其进行证明：a$x_{t}$+b$y_{t}$
a($x_{0}$+b't)+b($y_{0}$-a't)
a$x_{0}$+ab't+b$y_{0}$+a'bt
a$x_{0}$+b$y_{0}$+at$*$
$\displaystyle\frac{b}{g}$-bt$*$
$\displaystyle\frac{a}{g}$
a$x_{0}$+b$y_{0}$+$\displaystyle\frac{abt}{g}$-$\displaystyle\frac{abt}{g}$
a$x_{0}$+b$y_{0}$
因为lcm(a,b)=$\displaystyle\frac{ab}{g}$，x,y的系数分别为a,b，所以x,y的最小增减幅度分别为$\displaystyle\frac{b}{g}$，$\displaystyle\frac{a}{g}$即b',a'。
**线性同余方程：**
同理于二元一次方程。
------------
### **暴力解决本题**
观察本题，发现对于给出的条件有三种情况：无整数解、无正整数解但有整数解、有正整数解。
**无整数解：**
要求输出-1。根据扩展欧几里得算法的定义可得，当c%gcd(a,b)$\ne$0时方程无整数解，输出-1即可。
**判断有无正整数解：**
根据二元一次方程的通解可得，x,y的增减性相反，所以若存在正整数解，则当x取得最小正整数解时，y取得最大正整数解，若此时y$\le$0,则说明该方程无正整数解。
**无正整数解但有整数解：**
要求输出x,y的最小正整数解。通过整除运算计算出x取得最小正整数解时的k值，继而求出x的最小正整数解，同理可以求出y的最小正整数解。
**有整数解：**
此时要求最复杂，要求输出正整数解的数量、x,y的最小正整数解、x,y的最大正整数解。因为当x取得最小正整数解时，y取得最大正整数解，可以利用上述方法求出x取得最小正整数解时的k值，x的最小正整数解与y的最大正整数解，该方法同样可以用于求解x的最大正整数解与y的最小正整数解。正整数解的数量=$\displaystyle\frac{x_{max}-x_{min}}{b'}$+1或$\displaystyle\frac{y_{max}-y_{min}}{a'}$+1。
**代码:**
```cpp
#include<cstdio>
using namespace std;
int read()
{
int res=0;
char c;
c=getchar();
while(c<'0'||c>'9')
c=getchar();
while(c>='0'&&c<='9')
{
res=res*10+c-48;
c=getchar();
}
return res;
}//因为题目描述“读入输出量较大，注意使用较快的读入输出方式”，所以使用快读，实际上不用也可以通过本题。
long long gcd(long long n,long long m)
{
return (n%m==0)?m:gcd(m,n%m);
}//gcd，不做说明。
void exgcd(long long a,long long b,long long &x,long long &y)
{
if(!b)
{
x=1;
y=0;
return;
}
long long p;
exgcd(b,a%b,x,y);
p=x;
x=y;
y=p-(a/b)*y;
return;
}//exgcd，本题核心代码，在上述内容中已做说明 。
int t;
int main()
{
t=read();
while(t--)
{
long long x=0,y=0,a,b,c,g,xin,yin,xax,yax,npa=0,k;//t,x,y,a,b,c,g同先前描述，xin,yin为x,y的最小正整数解,xax,yax为x,y的最大正整数解，npa为正整数解数。
a=read();
b=read();
c=read();//读入部分。
g=gcd(a,b);
if(c%g!=0)
printf("-1\n");//判断有无整数解。
else
{
a/=g;
b/=g;
c/=g;
exgcd(a,b,x,y);
x*=c;
y*=c;//求出一组特解。
for(int i=1;i<=2;i++)//因为在求解y取得最小正整数解时的k值可能使x的值小于等于0，所以要再求解一遍以求出x的最小正整数解与y的最大正整数解。
if(x<=0)
{
k=-(x/b)+1;
x+=k*b;
y-=k*a;
xin=x;
yax=y;
}//求出x取得最小正整数解时的k值，x的最小正整数解与y的最大正整数解。
else if(y<=0)
{
k=-(y/a)+1;
x-=k*b;
y+=k*a;
yin=y;
xax=x;
}//求出y取得最小正整数解时的k值，x的最大正整数解与y的最小正整数解。
if(x>0&&y>0)//判断方程有无正整数解。
{
if(x%b!=0)
{
xin=x%b;
yax=y+a*(x/b);
}
else
{
xin=x%b+b;
yax=y+a*(x/b-1);
}
if(y%a!=0)
{
yin=y%a;
xax=x+b*(y/a);
}
else
{
yin=y%a+a;
xax=x+b*(y/a-1);
}//用于求解x,y均大于等于0时x,y的最小最大正整数解。
npa=(xax-xin)/b+1;//求出正整数解的数量。
}
if(!npa)
printf("%lld %lld\n",xin,yin);
else printf("%lld %lld %lld %lld %lld\n",npa,xin,yin,xax,yax);//输出部分。
}
}
return 0;
}
```
------------
### **对暴力进行优化**
由于以上代码不优化显得又长又慢常数巨大~~而且显得丑~~，所以尝试对暴力进行优化。
运用大眼观察法，可以发现我们能用模运算直接求出x,y的最小正整数解而不用求k值，对于x,y的最大正整数解则可以直接代入二元一次方程求解，即
$$\left\{\begin{aligned}x_{max}=\displaystyle\frac{c-x_{min}}{a}\\y_{max}=\displaystyle\frac{c-y_{min}}{b}\\\end{aligned}\right.$$
**代码：**
```cpp
#include<cstdio>
using namespace std;
int read()
{
int res=0;
char c;
c=getchar();
while(c<'0'||c>'9')
c=getchar();
while(c>='0'&&c<='9')
{
res=res*10+c-48;
c=getchar();
}
return res;
}
long long gcd(long long n,long long m)
{
return (n%m==0)?m:gcd(m,n%m);
}
void exgcd(long long a,long long b,long long &x,long long &y)
{
if(!b)
{
x=1;
y=0;
return;
}
long long p;
exgcd(b,a%b,x,y);
p=x;
x=y;
y=p-(a/b)*y;
return;
}
int t;
int main()
{
t=read();
while(t--)
{
long long x=0,y=0,a,b,c,g,xin,yin,xax,yax,npa=0,k;
a=read();
b=read();
c=read();
g=gcd(a,b);
if(c%g!=0)
printf("-1\n");
else
{
a/=g;
b/=g;
c/=g;
exgcd(a,b,x,y);
x*=c;
y*=c;
xin=x>0&&x%b!=0?x%b:x%b+b;//若x>0且x%b!=0，x的最小正整数解就等于x%b，否则则需要在x%b的基础上加b。
yax=(c-xin*a)/b;//求出对应的y的最大正整数解。
yin=y>0&&y%a!=0?y%a:y%a+a;
xax=(c-yin*b)/a;//同理于x。
if(xax>0)//判断有无正整数解。
npa=(xax-xin)/b+1;//求出正整数解的数量。
if(!npa)
printf("%lld %lld\n",xin,yin);
else printf("%lld %lld %lld %lld %lld\n",npa,xin,yin,xax,yax);
}
}
return 0;
}
```
再开个O2，交一发上去，发现跑得还挺快，能拿[rank6](https://cdn.luogu.com.cn/upload/pic/1024.png)(主要是因为人少)。
最后一点：**记开long long**。
~~为什么我会这么不熟......~~