CRT，中国剩余定理

发表于 更新于 阅读次数： 本文字数： 3.6k 阅读时长 ≈ 6 分钟

介绍

中国剩余定理 (Chinese remainder theorem, CRT)，简单而言，给定一个未知数 xx 和 nn 个关于 xx 的线性同余方程 (Congruence Equation) 组成的线性同余方程组 (Congruence Equation System)，其中任意两模数互质，求满足该方程组的最小非负整数解。

形如：

{x≡b1 (mod a1)x≡b2 (mod a2)...x≡bn (mod an)\begin{cases} x \equiv b_1\ ({\rm mod}\ a_1) \\ x\equiv b_2\ ({\rm mod}\ a_2) \\ ... \\ x \equiv b_n\ ({\rm mod}\ a_n)\end{cases}

的方程组称为线性同余方程组。

其中 ≡\equiv 为同余符号。如果给定一个正整数 aia_i （即模数），如果两个整数 xx 和 bib_i 满足 x−bix-b_i 能被 aia_i 整除，即 (x−bi) mod  ai=0(x-b_i)\rm \bmod\ a_i=0 ，那么就称整数 xx 与 bib_i 对模 aia_i 同余，记作 x≡bi ( mod  ai)x \equiv b_i\ ({\rm \bmod}\ a_i) 。

而 P1495则是 CRT 的模板题，题面较长就不放出来了，但本质就是上面的问题。

原理

原理这种东西，又怎么扯得清楚呢？

我们知道，在 CRT 中，有一个特别的限制：其中任意两模数互质，我们对问题的求解就是需要这个条件。（不信的看看扩中吧）

我们先就 P1495 给的样例来看吧，简单来说，这个问题就是在问：

找到一个最小正整数 xx ，使得：

{x mod  3=1x mod  5=1x mod  7=2\begin{cases} x\rm \bmod\ 3=1 \\ x\rm \bmod\ 5=1 \\ x \rm \bmod\ 7=2\end{cases}

以上方程组也可以转化为同余方程组的形式，不过写成这样可能好理解一些。

一下子看三个方程有一点晕吧？没关系，我们试着化大为小：

尝试找出这样的 x1,x2,x3x_1,x_2,x_3，分别使得：

x1 mod 3=1,x1 mod 5=0,x1 mod 7=0,x_1\rm \bmod3=1,x_1\rm \bmod5=0,x_1\rm \bmod7=0,

x2 mod 3=0,x2 mod 5=1,x2 mod 7=0,x_2\rm \bmod3=0,x_2\rm \bmod5=1,x_2\rm \bmod7=0,

x3 mod 3=0,x3 mod 5=0,x3 mod 7=2x_3\rm \bmod3=0,x_3\rm \bmod5=0,x_3\rm \bmod7=2

不难发现，x=x1+x2+x3x=x_1+x_2+x_3。

当然，我们不妨再进行问题的分解。

对于 x3 mod 7=2x_3\rm \bmod7=2 ，我们不难发现它复杂在于不像 x1,x2x_1,x_2 的方程那样，余数为 11 。那既然这样，我们就让它余数为 11 ！

于是转化为：

x3′ mod 3=0,x3′ mod 5=0,x3′ mod 7=1x'_3\rm \bmod3=0,x'_3\rm \bmod5=0,x'_3 \rm \bmod7=1

只需求出 x3′x'_3 ，再乘二即可。

那么怎么求 x3′x'_3 呢？对于 x3′x'_3 ，因为它一定是 33 和 55 的倍数，故设 x3′=15zx'_3=15z ，那么就有：

15z≡1( mod 7)15z \equiv 1(\rm \bmod 7)

这里要求 zz 的话需要对乘法逆元有所掌握，那么让我手算就有一点刁难我了，所以待会儿我直接用 exgcd 来解这个同余方程，详见下文！

当然当然，对于最后的结果需要为最小正整数解，只需要求：

x mod Mx \rm \bmod M

MM 为模数之积（其实为最小公倍数）。

大家可以自己手算一下。

推导到一般求法就很简单了，还不懂的可以对照下面的 “过程” 再看看。

过程

  1. 求模数的最小公倍数 MM ，由其中任意两模数互质得：

    M=∏ai[1≤i≤n]M = \prod a_i [1\le i \le n]

  2. 设数 ki=M/aik_i=M/a_i 。

  3. 找到一个最小正整数数 pip_i 使得：

    (piki) mod  ai=1(p_i k_i) \rm \bmod\ a_i = 1

  4. 则:

    x=(∑pikibi[1≤i≤n]) mod  Mx=(\sum p_ik_ib_i[1\le i\le n])\rm \bmod\ M

  5. 注意 3. 应如何找出满足条件的 pip_i ，如果直接由小到大枚举肯定时间复杂度过大，此时应观察该方程的形式，发现其可以写成：

    piki≡1(mod ai)p_ik_i\equiv 1(\rm mod\ a_i )

    的形式。熟悉的同学可以看出来：这不就是求乘法逆元吗？~~ 好吧，其实熟悉的同学在第三步就直接在求乘法逆元了。~~ 对于求乘法逆元，这里不再详细讲解，详见这篇文章。我们这里使用扩展欧几里得算法 (Extended Euclidean algorithm,ExGcd)，算法如下：

    上面的方程可以写成：

    kipi+aiy=1k_ip_i+a_iy=1

    其中 ki,aik_i,a_i 为常数，pip_i 就是需要求的未知数，yy 为辅助解。

    对于求解以上不定方程的代码模板如下：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    /*求同余方程ax≡1(\b\bmod b)的最小正整数解。*/
    #include<bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    int x,y;//存储解
    void exgcd(int a,int b){
        if(b==0){
            x=1,y=0;
            return;
        }//b=0时，方程有一组特解x=1,y=0，作为递归边界
        exgcd(b,a%b);//递归
        int tmp=x;
        x=y;
        y=tmp-(a/b)*y;
        return;
    }//辗转相除，与gcd写法类似
    int main()
    {
        int a,b;
        cin>>a>>b;
        exgcd(a,b);
        cout<<(x+b)%b;//这里不排除为负数解或非最小正整数解
    }

  6. 代码实现如下：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    #include<bits/stdc++.h>
    using namespace std;
    typedef long long ll;
    ll x,y;//存储解
    void exgcd(ll a,ll b){
        if(b==0){
            x=1,y=0;
            return;
        }//b=0时，方程有一组特解x=1,y=0，作为递归边界
        exgcd(b,a%b);//递归
        ll tmp=x;
        x=y;
        y=tmp-(a/b)*y;
        return;
    }//辗转相除，与gcd写法类似
    int main()
    {
        ll n;
        cin>>n;
        ll a[n],b[n];
        ll M=1;
        for(ll i=0;i<n;i++){
            cin>>a[i]>>b[i];
            M*=a[i];//求M
        }//输入不解释
        ll k[n];
        ll ans=0;
        for(ll i=0;i<n;i++) k[i]=M/a[i];//求ki
        for(ll i=0;i<n;i++){
            exgcd(k[i],a[i]);//扩欧解方程
            x=(x+a[i])%a[i];//这里不排除为负数解或非最小正整数解
            ans+=x*k[i]*b[i]%M;//求ans
        }
        cout<<ans%M;
    }