浅记线性同余方程（组）

技术分享 4个月前 (10-31) 0 999+

关注

线性同余方程

定义

线性同余方程就是形如 (axequiv bpmod m) 其中 (a,b,m) 是给定的整数。

解法

由同余的性质可知 (mmid ax-b) 即 (ax-b=km) 其中 (kin mathbb{Z})。

如果我们设 (k=-y) 的话，就有 (ax+my=b)，发现了吗？其实这就是 Bézout 定理。

由Bézout 定理我们可以得到，这个同余方程有解当且仅当 (gcd(a,m)mid b)。

我们考虑在有解的情况下使用扩展欧几里得算法先求解出 (ax+my=gcd(a,m)) 的一组特解 (left{begin{matrix}x=x_0 \y=y_0end{matrix}right.)，然后呢，我们就可以得到 (x=frac{x_0times b}{gcd(a,m)}) 就是原方程的一组解。

关于扩展欧几里得算法的说明

内容

我们考虑不定方程 (ax+by=gcd(a,b)) 的一组特解，我们可以采用递归的方法来求解，实际上这也就是扩展欧几里得算法：

显然当 (b=0) 时，有 (left{begin{matrix}x=1\y=0end{matrix}right.) 满足条件。
当 (bne 0) 时，我们根据欧几里得算法有 (gcd(a,b)=gcd(b,abmod b)) 于是，我们就有 ((abmod b)y+bx=gcd(b,abmod b)) 又由于 ((abmod b)y+bx=left(a-btimeslfloorfrac{a}{b}rfloorright)times y+bx) 将 (operatorname{RHS}) 展开，合并同类项后有 ((abmod b)y+bx=ay+btimes left(x-lfloorfrac{a}{b}rfloor yright)) 于是，我们令 (x_0=y,y_0=x-lfloor frac{a}{b}rfloor y) 就有 (ax_0+by_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-z*y; 	return d; }

功能介绍

以上函数的返回值为 (gcd(a,b))，注意到参数 (x,y) 均在前面加上了取地址符，表示在函数中可以改变 x 与 y 的值，而函数运行完成后 x 与 y 所保存的值就是 (ax+by=gcd(a,b)) 的一组特解。

一道模板题

洛谷 P1082：

这道题目就是模板题，方程可以写成 (ax+by=1) 的形式，于是我们使用扩展欧几里得算法，可以求出特解 (x_0) 然后 (x_0bmod b) 就是原方程的最小正整数解了。

#include<bits/stdc++.h> #define int long long using namespace std; int 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; } signed main(){ 	cin>>a>>b; 	int x,y; 	exgcd(a,b,x,y); 	cout<<(x%b+b)%b; 	return 0; }

线性同余方程组

作者太懒了，这里先讲解更加宽泛的扩展中国剩余定理吧，等以后再讲解特殊的中国剩余定理，顺便宣传一下博客：link。

问题简述

给定一个 (k) 个方程的线性同余方程组：

[left{begin{matrix} xequiv a_1pmod {m_1}\xequiv a_2pmod {m_2} \vdots \xequiv a_kpmod {m_k} end{matrix}right.]

其中 (m_1,m_2,dots,m_k) 不一定两两互质。

解题方法

我们的大致解题思路为将 (2) 个方程合并为一个新的方程，以此类推，最终我们会得到一个 (xequiv ypmod z) 的一个方程，易见上面的方程组的最小正整数解就是 (y)。

接下来我们来解决合并方程的问题，我们考虑如下两个方程：

[left{begin{matrix} xequiv a_1pmod {m_1}\xequiv a_2pmod {m_2} end{matrix}right.]

我们根据第一个式子可以写出 (x) 的通解 (x=a_1+m_1times k) 其中 (k) 为任意整数，我们将这个通解带入第二个式子就可以得到 (a_1+m_1times kequiv a_2pmod {m_2}) 我们移一下项就可以得到 (m_1times kequiv a_2-a_1pmod {m_2})，这就是上面的方程组合并后的结果。

而这个方程有解的充要条件是 (gcd(m_1,m_2)mid a_2-a_1)，这个其实就是裴蜀定理，这里不再概述。

我们继续讲，我们得到这个充要条件后我们可以判断这个方程是否有解，如果有解我们就继续进行接下来的操作。

我们设 (d=gcd(m_1,m_2))，然后将我们合并的方程变换一下就是：

[frac{m_1times k}{d}equiv frac{a_2-a_1}{d}pmod {frac{m_2}{d}} ]

然后，我们设 (m_1'=frac{m_1}{d},c=frac{a_2-a_1}{d},m_2'=frac{m_2}{d}) 于是我们就有：

[m_1'times kequiv cpmod {m_2'} ]

注意到此时 (m_1',m_2') 互质，所以 (m_1') 在模 (m_2') 的意义下存在乘法逆元，我们可以使用扩展欧几里得算法来求出逆元，即求出整数 (inv) 使得 (m_1'times invequiv 1pmod {m_2'})，所以我们继续将这个方程变换就变成了：

[kequiv ctimes invpmod {m_2'} ]

如果我们记 (k_0=ctimes inv) 则 (k) 的通解为 (k_0+m_2'times t) 其中 (t) 为任意整数。

然后我们将这个 (k) 带回一开始的式子就可以得出：

[begin{aligned} x&=a_1+m_1times(k_0+m_2'times t)\ &=(a_1+m_1times k_0)+(m_1times m_2')times t\ &=(a_1+m_1times k_0)+frac{m_1times m_2}{d}times t\ &=(a_1+m_1times k_0)+mathrm{lcm}(m_1,m_2)times tend{aligned}]

我们设 (x_0=a_1+m_1times k_0,L=mathrm{lcm}(m_1,m_2)) 所以我们就愉快地得出了：

[left{begin{matrix} xequiv a_1pmod {m_1}\xequiv a_2pmod {m_2} end{matrix}right.Longleftrightarrow xequiv x_0pmod L]

于是，我们完成了合并方程的使命！

最后其实就是一个递推的过程我们一次合并前 (2) 个方程，最后就能得到答案！

代码实现

#include<bits/stdc++.h> #define LL __int128 #define R register using namespace std; namespace fastIO{char *p1,*p2,buf[100000]; 	#define nc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,100000,stdin),p1==p2)?EOF:*p1++) 	inline void read(LL&n){LL x=0,f=1;char ch=nc();while(ch<48||ch>57){if(ch=='-'){f=-1;}ch=nc();}while(ch>=48&&ch<=57){x=(x<<3)+(x<<1)+(ch^48),ch=nc();}n=x*f;}inline void read(string&s){s="";char ch=nc();while(ch==' '||ch=='n'){ch=nc();}while(ch!=' '&&ch!='n'){s+=ch,ch=nc();}}inline void read(char&ch){ch=nc();while(ch==' '||ch=='n'){ch=nc();}}inline void write(LL x){if(x<0){putchar('-'),x=-x;}if(x>9){write(x/10);}putchar(x%10+'0');return;}inline void write(const string&s){for(R LL i=0;i<(int)s.size();i++){putchar(s[i]);}}inline void write(const char&c){putchar(c);} }using namespace fastIO; inline LL mul(LL a,LL b,const LL&mod){     a=(a%mod+mod)%mod;      b=(b%mod+mod)%mod;     LL res=0;     while(b){         if(b&1)res=(res+a)%mod;         a=(a+a)%mod;         b>>=1;     }     return res; } void exgcd(LL a,LL b,LL&x,LL&y){     if(b==0){         x=1;         y=0;     } 	else{         exgcd(b,a%b,y,x);         y-=x*(a/b);     } } LL inv_mod(LL a,LL m){     LL x,y;     exgcd(a,m,x,y);     return (x%m+m)%m; } LL gcd(LL a,LL b){     return b?gcd(b,a%b):a; } LL n,a[100005],b[100005]; signed main(){     read(n);     for(int i=0;i<n;i++){     	read(a[i]);     	read(b[i]);     }     LL a0=a[0];     LL b0=(b[0]%a0+a0)%a0;     for(int i=1;i<n;i++){         LL ai=a[i];         LL bi=(b[i]%ai+ai)%ai;         LL d=gcd(a0,ai);         LL dif=bi-b0;         LL a0_=a0/d;         LL ai_=ai/d;         LL dif_=dif/d;         LL c=(dif_%ai_+ai_)%ai_;         LL inv=inv_mod(a0_,ai_);         LL t0=mul(inv,c,ai_);         LL a0__=(a0/d)*ai;         LL mod__=a0__;         LL p=mul(a0,t0,mod__);         LL b0__=(b0+p)%mod__;         a0=mod__;         b0=b0__;     } 	write(b0);     return 0; }