中国剩余定理(CRT)

前置知识(逆元)

定义

对于一个线性同余方程 \(ax \equiv 1 \pmod b\),则 \(x\) 称为 \(a \bmod b\) 的逆元,记作 \(a^{-1}\)
\(a \bmod b\) 的逆元存在当且仅当 \(a\perp p\)

费马小定理求逆元

对于\(p\)为质数,由费马小定理

\[a^{p-1}\equiv 1\pmod b
\]

\[\begin{align*}
ax &\equiv a^{b-1} \pmod b\\
x &\equiv a^{b-2} \pmod b
\end{align*}
\]

于是可以用快速幂求解

il int qpow(int a,int b,int p){
    ri int as=1;
    while(b>0){
        if(b&1) as=as*a%p;
        a=a*a%p,b>>=1;
    }
    return as;
} 
il int Inv(int a,int p){
    return qpow(a,p-2,p);
}

扩展欧几里得法求逆元

递归解线性同余方程 \(ax \equiv 1 \pmod b\)即可

il void exgcd(cs int a,cs int b,int &x,int &y,int &g){
    if(!b) return x=1,y=0,g=a,void();
    return exgcd(b,a%b,y,x,g),y-=x*(a/b),void();
}
il int Inv(int a,int p){
    int x=0,y=0,g=0;
    exgcd(a,p,x,y,g);
    return (x%p+p)%p;
}

线性求1~n的逆元

\[x^{-1}\equiv
\begin{cases}
1, x=1
\\
- \lfloor \frac{p}{x} \rfloor \cdot (p\bmod x)^{-1} ,otherwise
\end{cases}
\pmod p
\]

il void solve(int n,int p){
    inv[1]=1;
    for(ri int i=2;i<=n;++i){
        inv[i]=(-p/i*inv[p%i]%p+p)%p;
    } 
    return;
}

线性求任意n给数的逆元

首先计算 \(n\) 个数的前缀积,记为 \(s_i\),然后使用快速幂或扩展欧几里得法计算 \(s_n\) 的逆元,记为 \(sv_n\)
因为 \(sv_n\)\(n\) 个数的积的逆元,所以当我们把它乘上 \(a_n\) 时,就会和 \(a_n\) 的逆元抵消,于是就得到了 \(a_1\)\(a_{n-1}\) 的积逆元,记为 \(sv_{n-1}\)
同理我们可以依次计算出所有的 \(sv_i\),于是 \(a_i^{-1}\) 就可以用 \(s_{i-1} \times sv_i\) 求得。
所以我们就在 \(O(n + \log p)\) 的时间内计算出了 \(n\) 个数的逆元。

il void solve(int n,int p){
    s[0]=1;
    for(ri int i=1;i<=n;++i){
        s[i]=s[i-1]*a[i]%p;
    } 
    sv[n]=Inv(s[i],p);
    for(ri int i=n;i;--i){
        sv[i-1]=sv[i]*a[i]%p;
    }
    for(ri int i=1;i<=n;++i){
        inv[i]=s[i-1]*sv[i]%p;
    }    
    return;
}

CRT

中国剩余定理 (Chinese Remainder Theorem, CRT) 可求解如下形式的一元线性同余方程组

\[\begin{cases} x \equiv r_1\ ({\rm mod}\ p_1) \\ x\equiv r_2\ ({\rm mod}\ p_2) \\ ... \\ x \equiv r_n\ ({\rm mod}\ p_n)\end{cases}
\]

其中 \(p_1, p_2 ... p_n\)两两互质

结论

1 . 计算所有模数的积 \(n\)

2 . 对于第 \(i\) 个方程:

  • 计算 \(m_i=\frac{n}{n_i}\)
  • 计算 \(m_i\) 在模 \(n_i\) 意义下的逆元 \(m_i^{-1}\)
  • 计算 \(c_i=m_i \cdot m_i^{-1}\)(不要对 \(n_i\) 取模)。

3 . 方程组在模 n 意义下的唯一解为:\(x=\sum_{i=1}^k a_i \cdot c_i \pmod n\)

