裴蜀定理与扩展欧几里得

裴蜀定理与扩展欧几里得

裴蜀定理

定理

对于任意整数\(a,b\)，设\(d=\gcd(a,b)\)，则方程\(ax+by=d(s\in \mathbb{Z})\)必定有无数组整数解\((x,y)\)。

且对于\(d\)的任意整数倍也有解。

证明

设\(d=\gcd(a,b),s\)为\(ax+by\)的最小正值，我们需要证明\(d=s\)。

设\(r=b\bmod s\)，则\(r=b-qs(q=\left\lfloor\frac{b}{s}\right\rfloor\in \mathbb{Z})\)

展开得\(r=b-q(ax+by)=a(-qx)+b(1-qy)\)，这时令\(r=ax'+by',x'=-qx,y'=1-qy\)，则\(r\)也是\(ax+by\)的一个值。

根据模运算的性质有：\(0\le r<s\)，由于\(s\)是\(ax+by\)的最小正值，则\(r=0\)，故\(s|b\)。

类似的，

设\(r=a\bmod s\)，则\(r=a-qs(q=\left\lfloor\frac{a}{s}\right\rfloor\in \mathbb{Z})\)

展开得\(r=a-q(ax+by)=a(1-qx)+b(-qy)\)，这时令\(r=ax'+by',x'=1-qx,y'=-qy\)，则\(r\)也是\(ax+by\)的一个值。

根据模运算的性质有：\(0\le r<s\)，由于\(s\)是\(ax+by\)的最小正值，则\(r=0\)，故\(s|a\)。

\(s|a,s|b\implies s|\gcd(a,b),s|d\)。

再者，设\(a=md,b=nd\)。则\(ax+by=d(mx+ny),\gcd(m,n)=1\)。

由于\(s=ax+by=d(mx+ny)\ge 1\)，所以可以得到\(d|s\)。

故有\(s|d,d|s\implies d=s\)。

\(QED.\)

裴蜀定理可以扩展到\(n\)元形式。

也即\(n\)元一次方程\(\sum_{i=1}^na_ix_i=k\)的充要条件是\(\gcd(a_1,a_2…a_n)|k\)。

证明过程模仿\(2\)元形式的证明，应用数学归纳法即可。

应用

裴蜀定理

如上文所说，只需要打一个GCD即可。

#include<iostream>
#include<cstdio>
using namespace std;
#define int long long
int gcd(int a,int b){
    return b==0?a:gcd(b,a%b);
}
signed main(){
    int n;cin>>n; 
    int a,ans;cin>>a;ans=a;
    for(int i=2;i<=n;i++){
        cin>>a;ans=gcd(ans,a);
        if(ans<0)ans=-ans;
    }
    cout<<ans<<"\n";
}

向量

首先，减去\((a,b)\)等价于加上\((-a,-b)\)，其余同理，所以只需要考虑四个向量：\((a,b),(b,a),(-a,b),(-b,a)\)。

所以可以得到方程组：

