题解 P5656 【【模板】二元一次不定方程(exgcd)】
2020/1/25 UPD:更正了证明过程中符号的问题。(By weapon_zipper)
2019/12/11 UPD:更正了在题解页面Latex代码漏出的问题。
该题解主要由以下内容组成:
- 扩展欧几里得算法(exgcd)
- 求二元一次方程/线性同余方程特解
- 求二元一次方程/线性同余方程通解
- 暴力解决本题
- 对暴力进行优化
扩展欧几里得算法(exgcd)
扩展欧几里得算法:用于求解形如ax+by=gcd(a,b)的不定方程特解。
求解方法&证明: 当b=0时,可以看出gcd(a,b)=a,而
是方程的一组特解。
当b
可以求得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')
此时就得到了
代码:
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
此时gcd(a',b')=1,可以利用exgcd求出a'x'+b'y'=1的一组特解,继而得出
我们便求得了ax+by=c的一组特解。
线性同余方程:
对于形如ax
求二元一次方程/线性同余方程通解
二元一次方程: 对于形如ax+by=c的二元一次方程,整数解通解为
(t
对其进行证明:a
a(
a
a
a
a
因为lcm(a,b)=
线性同余方程: 同理于二元一次方程。
暴力解决本题
观察本题,发现对于给出的条件有三种情况:无整数解、无正整数解但有整数解、有正整数解。
无整数解:
要求输出-1。根据扩展欧几里得算法的定义可得,当c%gcd(a,b)
判断有无正整数解:
根据二元一次方程的通解可得,x,y的增减性相反,所以若存在正整数解,则当x取得最小正整数解时,y取得最大正整数解,若此时y
无正整数解但有整数解: 要求输出x,y的最小正整数解。通过整除运算计算出x取得最小正整数解时的k值,继而求出x的最小正整数解,同理可以求出y的最小正整数解。
有整数解:
此时要求最复杂,要求输出正整数解的数量、x,y的最小正整数解、x,y的最大正整数解。因为当x取得最小正整数解时,y取得最大正整数解,可以利用上述方法求出x取得最小正整数解时的k值,x的最小正整数解与y的最大正整数解,该方法同样可以用于求解x的最大正整数解与y的最小正整数解。正整数解的数量=
代码:
#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的最大正整数解则可以直接代入二元一次方程求解,即
代码:
#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(主要是因为人少)。 最后一点:记开long long。
为什么我会这么不熟......