CRT&EXCRT学习笔记

学\(EXCRT\)的时间挺久了,有点忘了.
写一篇博客记录一下.

---

CRT

首先,我们要知道中国剩余定理是什么
它是用来求解这样一个同余方程的
\[ x\equiv a_1(mod\ p_1)\\ x\equiv a_2(mod\ p_2)\\ ...\\ x\equiv a_n(mod\ p_n) \]
其中,\(p_1,p_2...p_n\)两两互质.
我们先规定一些数.\(M=\prod_{i=1}^np_i,L=lcm(p_1,p_2...p_n)\)
\(CRT\)的流程是这样的.
假设现在执行到\(i\)
令\(m=\frac{M}{p_i}\)
找到一个\(x\),使得\(mx\equiv 1(mod\ p_i)\)
找\(x\)可以用逆元的方法实现.
那么,显然地,\(a_imx\equiv a_i(mod\ p_i)\)
我们令\(t_i=a_imx\)
那么,答案就是\((\sum_{i=1}^nt_i)\% L\)
这个证明显然吧

---

EXCRT

有的时候,出题人并不保证所有的\(p_i\)互质.
那么我们该怎么办呢?
其实也不难.
假设有\(2\)个方程\(x\equiv c_1(mod\ p_1),x\equiv c_2(mod\ p_2)\)
我们可以列出相应的方程.
\(x=c_1+k_1p_1,x=c_2+k_2p_2\)
那么\(c_1+k_1p_1=c_2+k_2p_2\)
\(k_1p_1=c_2-c_1+k_2p_2\)
令\(g=gcd(p_1,p_2)\),那么\(\frac{p_1}{g}k_1=\frac{c_2-c_1}{g}+\frac{p_2}{g}k_2\).
因此只有\(g\)能整除\(c_2-c_1\)时才可以合并,否则无解.
此时,可以化简方程为\(\frac{p_1}{g}k_1\equiv\frac{c_2-c_1}{g}(mod\ \frac{p_2}{g})\).
由于\((\frac{p_1}{g},\frac{p_2}{g})=1\),我们可以用\(exgcd\)求出\(\frac{p_1}{g}\)在模\(\frac{p_2}{g}\)意义下的逆元,假设为\(inv\).
那么\(k_1\equiv\frac{c_2-c_1}{g}*inv(mod\ \frac{p_2}{g})\).
我们把这个方程回代到第一个方程中,得到
\(x\equiv c_1+(\frac{c_2-c_1}{g}*inv\% \frac{p_2}{g})(mod\ \frac{p_1*p_2}{g})\)
但是,这个方程中所有的量都是已知的,因此我们就一直合并下去即可.

代码如下

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<vector>
#define N (100010)
#define inf (0x7f7f7f7f)
#define rg register int
#define Label puts("NAIVE")
#define spa print(' ')
#define ent print('\n')
#define rand() (((rand())<<(15))^(rand()))
typedef long double ld;
typedef long long LL;
typedef unsigned long long ull;
using namespace std;
inline char read(){
    static const int IN_LEN=1000000;
    static char buf[IN_LEN],*s,*t;
    return (s==t?t=(s=buf)+fread(buf,1,IN_LEN,stdin),(s==t?-1:*s++):*s++);
}
template<class T>
inline void read(T &x){
    static bool iosig;
    static char c;
    for(iosig=false,c=read();!isdigit(c);c=read()){
        if(c=='-')iosig=true;
        if(c==-1)return;
    }
    for(x=0;isdigit(c);c=read())x=((x+(x<<2))<<1)+(c^'0');
    if(iosig)x=-x;
}
inline char readchar(){
    static char c;
    for(c=read();!isalpha(c);c=read())
    if(c==-1)return 0;
    return c;
}
const int OUT_LEN = 10000000;
char obuf[OUT_LEN],*ooh=obuf;
inline void print(char c) {
    if(ooh==obuf+OUT_LEN)fwrite(obuf,1,OUT_LEN,stdout),ooh=obuf;
    *ooh++=c;
}
template<class T>
inline void print(T x){
    static int buf[30],cnt;
    if(x==0)print('0');
    else{
        if(x<0)print('-'),x=-x;
        for(cnt=0;x;x/=10)buf[++cnt]=x%10+48;
        while(cnt)print((char)buf[cnt--]);
    }
}
inline void flush(){fwrite(obuf,1,ooh-obuf,stdout);}
int n; LL c,m,c2,mo;
LL mult(LL a,LL b,LL mod){//a*b
    LL res=0,fu=1;
    if(a<0)fu=-fu,a=-a;
    if(b<0)fu=-fu,b=-b;
    while(b){
        if(b&1)res=(res+a)%mod;
        a=(a+a)%mod,b>>=1;
    }
    res*=fu;
    if(res<0)(res+=((-res-1)/mod+1)*mod);
    return res;
}
void exgcd(LL a,LL b,LL &x,LL &y){
    if(!b){x=1,y=0;return;}
    exgcd(b,a%b,y,x),y-=(a/b)*x;
}
LL gcd(LL a,LL b){return b?gcd(b,a%b):a;}
LL inv(LL a,LL b){
    LL x,y;
    exgcd(a,b,x,y);
    return (x<=0)?(x+b):x;
}
int main(){
    read(n),read(m),read(c);
    for(int i=1;i<n;i++){
        read(mo),read(c2);
        LL t=gcd(m,mo),s=inv(m/t,mo/t),tc=c,tm=m;
        m=(mo/t*tm),c=(tc+mult(tm,mult(s,(c2-c)/t,(mo/t)),m))%m;
    }
    if(c<0)c+=((-c-1)/m+1)*m;
    printf("%lld\n",c%m);
}

转载于:https://www.cnblogs.com/Romeolong/p/10076261.html