\[\left\{ \begin{aligned} (A-C)a+(B-D)b&=x \\ (B+D)a+(A+C)b&=y \\ \end{aligned} \right. \]

设\(\gcd(a,b)=d\)，则两个方程有解的充要条件是\(d|x,d|y\)（裴蜀定理）。此时即得到了\(A-C,A+C,B-D,B+D\)的值，不妨记为\(f_1,f_2,f_3,f_4\)。则：\(A=\frac{f_1+f_2}{2},C=\frac{f_2-f_1}{2},B=\frac{f_3+f_4}{2},D=\frac{f_4-f_3}{2}\)。

所以说，\(f_1,f_2\)需要奇偶性相同，\(f_3,f_4\)同理，现在考虑如何保证奇偶性相同呢？不妨分类讨论：

1. \(f_1,f_2,f_3,f_4\)均为偶数，此时有：\(x=f_1a+f_3b\)，提取一个公约数\(2\)，变成：\(d|\frac{x}{2}\implies 2d|x\)，同理可以得到：\(2d|y\)。

\(2d|x,2d|y\)

2. \(f_1,f_2,f_3,f_4\)均为奇数，此时有：\(x+a+b=(f_1+1)a+(f_3+1)b\)，再提取一个公约数\(2\)，变成\(2d|(x+a+b),2d|(y+a+b)\)。

\(2d|(x+a+b),2d|(y+a+b)\)

3. \(f_1,f_2\)为偶数，\(f_3,f_4\)为奇数，此时有：\(x+b=f_1a+(f_3+1)b,y+a=(f_4+1)a+f_2b\)，类似的我们可以得到：\(2d|(x+b),2d|(y+a)\)。

\(2d|(x+b),2d|(y+a)\)

4. \(f_1,f_2\)为奇数，\(f_3,f_4\)为偶数，此时有：\(x+a=(f_1+1)a+f_3b,y+b=f_4a+(f_2+1)b\)，类似的我们可以得到：\(2d|(x+a),2d|(y+b)\)。

\(2d|(x+a),2d|(y+b)\)

四种情况满足一种即可，前提是\(d|x，d|y\)。

#include<iostream>
#include<cstdio>
using namespace std;
#define int long long//警钟长鸣
int gcd(int a,int b){
	return b==0?a:gcd(b,a%b);
}
signed main(){
	int n,a,b,x,y;cin>>n;
	while(n--){
		cin>>a>>b>>x>>y;
		int d=gcd(a,b);
		d<<=1;
		if(x%d==0&&y%d==0)cout<<"Y\n";
		else if((x+a)%d==0&&(y+b)%d==0)cout<<"Y\n";
		else if((x+b)%d==0&&(y+a)%d==0)cout<<"Y\n";
		else if((x+a+b)%d==0&&(y+a+b)%d==0)cout<<"Y\n";
		else cout<<"N\n";
	}
	return 0;
}

模

这个题是真的得手玩了。

引理：一个数\(m\)在\(k\)进制下的数位和可以表示为\(m+x(1-k)\)。

证明：

由于一个数\(m\)总可以表示为\(\sum_{i=0}^{\infty}x_ik^i\)，\((0\le x_i<k)\)

分类讨论：

1. \(m<k\)，此时\(x=0\)
2. \(m=k^i,(i\in\mathbb{Z})\)

此时\(k\)进制下数位和\(=1=x+(1-x)=x+(1-k^i)\)，而由于当\(k=1\)时，\((1-k^i)=0\)，所以根据因式分解的知识，其必然存在一个因式为\((1-k)\)。此时说明\((1-k^i)\)也符合引理。

3. \(m\in \mathbb{Z}\)

尝试对\(2.\)中的方法进行扩展。设其有\(n\)位，则\(m=\sum_{i=0}^nx_ik^i,\text{数位和为}\sum_{i=0}^nx_i\)。

接着，将其拆开\(\sum_{i=0}^nx_i=\sum_{i=0}^nx_i·1=\sum_{i=0}^nx_i(k^i+(1-k^i))=\sum_{i=0}^nx_ik^i+\sum_{i=0}^nx_i(1-k^i)=m+t(1-k)\)

得证。

\(QED.\)

这样的话，我们就只需要判断方程：\(an+t(1-k)\equiv c(\bmod b)\)是否有解。

也即\(ax+(1-k)y+bz=c\)有无解，裴蜀定理即可。

#include<bits/stdc++.h>
using namespace std;
#define int long long//注意
int gcd(int a,int b){
	return b==0?a:gcd(b,a%b);
}
signed main(){
	int a,b,c,k,T;cin>>T;
    while(T--){cin>>a>>b>>c>>k;
		puts(c%(gcd(gcd(a,b),1-k))?"No":"Yes");	
	}
}

扩展欧几里得算法(exgcd)

应用：求解\(ax+by=\gcd(a,b)\)的不定方程，并给出一组解。

算法思想

根据普通\(gcd\)，我们知道\(\gcd(a,b)=\gcd(b,a\bmod b)\)。

这就意味着\(ax+by=\gcd(a,b)\implies bx+(a\bmod b)y=\gcd(a,b)\)有解。

这启发我们运用这个过程去解决方程。

显然我们可以根据\(\gcd\)一直化简到方程\(a'x+0y=a\)。这是欧几里得算法的终末态。

此时显然有一组解：\(x=1,y=0\)。(\(y=0\)是为了防爆精度)。

然后考虑回推。

设\(a'x+b'y=\gcd(a,b)\)是由\(ax+by=\gcd(a,b)\)推出，这个不定方程的解为\(x_0,y_0\)。则显然\(a'=b,b'=a\bmod b=a-\left\lfloor\frac{a}{b}\right\rfloor b\)。

带回式子，可以得到\(bx_0+(a-\left\lfloor\frac{a}{b}\right\rfloor b)y_0=\gcd(a,b)\)

整理一下，可以得到\(ay_0+b\left(x_0-\left\lfloor\frac{a}{b}\right\rfloor y_0\right)=\gcd(a,b)\)。

所以\(ax+by=\gcd(a,b)\)的解\(x_1=y_0,y_1=x_0-\left\lfloor\frac{a}{b}\right\rfloor y_0\)。

显然，求解\(\gcd(a,b)\)是一个“递”的过程，得出不定方程解是一个“归”的过程，可以考虑使用递归方法实现。

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,x,y);
	int z=x;x=y,y=z-(a/b)*y;
	return d;
} 

