标签：复习 CRT int a1 pmod m1 m2 时隔多年 equiv

1.CRT

中国剩余定理，用来求解同余方程组

\begin{cases} x\equiv a_1\pmod {m_1} \\ x\equiv a_2\pmod {m_2} \\ x\equiv a_3\pmod {m_3} \\ …………\\ x\equiv a_n\pmod {m_n}\\ \end{cases}⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧​x≡a1​(modm1​)x≡a2​(modm2​)x≡a3​(modm3​)…………x≡an​(modmn​)​

其中，mm 两两互质。

显然，这个柿子一定有无数个解，我们只需要找到最小解即可。

伟大的古人研究出来的算法

中国剩余定理告诉我们，这个柿子存在一个通解可以计算。

x=\sum_{i=1}^n a_iM_it_i,(M_i=(\prod_{j=1}^nm_j)/mi,t_i=inv(M_i,m_i))x=i=1∑n​ai​Mi​ti​,(Mi​=(j=1∏n​mj​)/mi,ti​=inv(Mi​,mi​))

对于逆元的计算，我们可以使用扩展欧几里德算法来计算。

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

为什么，因为我们构造一个特殊的方程：

M_it_i+m_iy=\gcd(m_i,M_i)Mi​ti​+mi​y=gcd(mi​,Mi​)

得到这个方程的解,因为两两互质，显然 \gcd(m_i,M_i)=1gcd(mi​,Mi​)=1

于是方程化为了：

M_it_i+m_iy=1Mi​ti​+mi​y=1

我们要求的逆元是 inv(M_i,m_i)inv(Mi​,mi​)

于是我们化简简单移项原方程：

M_it_i=1-m_iyMi​ti​=1−mi​y

再次化简得：

M_it_i\equiv 1 \pmod {m_i}Mi​ti​≡1(modmi​)

于是我们就可以得知：

M_it_i+m_iy=\gcd(m_i,M_i)Mi​ti​+mi​y=gcd(mi​,Mi​)

的特解中的t_iti​，就是我们要求的逆元。

那么我们反观我们的定理：

\begin{cases} x\equiv a_1\pmod {m_1} \\ x\equiv a_2\pmod {m_2} \\ x\equiv a_3\pmod {m_3} \\ …………\\ x\equiv a_n\pmod {m_n}\\ \end{cases}⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧​x≡a1​(modm1​)x≡a2​(modm2​)x≡a3​(modm3​)…………x≡an​(modmn​)​

x=\sum_{i=1}^n a_iM_it_i,(M_i=(\prod_{j=1}^nm_j)/mi,t_i=inv(M_i,m_i))x=i=1∑n​ai​Mi​ti​,(Mi​=(j=1∏n​mj​)/mi,ti​=inv(Mi​,mi​))

是不是感觉清爽很多了？

那么我们尝试证明一下。

证明：

因为 M_i=(\prod_{j=1}^nm_j)/miMi​=(∏j=1n​mj​)/mi ,即除了 m_imi​ 以外所有模数的倍数。

所以 \forall j,a_iM_it_i\equiv0\pmod {m_k}∀j,ai​Mi​ti​≡0(modmk​) ,因此不会对其他的方程造成影响，这很好，给我们的定理证明创造了机会。

又因为 t_i=inv(M_i,m_i)ti​=inv(Mi​,mi​),即 M_it_i\equiv1\pmod {m_i}Mi​ti​≡1(modmi​)

所以 a_iM_it_i\equiv a_i\pmod {m_i}ai​Mi​ti​≡ai​(modmi​)

因此累加每一个方程组的答案即可得到方程组的通解。

若答案过大，对 xx 取模 \prod_{j=1}^nm_j∏j=1n​mj​ 即可。

时间复杂度 O(n\log n)O(nlogn) 主要是扩展欧几里得算法和循环带来的复杂度。

附上模板代码：

#include<cstdio>
#include<algorithm>
#define int long long
#define N 1919
using namespace std;
int n,a[N],b[N];
int exgcd(int a,int b,int &x,int &y){
    if(b==0){x=1,y=0;return a;}
    int ans=exgcd(b,a%b,x,y);
    int t=x;x=y,y=t-a/b*y;
    return ans;
}
int crt(int n,int *a,int *m){
    int M=1,ans=0;
    for(int i=1;i<=n;i++)M*=m[i];
    for(int i=1;i<=n;i++){
        int mi=M/m[i],ti=0,y=0;
        exgcd(mi,m[i],ti,y);
        ans=(ans+(a[i]*mi*ti)%M)%M;
    }
    return (ans+M)%M;
}
signed main(){
    scanf("%lld",&n);
    for(int i=1;i<=n;i++)scanf("%lld%lld",&a[i],&b[i]);
    printf("%lld",crt(n,b,a));
    return 0;
}

2.扩展CRT

