保洁阿姨都能看懂的中国剩余定理和扩展中国剩余定理

背景：

这个定理又叫孙子定理，问题的提出最早在南北朝的《孙子兵法算经》里，“有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？”，
即，一个整数除以三余二，除以五余三，除以七余二，求这个整数。
 
 

问题：

求解同余方程组，形如：（盗图，版权意识薄弱，读书人的事情那能叫偷吗？）
解决这个问题前，我们要清楚的一点。令$M=lcm(m_1,m_2,\dots ,m_n)$，如果方程有解的话，答案一定是$x_0+kM$的形式，也就是解在模M意义下唯一。
首先$x_0$是方程的一组特解，其次不管代入哪个方程$kM$总是会被消去。
 
 
 

中国剩余定理：传送门

首先提出问题的简单版本，即方程的模数两两互质。
这就是裸裸 的孙子定理了。
孙子的做法是：我们只要想办法去构造一组解就好了。
那么怎么构造？
 
这个构造方法和**牛顿插值法拉格朗日插值法 **类似（感谢高中数学老师）。
思考序列$1,2,3,\pi$，构造一个通项公式，使得$f(1)=1,f(2)=2,f(3)=3,f(4)=\pi$。
利用拉格朗日插值法构造出函数：
$f(x)=1\times \frac{(x-2)(x-3)(x-4)}{(1-2)(1-3)(1-4)}+$$2\times \frac{(x-1)(x-3)(x-4)}{(2-1)(2-3)(2-4)}+$$3\times \frac{(x-1)(x-2)(x-4)}{(3-1)(3-2)(3-4)}+$$\pi \times \frac{(x-1)(x-2)(x-3)}{(4-1)(4-2)(4-3)}$
这个神仙原理比较显而易见，每一项由两部分构成，$\times$左边是每一项的数，$\times$右边是决定项是否存在的控制因子，当$x$为$1$时，除第$1$项外的其他三项控制因子为$0$，而第一项的为$1$，这就保证了，对应的$x$可以得到我们想要的对应项。
 
然后回到我们的问题，应该如何构造一组解，满足所有的方程。
我们令$M=m_1{m}_2{m}_3\dots m_n$。
再令$M_i=M/m_i$，就是我们的控制因子了。
然后我们要开始构造特解$x_0$了。
按照前面的思路暴躁地 构造就好了。
我们先暴力地 把每一项都填到式子里。
$x_0=a_1+a_2+\dots +a_n$
然后我们把控制因子$M_i$填进去。
$x_0=a_1{M}_1+a_2{M}_2+\dots +a_n{M}_n$
此时把它带入到第$i$个方程中，除$a_i{M}_i$以外的其他项都会被消掉，只剩这一项。
然后我们要做的就是把$M_i$也消掉，只剩下$a_i$，所以我们乘上它的逆元，就变成了最后这个式子了。
$x_0=a_1{M}_1{M}_1^{-1}+a_2{M}_2{M}_2^{-1}+\dots +a_n{M}_n{M}_n^{-1}$
这里$M_i^{-1}$是模$m_i$意义下的，前面都懂了的话，应该是知道的。
 

最后的式子要再列一遍：

$x_0\equiv a_1{M}_1{M}_1^{-1}+a_2{M}_2{M}_2^{-1}+\dots +a_n{M}_n{M}_n^{-1}(modM)$
还有，记得前提是要两两互素，不然求逆元那步会出现没有逆元的情况。

模板：(修正，改用了快速乘防止溢出)

#include <bits/stdc++.h>
using namespace std;
typedef __int128 ll;

ll fmul(ll a,ll b,ll mod)
{
    
    ll sum=0,base=(a%mod+mod)%mod;
    while(b)
    {
    
        if(b%2)sum=(sum+base)%mod;
        base=(base+base)%mod;
        b/=2;
    }
    return sum;
}

ll read()
{
    
   int X=0,w=0; char ch=0;
   while(!isdigit(ch)) {
    w|=ch=='-';ch=getchar();}
   while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();
   return w?-X:X;
}
void print(ll x)
{
    
   if(x<0){
    putchar('-');x=-x;}
   if(x>9) print(x/10);
   putchar(x%10+'0');
}