那么回归正题：求解\(ax+by=t\)的方程。

根据裴蜀定理，其有解的充分必要条件是\(\gcd(a,b)|t\)。现在我们来考虑求解它。

方法一

先解出\(ax_0+by_0=\gcd(a,b)\)，然后两式同乘\(\frac{t}{\gcd(a,b)}\)，即可得到\(ax+by=t\)的解\(x_1=x_0\frac{t}{\gcd(a,b)},y_1=y_0\frac{t}{\gcd(a,b)}\)。

方法二

显然，对于不定方程\(ax+by=1\)，有解的充要条件是\(\gcd(a,b)=1\)，对于\(ax+by=t\)的方程，设\(d=\gcd(a,b)\)，将方程化为\(\frac{a}{d}x+\frac{b}{d}y=\frac{t}{d}\)，解出方程\(\frac{a}{d}x+\frac{b}{d}y=1\)，对答案乘\(\frac{t}{d}\)即可。

下面提供方法一的代码：

void solve(int a,int b,int t,int &x,int &y){  //ax+by=t
    int d=exgcd(a,b,x,y);
    if(t%d){  
        puts("无解");return ;
    }
    x*=t/d,y*=t/d;
}

无穷解的讨论

设\(ax+by=t\)有解，且一组解为\(x_0,y_0\)，下面证明它有无穷多组解。

设\(a(x_0+\Delta x)+b(y_0-\Delta y)=t(\Delta x,\Delta y\in \mathbb{Z})\)。

拆开得到\(a\Delta x=b\Delta y\implies \frac{a}{b}=\frac{\Delta y}{\Delta x}=\frac{m\gcd(a,b)}{n\gcd(a,b)}=\frac{m}{n}\)。

显然，满足这个关系式的最小整数对\((\Delta x,\Delta y)\)为\((n,m)\)。且其任意整数倍都是符合条件的。

故可以得到：

