中国剩余定理

发表于 分类于 ，

概述:

中国剩余定理 (Chinese Remainder Theorem, aka. CRT). 解决的是这样一个问题：求同余方程组:

\[ \begin{cases} x \equiv b_1(mod\ a_1) \\ x \equiv b_2(mod\ a_2) \\ \vdots \\ x \equiv b_n(mod\ a_n) \\ \end{cases} \]

的解，其中 \(a_1, a_2, \cdots, a_n\) 两两互质.

实现:

首先设 \(M = \prod_{i=1}^{n}a_i\), \(M_i = M/a_i\), 显然 \((M_i,a_i) = 1\).

然后求 \(M_i\) 在 \(mod\ a_i\) 意义下的逆元，记为 \(t_i\). 则 \(x\) 的其中一解为:

\[ x_0 = \sum_{i=1}^{n} b_it_iM_i \]

又因为 \(M|a_1,M|a_2,\cdots,M|a_n\), 所以 \(x_0+kM\) 同样是该同余方程组的解 (\(k \in \mathbb{Z}\)). 同时 \(M\) 是满足上述性质的最小的数 (最小公倍数), 所以 \(x_0+kM\) 其实是该方程组的通解， 最小的正整数解即为 \(x_0\ mod\ M\).

1
2
3
4
5
6
7
8
9
10
11
12
13
int crt(int n) {
    int mm = 1; \\大M
    for (int i(1);i<=n;++i) {
        mm *= ai[i];
    }
    int ret = 0;
    for (int i(1);i<=n;++i) {
        int mi = mm/ai[i];
        int ti = inv(mi,ai[i]); \\inv为求逆元
        ret = (ret+(bi[i]*mi)%mm*ti)%mm;
    }
    return ret;
}

证明:

需要证明的只有一点: \(x_0\) 是否真的是原同余方程组的解.

对于累加项 \(b_jt_jM_j\), 当 \(i \neq j\) 时，\(M_j|a_i\), 所以 \(b_jt_jM_j \equiv 0(mod\ a_i)\). 当 \(i = j\) 时，\(M_jt_j \equiv 1(mod\ a_i)\), \(b_jt_jM_j \equiv b_j(mod\ a_i)\). 所以:

\[ \begin{aligned} \sum_{j=1}^n b_jt_jM_j & \equiv \sum_{j=1,j \neq i}^{n}b_jt_jM_j + b_it_iM_i \\ & \equiv b_i(mod\ a_i) \end{aligned} \]

扩展:

同样由于逆元的原因，当 \(a_1,a_2,\cdots,a_n\) 不两两互质时没有 \((M_i,a_i)=1\), 原来的方法也就无效了.

假设我们已经求得了前 \(i-1\) 个方程的通解 \(x_0+kM\), 其中 \(M=[a_1,\cdots,a_{i-1}]\),\(k \in \mathbb{Z}\), 现在加上第 \(i\) 个方程 \(x \equiv b_i(mod\ a_i)\).

设存在 \(t \in \mathbb{Z}\), 使得 \(x_0+tM \equiv b_i(mod\ a_i)\), 这样这 \(i\) 个方程便同时满足了.

对上述方程进行变形:

\[ \begin{aligned} & x_0+tM \equiv b_i(mod\ a_i) \\ \iff & Mt \equiv b_i - x_0(mod\ a_i) \end{aligned} \]

以 \(t\) 为未知数，由裴蜀定理，当 \((b_i-x_0) \nmid (M,a_i)\) 时该线性同余方程无解，原方程组无解， 否则可以通过扩展欧几里得算法求出 \(t\). 此时 \(x_0+tM\) 即为这 \(i\) 个方程的特解。令 \(M' = [M,a_i]\), 则方程组的通解为 \((x_0+tM)+kM'\), 最小正整数解为 \((x_0+tM)\ mod\ M'\).

于是求解 n 个同余方程的方程组便化为了求 n 次线性同余方程.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
typedef long long LL;
typedef __int128 LLL;
LL excrt(LL a[],LL b[],int n) {
    LL mm(1),x(0);
    for(int i(1);i<=n;++i) {
        b[i] = (b[i]-x%a[i]+a[i])%a[i]; //处理负数
        LL t,tmp;
        LL g = exgcd(mm,a[i],t,tmp);
        if(b[i]%g != 0) {
            return -1;
        }
        t = (t+a[i]) % a[i];
        t = ((LLL)(t) * b[i]/g)%a[i];   //防止乘法溢出
        x = ((LLL)(x) + (LLL)(t)*mm)%lcm(mm,a[i]);
        mm = lcm(mm,a[i]);  //迭代
    }
    return x;
}