CRT&EXCRT 中國剩餘定理及其擴展

前言:

中國剩餘定理又名孫子定理。因孫子二字歧義，常以段子形式普遍流傳。

 

中國剩餘定理並非很好理解，我也理解了不少次。

CRT 中國剩餘定理

中國剩餘定理，就是一個解同餘方程組的算法。

求知足n個條件的最小的x。

看起來很麻煩。

先找一個特殊狀況：$m_1,m_2,...m_n$兩兩互質。

這個時候，構造$M=m_1*m_2*...m_n$;

令$M_i=M/m_i$;

因此，構造$n$個數，其中第$i$個數是除$i$以外的其餘全部數的倍數，而且第$i$個數$mod m_i =1$

即：$M_i x = 1 ( mod m_i ) $求出這樣一個x，就求出了 這個數。

由於$m$之間兩兩互質，因此對於$n$個這樣的方程，$x$本質上就是$M_i$在$m_i$意義下的乘法逆元。

（不會$exgcd$？左轉：EXGCD 擴展歐幾里得）

由於互質，必定有解的。

用擴展歐幾里得算就能夠。

同理，構造$n$個數。$b_1,b_2....b_n$

其中，$b_i=M_i \times x_i$

那麼，由於$b_i  = 1 (mod m_i)$，因此$ b_i * a_i = a_i (mod m_i)$

那麼，原題目中的這個x就是：$x=(a_1\times b_1+a_2\times b_2+...+a_n\times b_n) $驗證一下，是否是？

 

總得來講，

對於$mi$互質的狀況，$x=\sum_1^n M_i\times a_i\times inv_i$

其中，$inv_i$表示，$M_i$在$mod\space m_i$意義下的逆元。

固然，爲了保證$x$最小，要讓$x$和$lcm$作一些處理。固然，由於互質，因此$lcm$就是$M$了

x=(x%M+M)%M

例題：poj1006 生理週期 Biorhythms

在CRT上的小小變形。注意開始計算的時間d就行了。上述最小的$x$就是$n+d$

答案$n=(a_1\times b_1+a_2\times b_2+a_3\times b_3-d)%lcm$

#include<cstdio> #include<iostream> #include<algorithm> #include<cstdlib>
using namespace std; int p,e,l; int p1,e1,l1; int tt,d; void exgcd(int a,int b,int &x,int &y){ if(b==0){ x=1,y=0;return; } exgcd(b,a%b,y,x); y-=(a/b)*x; } int main(){ exgcd(23*28,33,l1,tt); l1=(l1%33+33)%33; l1*=23*28; exgcd(28*33,23,p1,tt); p1=(p1%23+23)%23; p1*=28*33; exgcd(23*33,28,e1,tt); e1=(e1%28+28)%28; e1*=23*33; int lcm=23*28*33; int cnt=0; while(1){ ++cnt; scanf("%d%d%d%d",&p,&e,&l,&d); if(p==-1) break; int op=(p1*p+e1*e+l1*l-d+lcm)%lcm; if(op==0) op=lcm; printf("Case %d: the next triple peak occurs in %d days.\n",cnt,op); } return 0; }

poj1006

 

EXCRT 擴展中國剩餘定理

可是，不是全部的方程，$m$都是互質的。

$m$不是互質的時候，咱們的$M_i x = 1(mod m_i) $可能就沒有解了。因此掛掉。

孫子解決不了了。

可是現代人不是孫子也解決了。

（你是孫子嗎）

如今，$m_1,m_2,...m_n$之間沒有任何關係。

只考慮兩個怎麼處理？

能夠獲得：$x=a_1+k_1*m_1 ; x=a_2+k_2*m_2$

因此 $a_1+k_1*m_1 = a_2+k_2*m_2$

$k_2*m_2-k_1*m_1=a_1-a_2$;

很像：$a\times x +b\times y=c$

 

設$m_1,m_2$ 的$gcd$ 爲 $g$

設$a_1-a_2=c$

當$c$不是$g$的倍數的時候，那就完了。（$exgcd$無解狀況）

若是是，就用$exgcd$求出$k_2\times m_2+k_1\times m_1=gcd(m_1,m_2)$的$k1$

由於$c$是$g$的倍數，因此，兩邊同時乘上$c/g$，即$k_1$乘上$c/g$

就獲得了$k_2\times m_2+k_1\times m_1=c$的解$k_1$。

固然，最好$k_1$ 再模一下$m_2$ （$k_1$,$k_2$作出調整）,防止爆$long long$

 

而後能夠反推x，

可是注意，咱們列的原方程是:$k_2\times m_2-k_1\times m_1=c$

差一個符號，因此$k_1$ 其實是 $-k_1$

$x=-k_1\times m_1+a_1$;

這樣就求出了$x$。（能夠把這個$x0$轉化成最小的非負數解）

這個$x$符合第一個方程，也符合第二個方程。設這個$x$爲$x_0$

因此，能夠獲得通解是：$x=x_0+k\times lcm(m_1,m_2)$

知足這個條件的x就知足第一第二兩個方程。知足第一第二兩個方程的全部的解也都是這個方程的解。

因此第一第二個方程和這個方程是等價的。

將這個方程轉化一下，能夠獲得新的同餘方程：$x=x_0 (mod lcm(m_1,m_2))$

這樣，咱們成功的把兩個方程轉化成了一個方程，以此類推。

