中国剩余定理(CRT)：求解模线性方程组

定义

 设有同余式组 (S)：

(S)：⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪x≡a1(modm1)x≡a2(modm2)x≡a3(modm3)⋮x≡ak(modmk)

即

x≡ai(modmi),i=1,2,..,k

 其中 ai 和 mi 已知，且当i≠j 时， 有 (mi,mj)=1 , 求x对模 m1m2...mk 的解。

解的情况

 有唯一解， 证明略。
 下面会构造出此解，然后证明其正确性和唯一性。
 设

M=∏ki=1mi//即所有m的积

Mi=Mmi//即除mi以外所有m的积

//M−1i为Mi的逆元(对模mi)，即MiM−1i≡1(modmi)

 则(S)的解为:

x=tM+∑i=1kaiMiM−1i

即

x=(∑i=1kaiMiM−1i)(modM)

 下证其正确性:

我们考查任一在[1, k]内的第i项: aiMiM−1i, 有

aiMiM−1i≡{ai⋅1≡(modmj),j=i0≡(modmj),j≠i

这里不难理解， 第二个为0是因为mj|Mi。
所以x满足如下条件:
∀i∈{1,2,..,k}, 有

x=aiMiM−1i+∑j≠i(ajMjM−1J)≡ai+∑j≠i0≡ai(modmi)

综上， 此解正确。

 下面证明其唯一性

假设x1和 x2都是方程组(S)的解，那么：
∀i∈{1,2,..,k}, 有 x1−x2≡0(modmi)
即 mi|(x1−x2)
那么显然M|(x1−x2)
即x1与x2之间相差M的整数倍，那么在模M的范围内，只有一解。

更通用的解法

 前面说的模线性方程组是k个m互素的情况， 那么更一般的呢？
 更一般的模线性方程也能求解，只不过需要从前往后遍历每一个式子， 每一个都与下一个式子合并，作为一个新的式子继续与下一个合并， 最后得出结果。

 重点！！！
 我们先假设k=2时的情况，然后再慢慢推广。
 此时(S)为

(S)：{x≡a1(modm1)(1)x≡a2(modm2)(2)

即

{x=k1m1+a1(3)x=k2m2+a2(4)

 我们将(3)(4)联立，得

k1m1=k2m2+a2−a1

, 即

k1m1≡a2−a1(modm2)(5)

//这个式子特别重要

 因为(5)式未知数只有k1，所以我们用扩展gcd可以求出k1的值， 不妨设求出的k1=t, 则对某t’, 有k1=t+t′m2。我们将它代入(3)式，得到

x≡(a1+tm1)(modlcm(m1m2))(6)//lcm为最小公倍数

 那么(6)式中的x即为同时满足(1)和(2)的x值，至此k=2的情况讨论完毕。
 至于k>2的情况， 我们可以继续将(6)式与(3)式像(1)式和(2)式那样执行上面的操作，这样到最后得到的x就满足所有的条件了。
 具体的操作可参考代码理解， 这里简单说一下中间量的维护：
 每次对2个式子进行以上操作时， 第二个式子和以前一样，不作讨论，而第一个式了如果要写成 x≡a1(modm1)的形式的话，a1和m1都是要变的，那么变成什么样子决定着上述算法的最终实现。

 首先来说m1的变化 ， 从(1)到(6)， m1变成了lcm(m1,m2), 如果接着往下走的话， m1又会变成lcm(m1,m2,m3),…, 直到最后， m1变成所有m的最小公倍数。
 而对于a1的变化 ，我们可以这样想，如果k=1的话， 结果ans是等于a1的，而k = 2的时候， 观察(6)我们会发现， a1变成了a1+tm1, 那结果呢， 结果ans(不看模)也变成了a1+tm1, 如果再往下写你会发现， ans 不是变成了a1+tm1， 而是变成了ans+t⋅lcm, 其中t 是由我们解(5)得到的。
 这里说的有点乱，具体的细节可以结合下面的实现代码去理解。

实现

 m[i] 与m[j] 互质:

int CRT(int a[], int m[], int k) {    int i, d, x, y, Mi, ans = 0, M = 1;    for (i = 0; i < k; i++) M *= m[i];     for (i = 0; i < k; i++) {        Mi = M / m[i];        d = extgcd(m[i], Mi, x, y);     // y 为逆元 -- Mi*y === 1 (% m[i])        ans = (ans + a[i]*y*Mi) % M;    }    if (ans > 0) return ans;    else return (ans + M);}

 m[i] 与m[j] 不一定互质:

LL mod(LL a, LL m) {    //处理取模    LL res = a%m;    if(res <= 0) res += m;    return res;}bool CRT(LL a[], LL m[], LL k, LL& ans) {    ans = a[0];    LL lcm = m[0], x, y, d;    if(ans == 0) ans = m[0];    for(LL i = 1; i < k; ++i) {        d = extgcd(lcm, m[i], x, y);    //求t: t = (a[i]-ans)/d*x;        if((a[i]-ans)%d) return 0;        ans = mod(ans + lcm*mod((a[i]-ans)/d*x, m[i]/d),(lcm/d*m[i]));        lcm = lcm/d*m[i];    }    return 1;}

例题

POJ - 1006 Biorhythms
HDU - 1573X问题
POJ - 2891 Strange Way to Express Integers

平生写过的最浮夸的博客。。