中国剩余定理(CRT)及其程序实现

中国剩余定理(CRT)及其程序实现

Q:今有物不知其数,三三数之剩二,五五数之剩三,七七数之剩二,问物几何?
A:三人同行七十稀,五树梅花廿一枝,七子团圆正半月,除百零五便得知。

可能高二时候我也没想到,这个当时我特别喜欢的算法书上的定理,再次见到是在数学专业课高代上,而不是计算机算法课上。
中国剩余定理又称孙子定理,中国古代求解一次同余式组(见韩信点兵和曹冲养猪)的方法,最早可见于南北朝的数学著作《孙子算经》。

本文给出中国剩余定理的证明及程序实现。

罗栗老师讲的中国剩余定理
需要注意的是,由于整数环和多项式环性质极其相似,该定理在多项式环上也有相似结论。

Bezout定理

[latex]\large \forall p,q\in \mathbb{N} ,\exist u,v\in \mathbb{N} \; s.t. \; up+vq=(p,q)[/latex]
其中[latex](p,q)[/latex]为p,q的最大公约数。
Bezout定理通过Euclidean算法易证。

Bezout定理推论

若[latex]p,q[/latex]互质,则有[latex]up\equiv 1\pmod{q}[/latex]

中国剩余定理(Chinese Remainder Theorem)

已知[latex]m_{i}(i\in \mathbb{N} _{+})[/latex]满足[latex]\forall i,j\in \mathbb{N} _{+} ,i\not= j, (m_{i}, m_{j} )=1[/latex],
则[latex]\exist x \in \mathbb{N}, x \equiv bi \pmod{m_{i}}[/latex]。
具体地说,一定有[latex]x[/latex]:
[latex]\large x \equiv \sum b_{i} M_{i} N_{i} \pmod{M}[/latex],
其中[latex]\large M = \prod m_{i} , M_{i} = \prod _{k \not= i} m_{k}, N_{i} M_{i} \equiv 1 \pmod {m_{i}}[/latex]。

[latex]N_{i}[/latex]的存在性

由于[latex]m_{i}[/latex]两两互质,而[latex]M_{i} = \prod _{k \not= i} m_{k}[/latex],总是有[latex](m_{i} , M_{i} )=1[/latex],从而由Bezout定理:
[latex]\large \exist u_{i},N_{i}\in \mathbb{N} \; s.t. \; u_{i} m_{i} + N_{i} M_{i} = 1 \implies N_{i} M_{i} \equiv 1 \pmod{m_{i}}[/latex]

解的合理性

注意到[latex]N_{i} M_{i} \equiv 1 \pmod{m_{i}}[/latex],
而[latex]N_{i} M_{i} \equiv 0 \pmod{m_{j}} (i \not= j)[/latex],
从而满足题设方程,而[latex]\forall k \in \mathbb{Z} , x + k M_{i} \equiv x \pmod{m_{i}}[/latex],故模[latex]M[/latex]意义下解合理。

思考

中国剩余定理的构造和Langrange插值公式极其相似,按罗栗老师的话:都有类似开关的东西在控制整个式子。精妙所在。

程序实现

洛谷 P1495 【模板】中国剩余定理(CRT)/曹冲养猪
其实整个算法关键在于求[latex]M_{i}[/latex]的乘法逆元,用扩展Euclidean算法就可以了。然后再按解的构造累加起来,模一下M得到最小正整数解,为了防止溢出我每一步都模了,但是题解上看似乎最终再模也没有问题。
第一次提交,听取[latex]\color{red}{\rm WA}[/latex]声一片
中国剩余定理WA
原因写在代码注释里了,修改后完美[latex]\color{green}{\rm AC}[/latex]
中国剩余定理AC

#include <cstdio>
#include <cctype>

typedef long long int int64;

int64 n, a[20], b[20], Mi[20], M = 1, X;

inline int getint()
{
    int x = 0;
    char ch = 0;
    while (!isdigit(ch)) ch = getchar();
    while (isdigit(ch)) x = (x << 1) + (x << 3) + (ch ^ 0x30), ch = getchar();
    return x;
}

int64 exgcd(int64 a, int64 b, int64 &x, int64 &y)
//扩展Euclidean算法求解Bezout方程
{
    if (b == 0) {
        x = 1; y = 0;
        return a;
    }
    int64 r = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return r;
}

int main()
{
    n = getint();
    for (int i = 1; i <= n; i++) {
        a[i] = getint();
        b[i] = getint();
        Mi[i] = 1;
        M *= a[i];
    }

    for (int i = 1; i <= n; i++)
        Mi[i] = M / a[i];

    for (int i = 1; i <= n; i++) {
        int64 Ni, y;
        exgcd(Mi[i], a[i], Ni, y);  //求解同余方程Ni * M[i] = 1(mod a[i])
        X = (X + b[i] * Mi[i] * (Ni % a[i] + a[i])) % M;
        //Ni % a[i] + a[i] 防止Ni为负数,注意这里Ni为模a[i]意义下的乘法逆元,不能将a[i]写为M
        //写成X = ((X + b[i] * Mi[i] * Ni) % M + M) % M应该也行;
    }
    printf("%lld\n", X);

    return 0;
}

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注