浅谈中国剩余定理（CRT）及其扩展

Sya_Resory · 2020-08-07 11:40:38 · 个人记录

中国剩余定理（CRT）是中国古代求解一次同余式组的方法，是数论中一个重要的定理。《孙子算经》中首次提到了同余方程组问题以及解法，因此在数学文献中也会将中国剩余定理称为孙子定理。

——摘自度娘百科

引入

有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？

此题源自《孙子算经》卷下第二十六题，大体意思是：有一个数，除以三余数为二，除以五余数为三，除以七余数为二，问这个数是多少。对于这个问题，我们可以依据一个歌谣解决：

三人同行七十稀，五树梅花廿一支，七子团圆正半月，除百零五便得知。

这句歌谣的意思是：把除以三的余数乘上70，除以五的余数乘上21，除以七的余数乘上15，他们的和除以105的余数就是答案。代入一开始的问题就可以得到，答案是(2 \times 70 + 3 \times 21 + 2 \times 15) \mod 105\ =\ 23。

简单的分析

那么为什么要乘上70、21和15这三个奇怪的数呢？我们通过观察发现：

可以发现，每加上一个数，就是为以除这个数余数为一的数为除数的余数贡献1，而对其他的除数没有影响！没错这就是个绕口令 又因为\operatorname{lcm}(3,5,7)=105，所以除以105的余数就是这个同余方程的最小整数解。

CRT入门

从简单到复杂

CRT能够解决的问题当然不止模数为3、5、7的问题。实际上，CRT可以求出如下同余方程的解：

我们令\forall i \in [1,n]\quad m_i=\dfrac{\operatorname{lcm}(p_1,p_2\dots,p_n)}{p_i}=\dfrac{\prod_{j=1}^n p_j}{p_i}，则答案即为\sum_{i=1}^n q_i\times m_i \times m_i^{-1}，其中m_i^{-1}为m_i在模p_i意义下的逆元。

证明？

即得易见平凡，仿照上例显然。留作习题答案略，读者自证不难。反之亦然同理，推论自然成立，略去过程QED，由上可知证毕。

可以知道，\forall i,j\in [1,n]

• 当i \ne j时有 q_j \times m_j \times m_j^{-1} \equiv 0 \pmod {p_i}

• 当i = j时有 q_j \times m_j \times m_j^{-1} \equiv q_j \pmod {p_i}

因此，\forall j \in [1,n] \sum_{i=1}^n q_i\times m_i \times m_i^{-1} \equiv q_j \pmod {p_j}，因此定理成立。

关于逆元的那些事

乘法逆元，主要用来求模p意义下\dfrac{a}{b}的值。（p一般是质数）

若a \cdot x \equiv 1 \pmod b，其中\gcd(a,b)=1，则称x为a的乘法逆元，用a^{-1}表示。

如何求逆元？

• 扩展欧几里得算法（exgcd）

  乘法逆元的本质是a \cdot x \equiv 1 \pmod b，我们就可以把这个柿子转化成这个亚子：ax+by=1，其中x,y \in \mathbb{Z}。而这个不定方程就可以用exgcd求解。关于exgcd的详细解释可以看 这道题。

  代码：

  int exgcd(int a,int b,int& x,int& y) { // exgcd模板
    if(!b) {
        x = 1,y = 0;
        return a;
    }
    int x1,y1,ans = exgcd(b,a % b,x1,y1);
    x = y1,y = x1 - a / b * y1;
    return ans;
  }
  
  int inv(int a,int b) { // 求逆元
    int x,y;
    exgcd(a,b,x,y); //exgcd
    return (x % b + b) % b; // 负数处理
  }