最後留下的這個方程，它的x_0的最小非負數解，就是咱們要的最終答案！！！！！！！！！！！！！！

例題：poj2891

這是一個模板題，直接上代碼：

#include<cstdio> #include<algorithm> #include<iostream>
using namespace std; typedef long long ll; const int N=100000+10; int n; ll a[N],r[N]; ll exgcd(ll a,ll b,ll &x,ll &y){ if(b==0){ x=1,y=0;return a; } ll ret=exgcd(b,a%b,y,x); y-=(a/b)*x; return ret; } ll excrt(){ ll M=a[1],R=r[1],x,y,d; for(int i=2;i<=n;i++){ d=exgcd(M,a[i],x,y); if((R-r[i])%d) return -1; x=(R-r[i])/d*x%a[i]; R-=M*x; M=M/d*a[i]; R%=M; } return (R%M+M)%M; } int main() { while(scanf("%d",&n)!=EOF){ for(int i=1;i<=n;i++){ scanf("%lld%lld",&a[i],&r[i]); } printf("%lld\n",excrt()); } return 0; }

poj 2891

例題：NOI2018屠龍勇士

1.exgcd轉化成同餘方程的通常形式

2.求解同餘方程

細節：

1.Pi=1的狀況特殊處理

2.開始的ai、ki能夠不用mod p[i]以避免出現0的狀況較爲麻煩

代碼：

#include<bits/stdc++.h>
#define reg register int
#define il inline
#define fi first
#define se second
#define mk(a,b) make_pair(a,b)
#define numb (ch^'0')
#define int long long
using namespace std; typedef long long ll; template<class T>il void rd(T &x){ char ch;x=0;bool fl=false; while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true); for(x=numb;isdigit(ch=getchar());x=x*10+numb); (fl==true)&&(x=-x); } template<class T>il void output(T x){if(x/10)output(x/10);putchar(x%10+'0');} template<class T>il void ot(T x){if(x<0) putchar('-'),x=-x;output(x);putchar(' ');} template<class T>il void prt(T a[],int st,int nd){for(reg i=st;i<=nd;++i) ot(a[i]);putchar('\n');} namespace Miracle{ const int N=1e5+5; int n,m; ll a[N],p[N],st[N],k[N]; multiset<int>s; multiset<int>::iterator it; ll ad(ll x,ll y,ll mod){ return x+y>=mod?x+y-mod:x+y; } ll qk(ll x,ll y,ll mod){ ll ret=0; while(y){ if(y&1) ret=ad(ret,x,mod); x=ad(x,x,mod); y>>=1; } return ret; } void clear(){ s.clear(); } ll exgcd(ll a,ll b,ll &x,ll &y){ if(!b){ x=1;y=0;return a; } ll ret=exgcd(b,a%b,y,x); y-=(a/b)*x; return ret; } ll wrk(){//warning!!! x!=0 // cout<<" kk "<<endl; // prt(k,1,n);
    for(reg i=1;i<=n;++i){//warning!!! // a[i]%=p[i]; // k[i]%=p[i]; // if(!k[i]) return -1;
        ll x=0,y=0; ll g=exgcd(k[i],p[i],x,y); x=(x%(p[i]/g)+(p[i]/g))%(p[i]/g); if(a[i]%g) return -1; x=qk(x,a[i]/g,p[i]/g); p[i]=p[i]/g; a[i]=x; } // cout<<" after 1 "<<endl; // prt(a,1,n); // prt(p,1,n);
    for(reg i=1;i<n;++i){ // cout<<" turn "<<i<<endl;
        ll p1=p[i],p2=p[i+1]; // cout<<" p1 "<<p1<<" p2 "<<p2<<" a1 "<<a[i]<<" a2 "<<a[i+1]<<endl;
        ll x=0,y=0; ll g=exgcd(p1,p2,x,y); ll tmp=a[i+1]-a[i]; ll mo=p1/g*p2;//lcm
        tmp=(tmp%mo+mo)%mo; // cout<<" mo "<<mo<<" tmp "<<tmp<<" x "<<x<<" y "<<y<<endl;
        if(tmp%g) return -1; x=(x%(mo/p1)+(mo/p1))%(mo/p1); x=qk(x,tmp/g,mo/p1); ll nw=(a[i]+qk(x,p1,mo))%mo; a[i+1]=nw; p[i+1]=mo; } if(!a[n]) a[n]+=p[n]; return a[n]; }//warning!!! x!=0
int main(){ int t; rd(t); while(t--){ clear(); rd(n);rd(m); bool fl=true; for(reg i=1;i<=n;++i) rd(a[i]); for(reg i=1;i<=n;++i) { rd(p[i]); if(p[i]!=1) fl=false; } for(reg i=1;i<=n;++i) rd(st[i]); ll x; for(reg i=1;i<=m;++i){ rd(x);s.insert(x); } for(reg i=1;i<=n;++i){ it=s.upper_bound(a[i]); if(it!=s.begin()){ --it; k[i]=*it; s.erase(it); }else{ it=s.begin(); k[i]=*it; s.erase(it); } s.insert(st[i]); } if(fl){ ll ans=0; for(reg i=1;i<=n;++i){ ans=max(ans,(a[i]+k[i]-1)/k[i]); } printf("%lld\n",ans); }else{ printf("%lld\n",wrk()); } } return 0; } } signed main(){ Miracle::main(); return 0; } /* Author: *Miracle* Date: 2019/4/2 19:09:28 */

View Code