「exgcd」学习笔记

定义

又名扩展欧几里得算法(辗转相除法)

是用来求 ax+by=gcd(a,b) 中未知数的算法


算法证明

拿到一组 a,b ,设 G=gcd(a,b)

目标:求出满足 ax+by=G1xy

如果 已知一组 x2,y2 ,满足 bx2+ (amodb)y2=G2

此时结合 (1)(2)
ax+by= bx2+ (amodb)y23

那么当 如果 满足时,目标就成了求满足 (3)x,y,其中 a,b,x2,y2 均已知

根据取模运算是 amodb=ab×ab
所以方程 (3) 实则是

ax+by=bx2+(ab×ab)y2

> ax+by=bx2+ay2b×aby2

> ax+by=ay2+b(x2aby2)

那么在根据方程 ax+by=G1 得出一组必然的解:

x=y2,y=x2aby24

可见只需求出 x2,y2 ,就能求出正确的 x,y,问题就转化成了求 x2,y2

将方程 (1)(2) 对比一下:

ax+by=G1

bx2+ (amodb)y2=G2

发现原方程的 a 变成了 b,原方程的 b 变成了 a 而已

所以新的方程也可以看做 ax+by=G 的形式,然后按照上面的程序进行下来,发现问题又变为求 x3,y3 再求 x4,y4 ...... 即一个递归的过程

递归过程中 a,b 不断被替换为 b,amodb,这个过程和普通的欧几里得是一样的,所以最后会出现 an=G,bn=0

那么这特就是最后一层,此时就要直接返回了,需要一组 xn,yn 满足 anxn+bnyn=G5,然而 an=G,bn=0,所以只要 (5)xn1 就必定满足了,yn 就随便取个 0 就行了

最后一层结束后,就开始返回,知道最上一层,每一层都可以通过下一层的 xk+1,yk+1 求出当前层的 xk,yk

整个过程就是:以辗转相除的方式向下递进,不断缩小系数,保证会出现有确定解的最后一层


exgcd板子

int exgcd(int a,int b,int &x,int &y)
{
if(b==0) {x=1;y=0;return a;}
int d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}

x,y


答案处理

对于更为一般的方程,ax+by=c,他有解当切仅当 cmodd=0(设 d=gcd(x,y))。通过 exgcd 求出一组特解 x0,y0,令他们同时 ×c/d,就得到了 ax+by=c 的一组特解 cdx0,cdy0

事实上,方程 ax+by=c 的通解可以表示为:

x=cdx0+kbd y=cdy0kad

其中 k 取遍正数集合,其余量上面已设


例题1:同余方程

题目链接 同余方程

问题处理

题目问的是满足 axmodb=1 的最小正整数(a,b均为正整数)

问题可以转化成求 ax+by=1x 的值,其中 y 是个引入的辅助数(y一定是一个负数,但是写成 + 的形式以便 exgcd 算法)

问题仍需转化,exgcd 求得是 ax+by=gcd(a,b),所以求方程 ax+by=m 必须满足的条件就是 mmodgcd(a,b)=0


简单证明一下:

由最大公因数的定义,a,b 均是 gcd(a,b) 的倍数,因为 m=ax+by,所以 m 一定是 gcd(a,b) 的倍数,即 mmodgcd(a,b)=0


可得这道题中,必须满足 1modgcd(a,b)=0,那么 gcd(a,b) 就必须等于 1 了,即 a,b 互质(数据一定是满足的)

之后通过上述的 exgcd 算法即可,但是题目要求 x 的最小值,仅仅 exgcd 无法满足,所以需要答案处理


答案处理

求出来的 x,y 仅仅满足 ax+by=1,而 x 不一定是最小正整数解。有可能太大,也可能是负数

解决方案是:x 批量地减去或加上 b,能保证 ax+by=1,因为:

ax+by=1

>ax+by+k×bak×ba=1

>a(x+kb)+(yka)b=1

1.显然这并不会把方程中的任何整数变为非整数

2.加上或减去 b 不会使 x 错过任何解,可以这么理解:


已经求出一组 x,y 使得 ax+by=1,也就是 (1ax)÷b=yy 是整数,可见目前 1axb 的倍数

现在给 x 加上 kb(k为整数),则变为:

(1a(x+kb))÷b

=(1ax)÷b+ak

仍满足其为 b 的倍数,由此当 x 的变化量为 b 的倍数时,等式仍满足


因此到最后,如果 x 太小就不断 +b 直到 >=0,反之则一直减直到最小正整数值,就是这么写

x=(x%b+b)%b;

总代码如下

#include<bits/stdc++.h>
#define endl '\n'
using namespace std;
template<typename Tp> inline void read(Tp&x)
{
x=0;register bool z=1;
register char c=getchar();
for(;c<'0'||c>'9';c=getchar()) if(c=='-') z=0;
for(;'0'<=c&&c<='9';c=getchar()) x=(x<<1)+(x<<3)+(c^48);
x=(z?x:~x+1);
}
int a,b,x,y;
int exgcd(int a,int b,int &x,int &y)
{
if(b==0) {x=1;y=0;return a;}
int d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
signed main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt","r",stdin);
freopen("out.txt","w",stdout);
#endif
read(a),read(b);
exgcd(a,b,x,y);
x=(x%b+b)%b;
cout<<x;
}