• 费马小定理&快速幂

  费马小定理：\forall a \in \mathbb{Z},p \in prime,gcd(a,p) = 1有a^{p - 1}\equiv 1 \pmod p

  我们注意到这个柿子的右边正好是个1，于是就可以把这个柿子进行亿些小小的变换：

  \left\{ \begin{array}{lr} ax \equiv 1 \pmod p & \\ a^{p-1} = a \cdot a^{p-2} \equiv 1 \pmod p \end{array} \right.

  于是我们可以得到:

  x \equiv a^{p-2} \pmod p

  而这个柿子可以用快速幂快速求解。

  代码:

  int quick_pow(int a,int b,int p) {
    if(!b) return 1;
    if(b == 1) return a % p;
    int rep = quick_pow(a,b >> 1,p);
    rep = rep * rep % p;
    if(b & 1) return rep * a % p;
    return rep;
  }
  
  int inv(int a,int p) {
    return quick_pow(a % p,p - 2,p);
  }

代码实现

别装了，我知道你们都在等着这个（大雾

// n表示同余方程组数，p数组表示模数，q数组表示余数
int crt() {
    int M = 1,ans = 0; // 乘积M，结果ans
    for(int i = 1;i <= n;i ++) M *= p[i];
    for(int i = 1;i <= n;i ++) { // 循环枚举1~n
        int m = M / p[i];
        int x = inv(m); // 求逆元
        ans += q[i] * m * x; //求和
    }
    return ans % M;
}

来个栗子

右转美食街栗子店炒锅

• 洛谷P1495 曹冲养猪

很模板的一道CRT，建立a_i个猪圈有b_i头猪无家可归就是指x \equiv b_i \pmod{a_i}，又因为\forall i,j\in[1,n] \gcd(a_i,a_j)=1，所以这道题可以用CRT来做。

代码：

#include <cstdio>

typedef long long ll; // 一定要记得开long long！

ll n,a[15],b[15];

ll exgcd(ll a,ll b,ll& x,ll& y) { // exgcd模板
    if(!b) {
        x = 1,y = 0;
        return a;
    }
    ll x1,y1,ans = exgcd(b,a % b,x1,y1);
    x = y1,y = x1 - a / b * y1;
    return ans;
}

ll inv(ll a,ll b) { // exgcd求逆元
    ll x,y;
    exgcd(a,b,x,y);
    return (x % b + b) % b;
}

ll crt() { // crt 模板
    ll M = 1,ans = 0;
    for(int i = 1;i <= n;i ++) M *= 1LL * a[i];
    for(int i = 1;i <= n;i ++) {
        ll m = M / a[i],x = inv(m,a[i]);
        ans += 1LL * b[i] * m * x;
    }
    return ans % M;
}

int main() {
    scanf("%d",&n);
    for(int i = 1;i <= n;i ++) scanf("%d%d",a + i,b + i);
    printf("%lld\n",crt());
    return 0;
}

• UVA756 Biorhythms

还是一道模板题，不过因为模数已知，我们可以用\mathcal{O}(1)的方式求解。

根据题意我们可以列出这样的方程式：

解得：

根据之前的求解柿子可以得出:

所以输出的答案就是(5544p + 14421e + 1288i - d + 21251) \mod 21252 + 1，其中括号内加上21251是为了避免负数，括号外加1是特判0的情况，此时应输出21252。

CRT拓展

扩展中国剩余定理（excrt）

你学会了中国剩余定理，非常高兴，于是去向\color{black}l\color{red}ndjy大佬吹牛皮。这时\color{black}l\color{red}ndjy大佬向你抛来了一个同余方程组：

（q_1,q_2,q_3\dots,q_n不两两互质）

这时中国剩余定理不管用了，因为当i,j\in [1,n],i \ne j时 q_j \times m_j \times m_j^{-1} \mod {p_i}不一定为0！

怎么办呢？多半是废了

欢迎我们的嘉宾——扩展中国剩余定理（excrt）！别问，问就是作者中二病又发作了

我们考虑分开求解然后合并，假设求解到了第i个方程，前i-1个方程求解的结果是ans，那么可以得到：

为了简便我们设a_1=ans,b_1=M,a_2=q_i,b_2=p_i，那么：

我们利用扩展欧几里得求解b_1x'-b_2y'=\gcd(b_1,b_2)，那么就可以化简得到：

这部分推柿子可能有些难理解，建议自己手推一遍。

例题：洛谷P4777 【模板】扩展中国剩余定理（EXCRT）

代码：

#include <cstdio>

typedef long long ll;

ll n,ans,p[100005],q[100005];

inline ll read() { // 快读
#define gc c = getchar()
    ll d = 0;
    int f = 0,gc;
    while(c < 48 || c > 57) f |= (c == '-'),gc;
    while(c > 47 && c < 58) d = (d << 1LL) + (d << 3LL) + (c ^ 48),gc;
#undef gc
    return f ? -d : d;
}

ll exgcd(ll a,ll b,ll& x1,ll& y1) { // exgcd模板
    if(!b) {
        x1 = 1,y1 = 0;
        return a;
    }
    ll x2,y2,ans = exgcd(b,a % b,x2,y2);
    x1 = y2;
    y1 = x2 - a / b * y2;
    return ans;
}

ll pmul(ll a,ll b,ll mod) { // 快速乘 不用的话会爆long long
    ll rep = 0;
    while(b) {
        if(b & 1) rep = (rep + a) % mod;
        a = (a << 1) % mod,b >>= 1;
    }
    return rep;
}

bool excrt(ll& ans) { // 返回值表示成功与否
    ll M = p[1];
    ans = q[1];
    for(int i = 2;i <= n;i ++) {
        ll k1,x,y;
        ll g = exgcd(M,p[i],x,y); // 扩欧
        ll c = (q[i] - ans % p[i] + p[i]) % p[i];
        if(c % g != 0) return false; // 不能整除，无解
        k1 = pmul(x,c / g,p[i] / g); // 见上面的公式 这里加了一些取模操作
        ans = M * k1 + ans;
        M = M / g * p[i];
        ans = (ans % M + M) % M;
    }
    ans = (ans % M + M) % M;
    return true;
}

int main() {
    n = read();
    for(int i = 1;i <= n;i ++) p[i] = read(),q[i] = read();
    excrt(ans);
    printf("%lld\n",ans);
    return 0;
}

扩展卢卡斯定理（exlucas）

你又学会了扩展中国剩余定理，非常高兴，于是又去向\color{black}l\color{red}ndjy大佬吹牛皮。这时\color{black}l\color{red}ndjy大佬向你抛来了一个柿子：

求这个柿子的值。

这个问题可以用到扩展卢卡斯定理（exlucas）。这个算法的思路大概是这样的：先把模数p分解质因数：

然后再列出一个同余方程组：

然后我们就可以对形如C_n^m \equiv q \pmod{p^k} \quad (p \in prime)的同余式分别求解，然后用CRT合并即可。求解单个方程式的过程比较繁琐，这里就不再赘述。

参考资料

• crt,excrt学习总结
• 中国剩余定理(CRT) & 扩展中国剩余定理(ExCRT)总结
• 中国剩余定理学习笔记
• 乘法逆元