扩展 CRT 和 CRT 的思想完全不一样，扩展 CRT 更像是在合并同余方程打到化简的效果。

\begin{cases} x\equiv a_1\pmod {m_1} \\ x\equiv a_2\pmod {m_2} \\ x\equiv a_3\pmod {m_3} \\ …………\\ x\equiv a_n\pmod {m_n}\\ \end{cases}⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧​x≡a1​(modm1​)x≡a2​(modm2​)x≡a3​(modm3​)…………x≡an​(modmn​)​

当 mm 不再满足互质的关系时，我们应该怎么做？

观察一个二元一次同余方程：

\begin{cases} x\equiv a_1\pmod {m_1} \\ x\equiv a_2\pmod {m_2} \\ \end{cases}{x≡a1​(modm1​)x≡a2​(modm2​)​

可以写作：

\begin{cases} x=a_1+y_1m_1 \\ x=a_2+y_2m_2 \\ \end{cases}{x=a1​+y1​m1​x=a2​+y2​m2​​

整理得：

x=a_1+y_1m_1=a_2+y_2m_2x=a1​+y1​m1​=a2​+y2​m2​

继续化简：

a_1+y_1m_1=a_2+y_2m_2a1​+y1​m1​=a2​+y2​m2​

y_1m_1+zm_2=a_2-a_1 ( z=-y_2 )y1​m1​+zm2​=a2​−a1​(z=−y2​)

这个柿子在 a_2-a_1a2​−a1​ 可以通过扩展欧几里得算法得到特解 (y,z)(y,z)。

很显然，对于原方程：

m_1y \frac{a_2-a_1}{\gcd(m_1,m_2)}+m_2z \frac{a_2-a_1}{\gcd(m_1,m_2)} = \frac{a_2-a_1}{\gcd(m_1,m_2)}\gcd(m_1,m_2)=a_2-a_1m1​ygcd(m1​,m2​)a2​−a1​​+m2​zgcd(m1​,m2​)a2​−a1​​=gcd(m1​,m2​)a2​−a1​​gcd(m1​,m2​)=a2​−a1​

因为 x=a_1+y_1m_1x=a1​+y1​m1​,现在让 xx 尽可能小，其实是让 y\frac{a_2-a_1}{\gcd(m_1,m_2)}ygcd(m1​,m2​)a2​−a1​​ 尽可能小

我们现在要最小化 y\frac{a_2-a_1}{\gcd(m_1,m_2)}ygcd(m1​,m2​)a2​−a1​​ 不妨设为 tt。

t=y_1+k\frac{m_2}{\gcd(m_1,m_2)}t=y1​+kgcd(m1​,m2​)m2​​

带回发现：

x=a_1+y_1m_1+k\times lcm(m_1,m_2)x=a1​+y1​m1​+k×lcm(m1​,m2​)

要删掉后面那堆，显然只用对 tt 进行取模即可：

t=y_1+k\frac{m_2}{\gcd(m_1,m_2)}\pmod {\frac{m_2}{\gcd(m_1,m_2)}}t=y1​+kgcd(m1​,m2​)m2​​(modgcd(m1​,m2​)m2​​)

于是我们得到和合并后的方程组：

x\equiv m_1y_1+a_1\pmod {lcm(m_1,m_2)}x≡m1​y1​+a1​(modlcm(m1​,m2​))

#include<cstdio>
#include<algorithm>
#include<cstdlib>
#include<cmath>
#define int __int128
using namespace std;
long long n;
int exgcd(int a,int b,int &x,int &y){
    if(b==0){x=1,y=0;return a;}
    int ans=exgcd(b,a%b,x,y);
    int t=x;x=y,y=t-a/b*y;
    return ans;
}
int lcm(int a,int b){return a*b/__gcd(a,b);}
pair<int,int> merge(int a2,int m2,int a1,int m1){
    int x,y,c;
    int d=exgcd(m1,m2,x,y);
    c=a2-a1;
    if(c%d){puts("-1");exit(0);}
    x=x*c/d%(m2/d);
    if(x<0)x+=(m2/d);
    int md=lcm(m1,m2);
    a1=(m1*x+a1)%md;
    if(a1<0)a1+=md;
    m1=md;
    return {a1,m1};
}
signed main(){
    scanf("%lld",&n);
    int a1=0,m1=0;
    for(int i=1;i<=n;i++){
        long long a,m;
        scanf("%lld%lld",&m,&a);
        if(i==1)a1=a,m1=m;
        else{
            pair<int,int> fi=merge(a1,m1,a,m);
            a1=fi.first,m1=fi.second;
        }
    }
    printf("%lld",(long long)(a1%m1));
    return 0;
}

标签：复习,CRT,int,a1,pmod,m1,m2,时隔多年,equiv
来源： https://www.cnblogs.com/masida-hall-LZ/p/16629495.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。