例题2:青蛙的约会

题目链接 青蛙的约会

问题处理

两点在数轴上,A,BAx,这个位置,By 这个位置,数轴长度为 lA 每次向右移 m 个长度,B 每次向右移 n 个长度,题意表明他是一个环,也就是说 Ak 下一次 Ak+1 到达 (Ak+m)modl 的位置,Bk 下一次 Bk+1 到达 (Bk+n)modl 的位置

设他跳 k 次后相遇,即 Ak=Bk

带进去就是 (x+km)modl=(y+kn)modl


左右括号内同时减(加)去一个数等式仍然成立:

如果 amodP=bmodP

(ac)modP=(bc)modP


那么原式就可以变为

>((xy)+km)modl=knmodl (左右两边括号内同时减去 y )

>((xy)+k(mn))modl=0 (左右两边同时减去 kn )

那么根据该方程,(xy)+k(mn) 一定是 l 的倍数,不妨设 (xy)+k(mn)=t×l ( t 为新设的倍数值)

>k(mn)+tl=yx (现在知道为什么设 t 了吧)

于是终于将推导式转化成了 exgcd 的形式,k(mn)+tl=gcd(mn,l) ,借鉴上一道题的证明,yx 一定是 gcd(mn,l) 的倍数

就需要思考一下数据是否满足,题目要求如果两点永远无法相遇则输出 Impossible,所以如果不满足 yxgcd(mn,l) 的倍数时,输出 Impossible 即可

exgcd 求出 k(mn)+tl=gcd(mn,l)k 之后,对他进行答案处理即可


答案处理

1.借鉴上题的答案处理,我们刚刚求出来的时第一层的结果,但他显然不一定是第一次相遇的结果,例题1答案处理的式子是 x=(xmodb+b)modb

2.我们导出来的是 k(mn)+tl=yx,而他求的是 k(mn)+tl=gcd(mn,l),满足 yxgcd(mn,l) 的倍数,所以真正的答案应该是求出的答案 ×(yx)÷gcd(mn,l),设 c=yx,d=gcd(mn,l),则答案应 ×c/d

借鉴这两点,推导出这道题的答案处理式子为

x=((x*(c/d))%(b/d)+(b/d))%(b/d);

x 表示求出来那个 kb 表示 lc,d 上面已设,之前的 b÷d


总代码如下

#include<bits/stdc++.h>
#define int long long
#define endl '\n'
using namespace std;
template<typename Tp> inline void read(Tp&x)
{
x=0;register bool z=1;
register char c=getchar();
for(;c<'0'||c>'9';c=getchar()) if(c=='-') z=0;
for(;'0'<=c&&c<='9';c=getchar()) x=(x<<1)+(x<<3)+(c^48);
x=(z?x:~x+1);
}
int x,y,n,m,l,xx,yy;
int exgcd(int a,int b,int &x,int &y)
{
if(b==0) {x=1;y=0;return a;}
int d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
signed main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt","r",stdin);
freopen("out.txt","w",stdout);
#endif
read(xx),read(yy),read(m),read(n),read(l);
int a=m-n,b=l,c=yy-xx;
if(a<0) a=-a,c=-c;
int d=exgcd(a,b,x,y);
if((c%d)!=0) cout<<"Impossible";
else cout<<(((x*(c/d))%(b/d)+(b/d))%(b/d));
}

注意事项

放到 exgcd 里的这个 a 需要保证是正的,同时 mn 对应的一定是 yx,反之,如果 (mn)<0,将他改为 nm 的同时,yx 也要改为 xy,即代码中的这一行

int a=m-n,b=l,c=yy-xx;
if(a<0) a=-a,c=-c;

例题3:倒酒

题目链接 倒酒

问题处理

问题比较简单,也不需要怎么证明,输入 a,b

显然,能够倒出的最小数为 d=gcd(a,b)

而倒 A 的次数 x 和倒 B 的次数 y 满足 byax=d

放到 exgcd 中即可


答案处理

难点在于答案处理,怎么搞出来最小解?

显然,求出的 x 为负的,所以输出的应为 x

那我们按照前两题的思路这么打

x=(x%m-m)%m;
y=(y%n+n)%n;

但是这样是错误的,Wa 90

根据我们上面打过的总的答案处理,x=cdx0+kbd y=cdy0kad,本道题的 c=d,所以 x=x0+kbd y=y0+kad

同时保证最后输出的是正的最小值,先让 xa 变成 x,a,因为我们求出来的这个 x 显然是负数

之后,在保证 x,y>0 的前提下不断地枚举这个 k,根据上面这个公式:

x<0 时,让他不断加上 bd ,当他大于 0 的一瞬间,就是正的最小值

