其他分享
首页 > 其他分享> > 时隔多年,再次复习 CRT(数论全家桶)

时隔多年,再次复习 CRT(数论全家桶)

作者:互联网

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