ll ex_gcd(ll a,ll b,ll& x,ll& y)
{
    
    if(b==0)
    {
    
        x=1;y=0;
        return a;
    }
    ll ans=ex_gcd(b,a%b,x,y);
    ll tmp=x;
    x=y;
    y=tmp-a/b*y;
    return ans;
}

ll inv(ll a,ll mod)//存在逆元条件：gcd(a,mod)=1
{
    
    ll x,y;
    ll g=ex_gcd(a,mod,x,y);
    if(g!=1)return -1;
    return (x%mod+mod)%mod;
}

ll a[100005],m[100005];
ll crt(ll *a,ll *m,ll n)//长度为0到n-1
{
    
    ll M=1;
    for(int i=0;i<n;i++)M=M*m[i];
    ll ans=0;
    for(int i=0;i<n;i++)
    {
    
        ll MM=M/m[i];
        ans=(ans+fmul(fmul(a[i],MM,M),inv(MM,m[i]),M))%M;
    }
    return ans;
}

int main()
{
    
    ll n;
    n=read();
    for(int i=0;i<n;i++)
    {
    
        m[i]=read();a[i]=read();
    }
    print(crt(a,m,n));
    return 0;
}

扩展中国剩余定理：传送门

现在要解决问题的升级版，也就是没有模数两两互素的条件下。
既然是扩展的，那么举一反三
这是一种与前面完全不同的做法了。
考虑，当前面$i-1$个方程的通解是$c+kM$（$M$是前面模数的最小公倍数）时,我们加入第$i$个方程$x\equiv a_i(mod{m}_i)$。$M$会变成加入新的模数后的最小公倍数，我们要求在模新$M$意义下的解。
我们把前面的通解带入新的方程，得到$c+kM\equiv a_i(mod{m}_i)$，变形为$kM\equiv a_i-c(mod{m}_i)$。
注意此时不能对两边同乘$M$的逆元（因为这个调了一下午的bug），因为此时不能保证$M$有逆元。我们要求解新的$k$，将方程转化成$Mk+m_{i}y=a_i-c$。
此时如果方程无解，则方程组无解。
我们解出$k$的值，然后加入第$i$个方程后新的$c$就变成了$c+kM$，而M则是变成加新数的最小公倍数。
我们可以设初解为模$1$意义下为$0$，然后将方程组逐一加入。

模板：

#include <bits/stdc++.h>
using namespace std;
typedef __int128 ll;

inline ll read()
{
    
   ll X=0,w=0; char ch=0;
   while(!isdigit(ch)) {
    w|=ch=='-';ch=getchar();}
   while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();
   return w?-X:X;
}
inline void print(ll x)
{
    
   if(x<0){
    putchar('-');x=-x;}
   if(x>9) print(x/10);
   putchar(x%10+'0');
}

ll fmul(ll a,ll b,ll mod)
{
    
    ll sum=0,base=(a%mod+mod)%mod;
    while(b)
    {
    
        if(b%2)sum=(sum+base)%mod;
        base=(base+base)%mod;
        b/=2;
    }
    return sum;
}

ll gcd(ll a,ll b)
{
    
    return b?gcd(b,a%b):a;
}

ll ex_gcd(ll a,ll b,ll& x,ll& y)
{
    
    if(b==0)
    {
    
        x=1;y=0;
        return a;
    }
    ll g=ex_gcd(b,a%b,x,y);
    ll tmp=x;
    x=y;
    y=tmp-a/b*y;
    return g;
}

ll a[100005],m[100005];
ll ex_crt(ll *a,ll *m,ll n)//长度为0到n-1
{
    
    ll M=1,c=0;
    for(int i=0;i<n;i++)
    {
    
        ll t=(a[i]%m[i]-c%m[i]+m[i])%m[i];
        ll x,y;
        ll g=ex_gcd(M,m[i],x,y);
        if(t%g!=0)return -1;
        ll tM=M;
        M*=m[i]/gcd(M,m[i]);
        x=fmul(t/g,(x%M+M)%M,M);
        c=(c%M+fmul(x,tM,M)+M)%M;
    }
    return (c%M+M)%M;
}

int main()
{
    
    ll n;
    n=read();
    for(int i=0;i<n;i++)
    {
    
        m[i]=read();a[i]=read();
    }
    print(ex_crt(a,m,n));
    putchar('\n');
    return 0;
}