同理,当 x>=0 时,y 做出对应调整,不断地减去 ad

如下

x*=-1;a*=-1;
while((x<0)||(y<0))
x+=(x<0)?b/d:0,
y-=(x>=0)?a/d:0;

需着重理解一下


总代码如下

#include<bits/stdc++.h>
#define int long long
#define endl '\n'
using namespace std;
template<typename Tp> inline void read(Tp&x)
{
x=0;register bool z=1;
register char c=getchar();
for(;c<'0'||c>'9';c=getchar()) if(c=='-') z=0;
for(;'0'<=c&&c<='9';c=getchar()) x=(x<<1)+(x<<3)+(c^48);
x=(z?x:~x+1);
}
int x,y,n,m;
int exgcd(int a,int b,int &x,int &y)
{
if(b==0) {x=1;y=0;return a;}
int d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
signed main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt","r",stdin);
freopen("out.txt","w",stdout);
#endif
read(n),read(m);
int d=exgcd(n,m,x,y);
x*=-1;n*=-1;
while((x<0)||(y<0))
x+=(x<0)?m/d:0,
y-=(x>=0)?n/d:0;
cout<<d<<endl<<x<<' '<<y;
}

例题4:一元二次不定方程(模板)

题目链接一元二次不定方程

问题处理

给定一式子 ax+by=c

  1. xmin,ymin
  2. xmax,ymax
  3. 求解的个数

xmin,ymin

通过 exgcd 很容易求出 x,y,在通过简单的答案处理即可求出 xmin,ymin


xmax,ymax

不难理解,当 x 为最大值时,即对应 y 的最小值,同理通过此方式可求出 y 的最大值


求解的个数

现拿到一组解 axmin+ymax=c

使 x 不断加上 b,使 y 不断减去 a,直到 axmax+bymin=c

这样就求出了不定方程的所有解,不难发现,此过程中,共出现了 xmaxxminb+1 或者说 ymaxymina+1 组解,即解的个数


答案处理

首先设 x=xmin,y=ymin,d=gcd(a,b)


处理最小值

借鉴例题2和总答案处理中的通解,xmin,ymin 的答案处理为

x=(x*(c/d)%(b/d)+(b/d))%(b/d);
y=(y*(c/d)%(a/d)+(a/d))%(a/d);

但同时,题目要求 x,y 均为正整数,所以 0 也不行,故还需要

x=(x?x:x+(b/d)),y=(y?y:y+(a/d));

处理最大值

根据我们问题处理中最大值的求法,得出最大值为

x_max=(c-y*b)/a,y_max=(c-x*a)/b;

处理解的个数

直接根据我们上面推出的式子

ans=(x_max-x)/(b/d)+1;

或者

ans=(y_max-y)/(a/d)+1;

最后注意

题目要求若不存在 x,y 均为正整数的解,只输出 xmin,ymin,只需要判断 xmax,ymax 中是否存在 0 的即可


总代码如下

#include<bits/stdc++.h>
#define int long long
#define endl '\n'
using namespace std;
template<typename Tp> inline void read(Tp&x)
{
x=0;register bool z=1;
register char c=getchar();
for(;c<'0'||c>'9';c=getchar()) if(c=='-') z=0;
for(;'0'<=c&&c<='9';c=getchar()) x=(x<<1)+(x<<3)+(c^48);
x=(z?x:~x+1);
}
int T,a,b,c,x,y;
int exgcd(int a,int b,int &x,int &y)
{
if(b==0) {x=1;y=0;return a;}
int d=exgcd(b,a%b,y,x);
y-=a/b*x;
return d;
}
signed main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt","r",stdin);
freopen("out.txt","w",stdout);
#endif
read(T);
while(T--)
{
read(a),read(b),read(c);
x=y=0;
int d=exgcd(a,b,x,y);
if(c%d)
{
cout<<-1<<endl;
continue;
}
x=((x*(c/d))%(b/d)+(b/d))%(b/d);
y=((y*(c/d))%(a/d)+(a/d))%(a/d);
x=(x?x:x+(b/d)),y=(y?y:y+(a/d));
int x_max=(c-y*b)/a,y_max=(c-x*a)/b,ans=(y_max-y)/(a/d)+1;
if(x_max<=0||y_max<=0) cout<<x<<' '<<y<<endl;
else cout<<ans<<' '<<x<<' '<<y<<' '<<(c-y*b)/a<<' '<<(c-x*a)/b<<endl;
}
}
posted @   卡布叻_周深  阅读(49)  评论(2编辑  收藏  举报
相关博文:
阅读排行:
· 微软正式发布.NET 10 Preview 1:开启下一代开发框架新篇章
· 没有源码,如何修改代码逻辑?
· NetPad:一个.NET开源、跨平台的C#编辑器
· PowerShell开发游戏 · 打蜜蜂
· 凌晨三点救火实录:Java内存泄漏的七个神坑,你至少踩过三个!
点击右上角即可分享
微信分享提示