「算法笔记」CRT 与 exCRT

一、扩展欧几里得

求解方程 \(ax+by=\gcd(a,b)\)。

int exgcd(int a,int b,int &x,int &y){
    if(!b) return x=1,y=0,a;
    int d=exgcd(b,a%b,x,y);
    int z=x; x=y,y=z-y*(a/b);
    return d;
} 

对于更为一般的方程 \(ax+by=c\)，设 \(d=\gcd(a,b)\)。我们可以求出 \(ax+by=d\) 的一组特解 \(x_0,y_0\)。这所以 \(\frac{c}{d}\) 也必须为整数，否则无解（所以该方程有解当且仅当 \(d\mid c\)）。我们令 \(x_0,y_0\) 同时乘上 \(\frac{c}{d}\)，得 \(a\frac{cx_0}{d}+b\frac{cy_0}{d}=c\)。那么 \(x=\frac{c}{d}\cdot x_0,y=\frac{c}{d}\cdot y_0\) 即为一组解。

事实上，方程 \(ax+by=c\) 的通解可以表示为：

\[x=\frac{c}{d}x_0+k\frac{b}{d},\,y=\frac{c}{d}y_0-k\frac{a}{d}\,(k\in \mathbb{Z}) \]

那么若求得一个 \(x\)，则该方程 \(x\) 的最小整数解就是 (x%(b/d)+b/d)%(b/d)。

二、中国剩余定理

1. 求解

考虑求解方程组 \(x\equiv r_i\pmod {m_i}\)。即：

\[\begin{cases} x\equiv r_1\pmod {m_1}\\ x\equiv r_2\pmod {m_2}\\ \ \ \ \ \vdots\\ x\equiv r_k\pmod {m_k} \end{cases} \]

其中 \(m_i\) 两两互质。中国剩余定理（CRT）的表述如下：

令 \(M=\prod_{i=1}^k m_i\)，则该方程组在 \([0,M)\) 范围内有唯一解，且显然是最小整数解。设 \(M_i=\frac{M}{m_i}\)，\(t_i=M_i^{-1}\) 表示 \(M_i\) 在模 \(m_i\) 意义下的逆元，即 \(M_it_i\equiv 1\pmod {m_i}\)。则最小解为：

\[x=\left(\sum_{i=1}^k r_it_iM_i\right)\mod M \]

2. 证明

\[x=\sum_{i=1}^k r_it_iM_i\equiv r_j \pmod {m_j} \]

当 \(i\neq j\) 时，\(M_i=\frac{m_1m_2\cdots m_k}{m_i}\) 必为 \(m_j\) 的倍数，则有 \(r_it_iM_i\equiv 0\pmod {m_j}\)。

否则，又因为 \(M_jt_j\equiv 1\pmod {m_i}\)，所以 \(r_jt_jM_j\equiv r_j\pmod {m_j}\)。

综上，对于任意 \(i\)，总有 \(x=\sum_{i=1}^k r_it_iM_i\equiv r_i\pmod {m_i}\)。

三、扩展中国剩余定理

考虑 \(m_i\) 不互质的情况。

设当前求解到第 \(i\) 个方程，前 \(i-1\) 个方程的最小解为 \(X\)，\(\text{lcm}(m_1,m_2,\cdots,m_{i-1})=M\)。则前 \(i-1\) 个方程的通解为 \(X+tM\)（\(t\) 为任意非负整数）。

考虑第 \(i\) 个方程，现在要找一个 \(t\)，使得 \(X+tM\equiv r_i\pmod {m_i}\)。

移项得：\(tM\equiv r_i-X\pmod {m_i}\)。进一步转化为：\(tM+pm_i=r_i-X\)。显然这是一个不定方程，可以用 exgcd 求解。

解出最小的 \(t\) 后，更新 \(X=X+tM\)，\(M=\text{lcm}(M,m_i)\)。

最后为保证 \(X\) 最小，还需对 \(M\) 取模。

//Luogu P4777
#include<bits/stdc++.h>
#define int long long
#define LLL __int128    //过程中某个地方会炸。或改用快速乘。
using namespace std;
const int N=1e5+5;
int n,m[N],r[N],M;
int exgcd(int a,int b,int &x,int &y){
    if(!b) return x=1,y=0,a;
    int d=exgcd(b,a%b,x,y);
    int z=x; x=y,y=z-y*(a/b);
    return d;
} 
int exCRT(int n){
    int M=m[1],X=r[1];
    for(int i=2;i<=n;i++){
        int a=M,b=m[i],c=r[i]-X,x,y,d=exgcd(a,b,x,y),t;
        if(c%d) return -1;
        a/=d,b/=d,c/=d,t=((LLL)x*c%b+b)%b;
        X+=t*M,M*=m[i]/d,X=(X+M)%M;
    } 
    return X;
}
signed main(){
    scanf("%lld",&n);
    for(int i=1;i<=n;i++)
        scanf("%lld%lld",&m[i],&r[i]); 
    printf("%lld\n",exCRT(n));
    return 0;
}