il int crt(int a[],int n[],int k){
    int x=0,N=1;
    for(ri int i=1;i<=k;++i) N*=n[i];
    for(ri int i=1,m,inv,y;i<=k;++i){
        m=N/n[i];
        exgcd(m,n[i],inv,y);
        inv=(inv%n[i]+n[i])%n[i];
        x=(x+a[i]*m*inv%N)%N;
    }
    return (x%N+N)%N;
}
AC·code
#include<bits/stdc++.h>
#define il inline
#define ri register
#define cs const
#define int __int128
using namespace std;
cs int MAXN=15;
il int rd(){
    ri int x=0,f=1;
    ri char c=getchar();
    while(c<'0'||c>'9') {if(c=='-') f=-f;c=getchar();}
    while(c>='0'&&c<='9') x=x*10+(c^48),c=getchar();
    return x*f;
}
il void wt(int x){
    if(x<0) x=-x,putchar('-');
    if(x>9) wt(x/10);
    return putchar(x%10+48),void();
}
il void exgcd(cs int a,cs int b,int &x,int &y){
    if(!b) return x=1,y=0,void();
    exgcd(b,a%b,y,x);
    return y-=x*(a/b),void();
}
il int crt(int a[],int n[],int k){
    int x=0,N=1;
    for(ri int i=1;i<=k;++i) N*=n[i];
    for(ri int i=1,m,inv,y;i<=k;++i){
        m=N/n[i];
        exgcd(m,n[i],inv,y);
        inv=(inv%n[i]+n[i])%n[i];
        x=(x+a[i]*m*inv%N)%N;
    }
    return (x%N+N)%N;
}
signed main(){
    int n=0,a[MAXN]={0},b[MAXN]={0};
    n=rd();
    for(ri int i=1;i<=n;i++){
        a[i]=rd(),b[i]=rd();
    }
    int ans=crt(b,a,n);
    wt(ans);
    return 0;
}

exCRT

对于\(p_1, p_2 ... p_n\)不一定互质的情况就需要用到 exCRT

结论

两个方程

1 . 设两个方程分别是 \(x\equiv a_1 \pmod {m_1}、x\equiv a_2 \pmod {m_2}\)
2 . 将它们转化为不定方程:\(x=m_1 \cdot p+a_1=m_2 \cdot q+a_2\) ,其中 \(p, q\) 是整数,则有 \(m_1 \cdot p-m_2 \cdot q=a_2-a_1\)
3 . 由裴蜀定理,当 \(a_2-a_1\) 不能被 \(\gcd(m_1,m_2)\) 整除时,无解;
4 . 其他情况下,可以通过扩展欧几里得算法解出来一组可行解 \((p, q)\)
5 . 则原来的两方程组成的模方程组的解为 \(x\equiv b\pmod M\),其中 \(b=m_1 \cdot p+a_1,M=\text{lcm}(m_1, m_2)\)

多个方程

用上面的方法两两合并即可。

il int excrt(int a[],int b[],int n){
    int a1=a[1],b1=b[1];
    for(ri int i=2,x,y;i<=n;++i){
        exgcd(b1,b[i],x,y);
        if((a[i]-a1)%gcd){printf("No solution");return 0;}
        x=((x*((a[i]-a1)/gcd))%(b[i]/gcd)+(b[i]/gcd))%(b[i]/gcd);
        a1=(a1+x*b1)%(b1*b[i]/gcd),b1*=b[i]/gcd;
    }
    return a1;
}
AC·code
#include<bits/stdc++.h>
#define il inline
#define ri register
#define cs const
#define F(s) freopen(#s".in","r",stdin),freopen(#s".out","w",stdout); 
#define int __int128
using namespace std;
il int rd(){
    ri int x=0;
    ri char c=getchar();
    while(c<'0'||c>'9') c=getchar();
    while(c>='0'&&c<='9') x=x*10+(c^48),c=getchar();
    return x;
}
il void wt(int x){
    if(x<0) putchar('-'),x=-x;
    if(x>9) wt(x/10);
    return putchar(x%10+48),void();
}
int gcd;
il void exgcd(cs int a,cs int b,int &x,int &y){
    if(!b) return x=1,y=0,gcd=a,void();
    return exgcd(b,a%b,y,x),y-=a/b*x,void();
}
il int excrt(int a[],int b[],int n){
    int a1=a[1],b1=b[1];
    for(ri int i=2,x,y;i<=n;++i){
        exgcd(b1,b[i],x,y);
        if((a[i]-a1)%gcd){printf("No solution");return 0;}
        x=((x*((a[i]-a1)/gcd))%(b[i]/gcd)+(b[i]/gcd))%(b[i]/gcd);
        a1=(a1+x*b1)%(b1*b[i]/gcd),b1*=b[i]/gcd;
    }
    return a1;
}
cs int N=1e5+5;
int n,a[N],b[N];
signed main(){
    n=rd();
    for(ri int i=1;i<=n;++i) a[i]=rd(),b[i]=rd();
    wt(excrt(b,a,n));
    return 0;
}

原文链接: https://www.cnblogs.com/windymoon/p/17077270.html

欢迎关注

微信关注下方公众号,第一时间获取干货硬货;公众号内回复【pdf】免费获取数百本计算机经典书籍

    中国剩余定理(CRT)

原创文章受到原创版权保护。转载请注明出处:https://www.ccppcoding.com/archives/314830

非原创文章文中已经注明原地址,如有侵权,联系删除

关注公众号【高性能架构探索】,第一时间获取最新文章

转载文章受原作者版权保护。转载请注明原作者出处!

(0)
上一篇 2023年2月16日 下午1:31
下一篇 2023年2月16日 下午1:32

相关推荐