2018/6/10 星期日

中国剩余定理

形式化定义为,求满足下列等式的最小的正整数N,其中$k_i,r_i$由题目给定:
$N\mod m_1=r_1\Leftrightarrow N=k_1m_1+r_1$
$N\mod m_2=r_2\Leftrightarrow N=k_2m_2+r_2$
$N\mod m_3=r_3\Leftrightarrow N=k_3m_3+r_3$
$\cdots$
$N\mod m_n=r_n\Leftrightarrow N=k_nm_n+r_n$
(其实质是求一组模线性方程的解,所以有的地方也叫它 模线性方程 )

之前做过两次中国剩余定理的题目,当时也只是似懂非懂的把题目给做了, 有些没有理解的地方并没有去深究.今天跟乐乐一起讨论这个定理的时候, 我给她讲了大概的步骤, 就是通过不断把两个模线性方程合并为一个模线性方程, 到最后合并到只剩下两个方程的时候, 通过欧几里得求解模线性方程组就能得到最后的答案, 实际上前面的合并过程也是用扩展欧几里得做到的. 但是之前没有弄懂两个等式到底是怎么合并成为一个式子的,没想到今天被乐乐给讲懂了. 下面是我的理解:

对于每一个式子$N=k_im_i+r_i$,都确定了一个解系,每个解系都是一个无上界的等差数列, 那么这n个式子就确定了n个等差数列,这n个等差数列的交集肯定符合这n个等式, 所以这n个等差数列的交集中最小的正整数就是答案, 如果交集为空集, 则无解. 这是一个容易理解的思路, 但不易于用程序求解, 如果真要按照这个思路用程序求解, 可以写出类似于素数筛的代码,其时间复杂度太高.

中国剩余定理的思路是, 通过将两个模线性方程合并为一个,从而减小问题的规模, 一直进行下去, 就可以得到最后的答案, 在理解了扩展欧几里得之后, 对于这个问题,主要的难点就在于理解如何将两个式子合并为一个式子了.

就拿前面两个式子来说.
$N = k_1m_1 + r_1$
$N = k_2m_2 + r_2$
注意在每一步中都要把$N$看成一个解系而非一个数字.

1. 如果$r_1=r_2=0$ 那么合并之后: $N=k\cdot lcm(m_1,m_2)+0$
2. 如果$r_1=0,r_2\not=0$, 那么合并之后: $N=k\cdot lcm(m_1,m_2)+r_2$
3. 如果$r_2=0,r_1\not=0$, 那么合并之后: $N=k\cdot lcm(m_1,m_2)+r_1$
4. 如果$r_1\not=0, r_2\not=0$, 那么合并之后: $N=k\cdot lcm(m_1,m_2)+R$
   (其中$R$是满足这两个式子的最小的$N$,可以通过扩展欧几里得求出来)

对第一种情况进行扩展, 如果所有的$r$都为0, 那么答案就是所有$m$的最小公倍数(k取1).例如如果只有三个模线性方程,那么答案就是这三个$m$的最小公倍数(k取1).实际上前面三种情况都是第四种情况的特例,我们只需要解决第四种情况就可以了. 现在来讨论如何计算出第四种情况中的$R$.我们把前面两个式子进行变换,转换成一个标准的扩展欧几里得问题:
$N = k_1m_1 + r_1$
$N = k_2m_2 + r_2$
$\Rightarrow k_1m_1 + r_1 = k_2m_2 + r_2$
$\Rightarrow k_1m_1 - k_2m_2 = r_2 - r_1$
令$a=m_1, x=k_1, b=m_2, y=-k_2, n=r_2-r_1$, 则上式变为了标准的扩展GCD问题:
$\Rightarrow ax+by=n$
如果该方程无解, 那么原问题也无解. 否则,我们可以计算出一组解$(x_0,y_0)$满足$ax_0+by_0=n$, 此时有
$k_1=x_0+t\cdot T_1$
$-k_2=y_0+t\cdot T_2$
将$k_1$代入$N=k_1m_1+r_1$(或者将$k_2$代入$N=k_2m_2+r_2$)即可得到满足前两个式子所确定的两个等差数列的交集,即满足前面两个式子的解系, 其中最小的正整数即为我们要求的$R$, 那么如何快速求出$R$呢, 很简单, 只需要将$k_1$调整到最小的正整数,那么$R=k_1m_1+r_1$,求出了$R$之后, 就完成了前面两个式子的合并, 合并后的式子$N=k\cdot lcm(m_1,m_2)+R$中, $R$和$lcm(m_1,m_2)$都是确定的,将式子写成$N=k\cdot M + R$, 然后将该式继续与剩下的式子进行合并, 最后$R$就是答案.

代码实现如下:

/* hihoCoder1303 */

#include<bits/stdc++.h>
using namespace std;
const int maxn = 1e3+5;
typedef long long LL;

LL m[maxn], r[maxn];

LL gcd( LL a, LL b )
{
    while(a && b) a>b? a%=b : b%=a;
    return a+b;
}

LL ext_gcd(LL a, LL b, LL &x, LL &y)
{
    if(b == 0)
    {
        x = 1; y = 0;
        return a;
    }
    LL c = ext_gcd(b, a%b, x, y);
    LL x0 = x;
    x = y;
    y = x0 - a / b * y;
    return c;
}

LL chain_remain(int n)
{
    LL M = m[0], R=r[0], k1, k2;
    for(int i=1; i<n; ++i)
    {
        LL d = gcd(M, m[i]);
        LL c = r[i] - R;
        if(c % d) return -1; //无解
        ext_gcd(M/d, m[i]/d, k1, k2);
        k1 *= c/d; //原方程的解
        LL T = m[i] / d;
        k1 = (k1 % T + T) % T;  //将k1调整到最小正整数解
        R = M * k1 + R;  //R = m[1] * k1 + r[1]
        M = M / d * m[i]; // M = lcm(M, m[i]) = M * m[i] / gcd(M, m[i])
    }
    return R;
}

int main()
{
    #ifdef WFX
    freopen("in.txt","r",stdin);
    #endif

    int N;
    scanf("%d", &N);
    for(int i=0; i<N; ++i) scanf("%lld%lld", &m[i], &r[i]);
    printf("%lld", chain_remain(N));

    return 0;
}