「算法笔记」线性同余方程

线性同余方程是由一系列同余式组成的，本文主要讲解一元一次同余方程的求解。

基本算法

丢番图方程

丢番图方程是形如 $ax+by=c$ 的方程，$ax\equiv b\pmod m$ 也是其常见形式之一。

裴蜀定理：丢番图方程 $ax+by=c$ 有解当且仅当 $\gcd(a,b)\mid c$。

扩展欧几里德算法

扩展欧几里德算法简称 $\text{exGCD}$，用于求解丢番图方程的整数解。

推导过程：

$\begin{cases} ax_1+by_1=\gcd(a,b) \\ bx_2+(a\bmod b)y_2=\gcd(b,a\bmod b) \end{cases}$

我们可以通过 $\gcd(a,b)=\gcd(b,a\bmod b)$ 得到：

$\begin{aligned} ax_1+by_1&=bx_2+(a\bmod b)y_2 \\ &=bx_2+(a-\left\lfloor\frac{a}{b}\right\rfloor\times b)y_2 \\ &=ay_2+b(x_2-\left\lfloor\frac{a}{b}\right\rfloor\times y_2) \end{aligned}$

所以：

$\begin{cases} x_1=y_2 \\ y_1=x_2-\left\lfloor\frac{a}{b}\right\rfloor\times y_2 \end{cases}$

很显然这是一个递归式，末状态为 $b=0,a=\gcd(a,b)$ 时，$x=1,y=0$。

扩展欧几里德的过程可以理解为从末状态向上不断回溯的过程，直到得到原方程的一组解。

方程通解

设 $(x_0,y_0)$ 为方程 $ax+by=c$ 的一组特解，那么方程的通解为（其中 $K$ 为一个系数）：

$(x_0+K\frac{b}{\gcd(a,b)},y_0-K\frac{a}{\gcd(a,b)})$

将这个通解代入方程易证。

求方程组的解

中国剩余定理（CRT）

求解

设 $M=\prod_{i=1}^n m_i$，并设 $M_i=\frac{M}{m_i}$，$t_i={M_i}^{-1}$ 为 $M_i$ 模 $m_i$ 的数论倒数（逆元），那么可以得到方程在模 $M$ 意义下的唯一解：

$x=\left(\sum_{i=1}^n r_it_iM_i\right)\bmod M$

时间复杂度：$O(n\log M)$

证明

$x=\sum_{i=1}^n r_it_iM_i\equiv r_j\pmod m_j$

当 $i\neq j$ 时，$\gcd(M_i,m_j)=1,\gcd(t_i,m_j)=1$，显然不满足条件。

当 $i=j$ 时，显然满足条件。

故中国剩余定理的本质思想就是让每一项只对自己有贡献。

代码

由于精度问题，这里的 $n$ 真的不能开得再大了……

#include <cstdio>
typedef long long LL;

const int N=105;
int n;
LL a[N],m[N];

LL exgcd(LL a,LL b,LL &x,LL &y) {
    if(!b) {x=1,y=0;return a;}
    LL d=exgcd(b,a%b,y,x);
    y-=a/b*x; return d;
}
LL China(int n) {
    LL M=1,ans=0;
    for(int i=1;i<=n;++i) M*=m[i];
    for(int i=1;i<=n;++i) {
        LL k=M/m[i],x,y;
        exgcd(k,m[i],x,y);
        ans=(ans+x*k%M*a[i]%M)%M;
    }
    return (ans+M)%M;
}
int main() {
    scanf("%d",&n);
    for(int i=1;i<=n;++i) scanf("%lld%lld",&m[i],&a[i]);
    printf("%lld\n",China(n));
    return 0;
}

扩展中国剩余定理（ExCRT）

求解

由于中国剩余定理有特殊限制：$m_i$ 必须两两互质。我们必须找到一个运用更广泛的算法来求出任意线性同余方程组的解。

我们这次直接考虑两个相邻的方程：

$\begin{cases} x\equiv r_1\pmod {m_1} \\ x\equiv r_2\pmod {m_2} \end{cases}$

我们把他们拆开得到：$x=k_1m_1+r_1=k_2m_2+r_2$（其中 $k_1$ 和 $k_2$ 均为系数）。

这个式子如果不看 $x$ 则变成：

$k_1m_1-k_2m_2=a_2-a_1$

我们把 $k_1$ 和 $k_2$ 看作未知数然后解该丢番图方程。

设 $k_1$ 的通解为 $k+K\frac{m_2}{\gcd(m_1,m_2)}$，那么带回到 $x=k_1m_1+r_1$ 中可以得到：

$x=km_1+r_1+K\frac{m_1m_2}{\gcd(m_1,m_2)}$

注意到上述式子中，$\frac{m_1m_2}{\gcd(m_1,m_2)}$ 就是 $\operatorname{lcm}(m_1,m_2)$，那么我们根据模运算的性质可以得到一组新的同余方程：

$x\equiv km_1+r_1\pmod {\operatorname{lcm}(m_1,m_2)}$

我们用这样的方法把相邻的同余方程两两合并，每次消去一个方程，最终只会剩下一个同余方程，这个同余方程的余数 $r’$ 一定满足所有的同余方程。

时间复杂度：$O(n\log m_i)$

代码

#include <cstdio>
typedef long long LL;

const int N=1e6+5;
int n;
LL r[N],m[N];

LL exgcd(LL a,LL b,LL &x,LL &y) {
    if(!b) {x=1,y=0;return a;}
    LL d=exgcd(b,a%b,y,x);
    y-=a/b*x; return d;
}
LL China(int n) {
    LL M=m[1],R=r[1];
    for(int i=2;i<=n;++i) {
        LL a=M,b=m[i],c=r[i]-R,x,y,d=exgcd(a,b,x,y);
        if(c%d) return -1;
        a/=d,b/=d,c/=d,x=(x*c%b+b)%b;
        R+=x*M,M*=m[i]/d,R=(R+M)%M;
    }
    return R;
}
int main(n) {
    scanf("%d",&n);
    for(int i=1;i<=n;++i) scanf("%lld%lld",&m[i],&r[i]);
    printf("%lld\n",exCRT());
    return 0;
}