1. 程式人生 > >中國剩餘定理(CRT)及其拓展(ExCRT)

中國剩餘定理(CRT)及其拓展(ExCRT)

中國剩餘定理 CRT

推導

給定\(n\)個同餘方程
\[ \left\{ \begin{aligned} x &\equiv a_1 \pmod{m_1} \\ x &\equiv a_2 \pmod{m_2} \\ &... \\ x &\equiv a_n \pmod{m_n} \end{aligned} \right. \]
\(m_1, m_2 , ... , m_n\)兩兩互質

令\(M = \prod_{i=1}^{n} m_i\),求\(x \mod M\)

解決該問題的方法是構造。

我們假定最終答案的形式是一個\(n\)個項的和式,對每個同餘方程的構造反應在對應項的係數上。

如果要對每一個項分別構造,就要求為每一項乘上一個合適的數,使得每項構造的係數對其他方程的結果沒有影響。

容易想到構造
\[ M_i = \frac{M}{m_i} \]

顯然該數僅在模\(m_i\)時不為\(0\),於是改變該項的係數將不會對其他方程造成影響。

現在我們希望該項模\(m_i\)意義下是\(a_i\),但上一次的構造殘留下了一個\(M_i\)。簡單粗暴的乘上\(M_i\)在模\(m_i\)意義下的逆元\(inv_{m_i}(M_i)\),讓該項在模\(m_i\)意義下變為\(1\),然後乘上\(a_i\)就構造出來了。

綜上,答案為
\[ \sum_{i=1}^{n} M_i inv_{m_i}(M_i) a_i \mod{M} \]

模數互質條件保證了\(M_i\)非\(0\),保證了\(inv_{m_i}(M_i)\)的存在。

實現

大部分題的所有\(m_i\)都是質數,求逆元\(QPow\)即可。

對於一般的情況,上exgcd就行。

板題:洛谷P1495 曹衝養豬

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<ctime>
#include<cstdlib>
#include<algorithm>
#include<queue>
#include<vector>
#include<map>
#include<set>
using namespace std;
typedef long long ll;

namespace ExGcd{
    ll x,y;
    ll ExGcd(ll a,ll b){
        ll ans;
        if(a%b==0){
            x=0;y=1;ans=b;
        }else{
            ans=ExGcd(b,a%b);
            ll x1=x,y1=y;
            x=y1;y=x1-a/b*y1;
        }
        return ans;
    }
    bool SolveEqu(ll a,ll b,ll c){
        ll d=ExGcd(a,b);
        if(c%d!=0) return 0;
        x*=c/d;y*=c/d;
        x=(x%b+b)%b;
        y=(c-a*x)/b;
        return 1;
    }
}
ll Inv(ll a,ll m){
    ExGcd::SolveEqu(a,m,1);
    return ExGcd::x;
}

const ll CRTN=20;
namespace CRT{
    ll N;
    ll m[CRTN],a[CRTN];
    ll Sol(){
        ll ans=0,M=1;
        for(ll i=1;i<=N;i++) M*=m[i];
        for(ll i=1;i<=N;i++){
            ll Mi=M/m[i];
            ans=(ans+Mi*Inv(Mi,m[i])*a[i])%M;
        }
        return ans;
    }
}
int main(){
    using namespace CRT;
    scanf("%lld",&N);
    for(ll i=1;i<=N;i++)
        scanf("%lld%lld",&m[i],&a[i]);
    printf("%lld",Sol());
    return 0;
}

拓展中國剩餘定理 ExCRT

ExCRT和CRT並沒有什麼關係,正如ExLucas和Lucas沒什麼關係一樣(

其實從純推理的角度來看,ExCRT可能還要好想一點(

推導

問題同CRT,但是模數是任意的,並不要求互質。

這時,我們就不能保證存在逆元了。那麼如何解決該問題呢?

考慮如何合併兩個方程。如果我們找到了合併的方法,就能如法炮製將\(n\)個方程依次合併起來,得到答案。
\[ \left\{ \begin{aligned} x &\equiv a_1 \pmod{m_1} \\ x &\equiv a_2 \pmod{m_2} \end{aligned} \right. \]
去掉同餘,化為不定方程
\[ \left\{ \begin{aligned} x &= m_1 y_1 + a_1 \\ x &= m_2 y_2 + a_2 \end{aligned} \right. \]

於是得到
\[ m_1 y_1 + a_1 = m_2 y_2 + a_2 \]

只要找到一組滿足該式的\(y_1\)和\(y_2\),就能反算出\(x\),實現合併。

而我們得到的是一個二元一次的不定方程,可以用exgcd求解。

化為標準式
\[ m_1 y_1 - m_2 y_2 = a_2 - a_1 \]

解就是了。如果沒有解,說明同餘方程組無解。

於是最後化得的合併式為
\[ x \equiv m_1 y_1 + a_1 \pmod{lcm(m_1,m_2)} \]

實現

唯一需要注意的地方是,本來解方程應該解\((m_1,-m_2,a_2-a_1)\),但ExGCD不好處理負數,所以把\(m_2\)改成了\(m_1\)。因為我們並不需要用到\(y_2\),所以顯然是沒有影響的。

板題:poj2891 Strange Way to Express Integers or 洛谷P4777 擴充套件中國剩餘定理(EXCRT)

會被卡乘法爆ll,懶得改

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<ctime>
#include<cstdlib>
#include<algorithm>
#include<queue>
#include<vector>
#include<map>
#include<set>
using namespace std;
typedef long long ll;

namespace ExGcd{
    ll x,y;
    ll ExGcd(ll a,ll b){
        ll ans;
        if(a%b==0){
            x=0;y=1;ans=b;
        }else{
            ans=ExGcd(b,a%b);
            ll x1=x,y1=y;
            x=y1;y=x1-a/b*y1;
        }
        return ans;
    }
    bool SolveEqu(ll a,ll b,ll c){
        ll d=ExGcd(a,b);
        if(c%d!=0) return 0;
        x*=c/d;y*=c/d;
        x=(x%b+b)%b;
        y=(c-a*x)/b;
        return 1;
    }
}
ll Gcd(ll a,ll b){
    if(a%b==0) return b;
    return Gcd(b,a%b);
}
namespace ExCRT{
    ll a1,m1;
    void Init(){
        a1=0;m1=1;
    }
    void Expand(ll a2,ll m2){
        ExGcd::SolveEqu(m1,m2,a2-a1);
        ll y1=ExGcd::x;
        ll mn=m1*m2/Gcd(m1,m2);
        a1=(m1*y1+a1)%mn;
        m1=mn;
    }
}
int main(){
    ll N;scanf("%lld",&N);
    ExCRT::Init();
    for(ll i=1;i<=N;i++){
        ll a,m;scanf("%lld%lld",&m,&a);
        ExCRT::Expand(a,m);
    }
    printf("%lld",ExCRT::a1);
    return 0;
}