\[\left\{ \begin{aligned} \Delta x=k\frac{b}{\gcd(a,b)}\\ \Delta y=k\frac{a}{\gcd(a,b)}\\ \end{aligned} \right. \]

其中\(k\)为任意整数。

那么借助它，我们可以得到一些方法：

对于不定方程\(ax+by=t\)中，满足\(x\)为最小正值，则将任意解中的\(x_0\)对\(\frac{b}{\gcd(a,b)}\)取模(记

得防负数)即可。同时就可以求出对应的\(y\)值。

同余方程

还有一个应用就是解同余方程\(ax\equiv t(\bmod b)\)

这个式子可以转化为求解不定方程\(ax-by=t\)

解出这个式子，所有的\(x\)都是合法的，无解说明不存在。

应用

青蛙的约会

根据题意，设\(p\)为跳的次数，则容易得到同余方程：

\(x+mp\equiv y+np(\bmod L)\)，移项变式得到：

\(p(m-n)\equiv y-x(\bmod L)\)，等价于求解\(p(m-n)+qL=y-x\)。

取\(q\)的最小正值解即可。

#include<iostream>
#include<cstdio>
using namespace std;
#define int long long
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,x,y);
	int z=x;x=y,y=z-(a/b)*y;
	return d;
} 
signed main(){
	int a,b,m,n,l,x,y;cin>>x>>y>>m>>n>>l;
	if(n<m){
		swap(n,m);swap(x,y);
	}
	int d=exgcd(n-m,l,a,b);
	int p=x-y;
	if(p%d){
		cout<<"Impossible\n";
	}
	else {
		a=(a*(p/d)%(l/d)+(l/d))%(l/d);cout<<a<<"\n";
	}
}

最大公约数问题

描述：求不定方程\(ax+by=\gcd(a,b)\)的一组解\((x_0,y_0)\)，满足\(|x_0|+|y_0|\)最小。

显然，由于\(\min(a,b)\ge \gcd(a,b)\ge 1\)，故这组解中\(x,y\)其中一个就是最小正值，把两种情况算出取\(\min\)即可。

#include<iostream>
#include<cstdio>
using namespace std;
inline int abs(int x){
    return (x<0?-x:x);
}
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,x,y);
	int z=x;x=y,y=z-(a/b)*y;
	return d;
}
int main(){
	int t,a,b,x,y;
	cin>>t;
	while(t--){
		cin>>a>>b;
		int d=exgcd(a,b,x,y);
		int ans=1e9;
		x=(x%b+b)%b,y=(d-a*x)/b;
		ans=x-y;int xx=x,yy=y;
		y=(y%a+a)%a,x=(d-b*y)/a;
		if(abs(xx)+abs(yy)<abs(x)+abs(y))cout<<xx<<" "<<yy<<"\n";
		else cout<<x<<" "<<y<<"\n";
	}
} 

荒岛野人

这个题注意观察到\(n\le 15,M\le 10^6\)。这启发我们直接暴力自大到小枚举答案，然后枚举每一对野人，以青蛙的约会的方法判断有无解。

注意枚举的下界就是野人居住洞穴的编号最大值。

注意不能二分/倍增，会爆掉的。（答案是不规则的）。

#include<iostream>
#include<cstdio>
using namespace std;
#define N 105
#define int long long
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,x,y);
	int z=x;x=y,y=z-(a/b)*y;
	return d;
} 
int c[N],p[N],l[N],n;
bool check(int mid){
	int x,y;
	for(int i=1;i<=n;i++){
		for(int j=i+1;j<=n;j++){
			int a=p[i]-p[j],b=mid,C=c[j]-c[i];
			int d=exgcd(a,b,x,y);
			if(C%d)continue;
			x=x*(C/d);b/=d;
			if(b<0)b=-b;
			int ti=(x%b+b)%b;
			if(ti>l[i]||ti>l[j])continue;
			return false;
		}
	}
	return true;
}
signed main(){
	cin>>n;
	int ans=0;
	for(int i=1;i<=n;i++){
		cin>>c[i]>>p[i]>>l[i];
		ans=max(ans,c[i]);
	}
	while(!check(ans))++ans;
	cout<<ans<<"\n";
}

conscious

观察数据范围\(1\le a,sub<m\)，可以发现本题就是让我们求出\(ax+ym=s\)的最小正值解\(x\)，且\(|s-sub|\)是最小的。

借助裴蜀定理，我们其实是可以确定\(|s-sub|\)的值的。

