//https://img2018.cnblogs.com/blog/1646268/201908/1646268-20190806114008215-138720377.jpg

浅谈中国剩余定理

中国剩余定理

定义

中国剩余定理(CRT)可以求解如下形式的一元线性同余方程组(其中 \(n_{1},n_{2},\dots,n_{k}\) 两两互质)

\[\left\{\begin{matrix} x\equiv a_{1}\pmod{n_{1}}\\ x\equiv a_{2}\pmod{n_{2}}\\ \vdots\\ x\equiv a_{k}\pmod{n_{k}} \end{matrix}\right. \]

过程

  • 计算所有模数的积 \(n\)

  • 对于第 \(i\) 个方程:

    • 计算 \(m_{i}=\frac{n}{n_{i}}\);

    • 计算 \(m_{i}\) 在模 \(n_{i}\) 意义下的逆元 \(m_{i}^{-1}\)

    • 计算 \(c_{i}=m_{i}m_{i}^{-1}\) (不要对 \(n_{i}\) 取模)

  • 方程组在模 \(n\) 意义下的唯一解为 \(x=\sum_{i=1}^{k}a_{i}c_{i}\pmod{n}\)

证明

我们需要证明上面算法计算所得的 \(x\) 对于任意 \(i=1,2,\dots,k\) 满足 \(x\equiv a_{i}\pmod{n_{i}}\)

\(i\ne j\) 时,有 \(m_{j}\equiv 0\pmod{n_{i}}\) ,故 \(c_{j}\equiv m_{j}\equiv 0\pmod{n_{i}}\) 又有 \(c_{i}\equiv m_{i}\times (m_{i}^{-1}\bmod{n_{i}})\equiv 1\pmod{n_{i}}\),所以我们有:

\[x\equiv \sum_{j=1}^{k}a_{j}c_{j} \pmod{n_{i}}\\ \equiv a_{i}c_{i}\pmod{n_{i}}\\ \equiv a_{i}\times m_{i}\times(m_{i}^{-1}\bmod{n_{i}})\pmod{n_{i}}\\ \equiv a_{i}\pmod{n_{i}} \]

其中 \(m_{i}^{-1}\) 表示 \(m_{i}\) 的逆元。

即对于任意 \(i=1,2,\ldots,k\),上面算法提到的 \(x\) 总是满足 \(x\equiv a_{i}\pmod{n_{i}}\),即证明了解同余方程组的算法的正确性。

因为我们没有对输入的 \(a_{i}\) 作特殊限制,所以任何一组输入都对应一个解 \(x\)

另外,若 \(x\ne y\),则总存在 \(i\) 使得 \(x\)\(y\) 在模 \(n_{i}\) 下不同余

故系数列表 \(\{a_{i}\}\) 与解 \(x\) 之间是一一映射关系,方程组总是有唯一解

P1495【模板】中国剩余定理(CRT)/ 曹冲养猪 - 洛谷

参考代码:

#include<bits/stdc++.h>
#define int long long
using namespace std;
int n,r[11],m[11],mi[11],M=1,ans;
inline void exgcd(int a,int b,int &x,int &y)//扩欧求逆元 
{
    if(b==0){x=1,y=0;return;}
    exgcd(b,a%b,x,y);
    int z=x;x=y;y=z-y*(a/b);
}
signed main()
{
    cin>>n;
    for(int i=1;i<=n;i++)
    {
        cin>>m[i];
        M*=m[i];
        cin>>r[i];
    }
    for(int i=1;i<=n;i++)
    {
        mi[i]=M/m[i];
        int x=0,y=0;
        exgcd(mi[i],m[i],x,y);//计算mi[i]在模m[i]意义下的逆元x 
        if(x<0)x+=m[i];//有可能是负数,加上模数即可 
        ans+=r[i]*mi[i]*x;//累加答案 
    }
    cout<<(ans%M)<<endl;
    return 0;
} 

Garner 算法

CRT的另一个用途是用一组比较小的质数来表示一个大的整数。

例如,若 \(a\) 满足如下线性方程组,且 \(a<\prod_{i=1}^{k}P_{i}\) (其中 \(P_{i}\) 为质数):

\[\left\{\begin{matrix} a\equiv a_{1}\pmod{p_{1}}\\ a\equiv a_{2}\pmod{p_{2}}\\ \vdots\\ a\equiv a_{k}\pmod{p_{k}} \end{matrix}\right. \]

我们可以用以下形式的式子(称作 \(a\) 的混合基数表示)表示 \(a\)

\[a=x_{1}+x_{2}p_{1}+x_{3}p_{1}p_{2}+\dots+x_{k}p_{1}\dots p_{k-1} \]

Garner 算法将用来计算系数 \(x_{1},\dots,x_{k}\)

\(r_{ij}\)\(p_{i}\) 在模 \(p_{j}\) 意义下的逆:

\[p_{i}\times r_{ij}\equiv 1\pmod{p_{j}} \]

\(a\) 带入我们得到的第一个方程:

\[a_{1}\equiv x_{1}\pmod{p_{1}} \]

代入第二个方程得到:

\[a_{2}\equiv x_{1}+x_{2}p_{1}\pmod{p_{2}} \]

方程两边减 \(x_{1}\),除 \(p_{1}\) 后得

\[a_{2}-x_{1}\equiv x_{2}p_{1}\pmod{p_{2}}\\ (a_{2}-x_{1})r_{1,2}\equiv x_{2}\pmod{p_{2}}\\ x_{2}\equiv (a_{2}-x_{1})r_{1,2}\pmod{p_{2}} \]

类似地,我们可以得到:

扩展CRT

上面的朴素 CRT 是求解

\[\left\{\begin{matrix} x\equiv a_{1}\pmod{n_{1}}\\ x\equiv a_{2}\pmod{n_{2}}\\ \vdots\\ x\equiv a_{k}\pmod{n_{k}} \end{matrix}\right. \]

一类的同余方程组,但是如果 \(n_{1},n_{2},\dots,n_{k}\) 不互质的话怎么办?

据某位学长说扩展中国剩余定理跟中国剩余定理没半毛钱关系,一个是用扩展欧几里得,一个是用构造

首先我们从简单入手,考虑只有两个式子的情况

\[x\equiv a_{1}\pmod{n_{1}}\\ x\equiv a_{2}\pmod{n_{2}} \]

将两个式子变形:

\[x=a_{1}+n_{1}k_{1}\\ x=a_{2}+n_{2}k_{2} \]

联立:

\[a_{1}+n_{1}k_{1}=a_{2}+n_{2}k_{2} \]

移项:

\[n_{1}k_{1}=a_{2}-a_{1}+n_{2}k_{2} \]

在这里需要注意就是这个方程有解的条件是 \(\gcd(n_{1},n_{2})\mid(a_{2}-a_{1})\),因为后面会用到 \(\frac{a_{2}-a_{1}}{\gcd(n_{1},n_{2})}\) 这一项,如果不能整除会出现小数。

对于上面的方程,两边同除 \(\gcd(n_{1},n_{2})\)

\[\frac{n_{1}k_{1}}{\gcd(n_{1},n_{2})}= \frac{a_{2}-a_{1}}{\gcd(n_{1},n_{2})}+ \frac{n_{2}k_{2}}{\gcd(n_{1},n_{2})}\\ \frac{n_{1}}{\gcd(n_{1},n_{2})}k_{1}= \frac{a_{2}-a_{1}}{\gcd(n_{1},n_{2})}+ \frac{n_{2}}{\gcd(n_{1},n_{2})}k_{2}\\ \]

转换一下:

\[\frac{n_{1}}{\gcd(n_{1},n_{2})}k_{1}\equiv \frac{a_{2}-a_{1}}{\gcd(n_{1},n_{2})}\pmod{ \frac{n_{2}}{\gcd(n_{1},n_{2})}} \]

此时我们把 \(k_{2}\) 搞没了

同余式两边同除 \(\frac{n_{1}}{\gcd(n_{1},n_{2})}\)

\[k_{1}\equiv inv(\frac{n_{1}}{\gcd(n_{1},n_{2})},\frac{n_{2}}{\gcd(n_{1},n_{2})})\times \frac{a_{2}-a_{1}}{\gcd(n_{1},n_{2})}\pmod{\frac{n_{2}}{\gcd(n_{1},n_{2})}} \]

\(inv(a,b)\) 表示 \(a\) 在模 \(b\) 意义下的逆元。

\[k_{1}\equiv inv(\frac{n_{1}}{\gcd(n_{1},n_{2})},\frac{n_{2}}{\gcd(n_{1},n_{2})})\times \frac{a_{2}-a_{1}}{\gcd(n_{1},n_{2})}+\frac{n_{2}}{\gcd(n_{1},n_{2})}\times y \]

我们把开头的式子拿过来,把 \(k_{1}\) 代回去。

\[x=inv(\frac{n_{1}}{\gcd(n_{1},n_{2})},\frac{n_{2}}{\gcd(n_{1},n_{2})})\times \frac{a_{2}-a_{1}}{\gcd(n_{1},n_{2})}\times n_{1}+y\frac{n_{1}n_{2}}{\gcd(n_{1},n_{2})} +a_{1} \]

\[x\equiv inv(\frac{n_{1}}{\gcd(n_{1},n_{2})},\frac{n_{2}}{\gcd(n_{1},n_{2})})\times \frac{a_{2}-a_{1}}{\gcd(n_{1},n_{2})}\times n_{1}+a_{1}\pmod{\frac{n_{1}n_{2}}{\gcd(n_{1},n_{2})}} \]

此时,整个式子中的元素我们都知道了

这个式子可以看作是

\[x\equiv a\pmod{n} \]

其中:

\[a=(inv(\frac{n_{1}}{\gcd(n_{1},n_{2})},\frac{n_{2}}{\gcd(n_{1},n_{2})}) \times \frac{a_{2}-a_{1}}{\gcd(n_{1},n_{2})})\%\frac{n_{2}}{\gcd(n_{1},n_{2})} \times n_{1}+a_{1} \]

\[n=\frac{n_{1}n_{2}}{\gcd(n_{1},n_{2})} \]

如果是多个式子可以逐个合并。

上面的式子如果你觉得又臭又长可以把 \(inv\) 及括号里的替换为用扩欧对两模数求逆元的 \(x\),他俩是一个数;\(n\) 也就是 \(\text{lcm}(n_{1},n_{2})\)

P4777【模板】扩展中国剩余定理(EXCRT) - 洛谷

#include<bits/stdc++.h>
#define int __int128
#define N 100010
using namespace std;
int n,x,y,d;
inline int read(){int x=0,f=1;char ch=getchar();while(!isdigit(ch)){f=ch!='-';ch=getchar();}while(isdigit(ch)){x=(x<<1)+(x<<3)+(ch^48);ch=getchar();}return f?x:-x;}
inline void print(int x){if(x>=10)print(x/10);putchar(x%10+48);}
inline void exgcd(int a,int b,int &x,int &y)//扩欧求逆元 
{
    if(b==0)return d=a,x=1,y=0,void();//d是当前的ab的最大公约数 
    exgcd(b,a%b,x,y);
    int z=x;x=y;y=z-y*(a/b);
}
inline int lcm(int a,int b){return a*b/__gcd(a,b);}//最小公倍数函数 
int a,b,A,B;
inline void merge()//合并两个同余方程组 
{
    exgcd(a,A,x,y);//先求两个值的一组特解和最大公约数 
    int c=B-b;//计算两个式子的余数的差值 
    if(c%d)puts("-1"),exit(0);//如果是gcd的倍数说明分母后面是0无解 
    x=x*c/d%(A/d);//计算右半部分式子的值 
    if(x<0)x+=A/d;//如果小于0就加上模数 
    int mod=lcm(a,A);//计算新模数,A,a的最小公倍数 
    b=(a*x+b)%mod;//计算余数也就是新b的值, 
    if(b<0)b+=mod;//如果b小于0就直接加上模数 
    a=mod;//新的模数 
}
signed main()
{
    n=read();a=read();b=read(); 
    for(int i=2;i<=n;i++)
    {
        int _A=read(),_B=read();//当前式子的模数,余数 
        A=_A,B=_B;
        merge();//合并 
    }
    print(b%a);
    return 0;
}

部分参考自 OI Wiki

posted @ 2023-05-26 17:30  北烛青澜  阅读(49)  评论(0编辑  收藏  举报