根据裴蜀定理，\(ax+ym\)可以得到所有\(\gcd(a,m)\)的整数倍。这个\(s\)也就转化了。

所以我们可以发现，可能的\(s\)只有两个：\(\left\lfloor\frac{sub}{\gcd(a,m)}\right\rfloor\gcd(a,m),\left\lceil\frac{sub}{\gcd(a,m)}\right\rceil\gcd(a,m)\)。

注意\(\left\lceil\frac{sub}{\gcd(a,m)}\right\rceil\gcd(a,m)\)有可能大于等于\(m\)，此时这个解就需要被舍弃。

这两个决策谁更优呢？设\(k=sub\bmod \gcd(a,m)\)，则:

1. \(2k<\gcd(a,m)\)，说明\(\left\lfloor\frac{sub}{\gcd(a,m)}\right\rfloor\gcd(a,m)\)更优
2. \(2k=\gcd(a,m)\)，说明\(\left\lfloor\frac{sub}{\gcd(a,m)}\right\rfloor\gcd(a,m),\left\lceil\frac{sub}{\gcd(a,m)}\right\rceil\gcd(a,m)\)一样优，此时对二者的\(x\)的最小正值取较小值。
3. \(2k>\gcd(a,m)\)，说明\(\left\lceil\frac{sub}{\gcd(a,m)}\right\rceil\gcd(a,m)\)更优。

求出最优方案中对应的\(x\)，然后直接计算答案：\(ax\bmod m\)。

注意若\(|sub -ax\bmod m|\ge sub\)，直接将答案设为\(0\)即可。

#include<bits/stdc++.h>
using namespace std;
#define int long long
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,x,y);
	int z=x;x=y,y=z-(a/b)*y;
	return d;
}
int solve(int a,int b,int k){//ax+by=k 
	int x,y;int d=exgcd(a,b,x,y);
	b/=d,x*=k/d;
	return (x%b+b)%b;  
}
signed main(){
	ios::sync_with_stdio(false);
	int T;cin>>T;
	while(T--){
		int sb,ans,res,a,m,sub,x,y;cin>>a>>m>>sub;int d=exgcd(a,m,x,y);
		int p=sub%d;
		if(sub+d-p>=m||p+p<d)
			ans=solve(a,m,sub-p);
		else if(p+p>d)
			ans=solve(a,m,sub+d-p);
		else {
			ans=min(solve(a,m,sub-p),solve(a,m,sub+d-p));
		}
		res=ans;ans=((res*a)%m+m)%m;
		if(abs(ans-sub)>=sub){
			cout<<"0 0\n";continue;
		}
		cout<<ans<<" "<<res<<"\n";
	}
	return 0;
}

AHOI2005洗牌

注意观察，手玩几组会发现一个规律：位置为\(x\)的数，交换一次就会变成\(2x\bmod (n+1)\)。

所以本题实际解决的是：\(L2^m\equiv x(\bmod n+1)\)。

进行移项变式，可以得到\(\frac{x}{2^m}+y(n+1)=L\)。求出最小正值解即可。

#include<iostream>
#include<cstdio>
using namespace std;
#define int __int128
#define ll long long
void read(int &x){
    x=0;
    char ch=getchar();
    while(ch>'9'||ch<'0')ch=getchar();
    while(ch>='0'&&ch<='9')x=x*10+(ch-'0'),ch=getchar();
}
int power(int a,int b,int p){
	int ans=1;
	while(b){
		if(b&1)ans=ans*a%p;
		a=a*a%p;
		b>>=1;
	}
	return ans%p;
}
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,x,y);
	int z=x;x=y,y=z-(a/b)*y;
	return d;
}
signed main(){
    int n,m,l;
	read(n),read(m),read(l);
	int x,y;int d=exgcd(power(2,m,n+1),n+1,x,y);
	x=(x%(n+1)+(n+1))%(n+1); 
	int inv=x;ll ans=inv%(n+1)*l%(n+1);
	cout<<ans; 
}