浅谈扩展中国剩余定理

洛谷p4777

前言

在洛谷上，这是一个连模板都是省选/NOI-的题目，我也钻研了好几天，也许还没有达到通透的境界，但我还是想写一些我在学习过程中的一些感想。

用途

这个定理，其实不能说是一个定理，我觉得更像一个算法。
给出如下同余方程组，让你求方程的最小正数解：
$x\equiv b_1$(mod $a_1)$
$x\equiv b_2$(mod $a_2)$
……
$x\equiv b_i$(mod $a_i)$

算法

我们先从简单的着手，一步一步来。先看前两个方程。
根据第一个方程，我们可以设$x=a_1\ast x_1+b_1$。同理，根据第二个方程，我们可以设$x=a_2\ast (-x_2)+b_2$。（$x_2$前面加个负号是为了方便下面的移项）
因为x是相等的，所以$a_1\ast x_1+b_1=a_2\ast (-x_2)+b_2$
移项得：$a_1\ast x_1+a_2\ast x_2=b_2-b_1$
$a_1\ast x_1^{\prime }+a_2\ast x_2^{\prime }=gcd(a_1,a_2)$的解是可以通过扩展欧几里得求出来的
如果$b_2-b_1$不是$gcd(a_1,a_2)$的倍数，那么无解。
然后有$a_1\ast x_1^{\prime }\ast (b_2-b_1)/gcd(a_1,a_2)+a_2\ast x_2^{\prime }\ast (b_2-b_1)/gcd(a_1,a_2)=b_2-b_1$
$x_1=x_1^{\prime }\ast (b_2-b_1)/gcd(a_1,a_2)，x_2=x_2^{\prime }\ast (b_2-b_1)/gcd(a_1,a_2)$
当然，用扩展欧几里得求出的只是一组解，而我们希望得到的是最小的那一组解，所以要调整一下$x_1$。
把$x_1$变成$(x_1\%a_2+a_2)\%a_2$，这样$x_1$就是最小正数解了。
那么$x=(a_1\ast x_1+b_1)\%(a_1\ast a_2/gcd(a_1,a_2))$
记最终的答案为ans，$a_1\ast a_2/gcd(a_1,a_2)$为m，那么$ans\equiv x$(mod $m$)
这样我们就构造出了一个新的同余方程，来代替这两个方程。

推向一般化，我们现在有这两个同余方程：
$ans\equiv x$(mod $m$)
$ans\equiv b_i$(mod $a_i$)
按照上面的步骤再走一遍。
根据第一个方程，设$ans=m\ast x_1+x$。同理，根据第二个方程，设$ans=a_i\ast (-x_2)+b_i$。
那么$m\ast x_1+x=a_i\ast (-x_2)+b_i$
移项得：$m\ast x_1+a_i\ast x_2=b_i-x$
令$z=b_i-x$，要是z很大或者是个大负数，那将会面临溢出的风险，所以要取模。但对什么取模呢？因为在后面$x_2$的值根本没用到，所以z可以对$a_i$取模，这样只对$x_2$的值有影响，对$x_1$并没有影响。
根据扩展欧几里得求出$m\ast x_1^{\prime }+a_i\ast x_2^{\prime }=gcd(m,a_i)$的一组解
如果$z$不是$gcd(m,a_i)$的倍数，那么无解。
然后有$m\ast x_1^{\prime }\ast z/gcd(m,a_i)+a_i\ast x_2^{\prime }\ast z/gcd(m,a_i)=b_i-x$
那么$x_1=x_1^{\prime }\ast z/gcd(m,a_i)$，$x_2=x_2^{\prime }\ast z/gcd(m,a_i)$（其实$x_2$无所谓，反正后面用不到）
再调整$x_1$，使之变为最小正数解：$x_1=(x_1\%a_i+a_i)\%a_i$
然后构造新的同余方程：$x=(m\ast x_1+x)\%(m/gcd(m,a_i)\ast a_i)$，$m=m/gcd(m,a_i)\ast a_i$
则新的同余方程组还是$ans\equiv x$(mod $m$)（ans只是存在于我们脑袋中的一个量，最终的答案就是x）
由于数据范围非常大，上面的带取模的乘法最好用快速乘，不然可能会溢出。

具体实现看我代码：

#include<iostream>
#include<cmath>
#include<cstdio> 
using namespace std;
int n;
long long a[100005],b[100005];
long long exgcd(long long a,long b,long long &x,long long &y)
{
    if (b==0)
    {
        x=1;
        y=0;
        return a;
    }
    long long gcd=exgcd(b,a%b,x,y);
    int t=x;
    x=y;
    y=t-a/b*y;
    return gcd;
}
long long cheng(long long x,long long y,long long z)//快速乘 
{
    long long t=0;
    x=(x%z+z)%z;//快速乘一定要加这两句话，因为x可能是负数
    y=(y%z+z)%z;
    while (x>0)
    {
        if (x%2==1) t=(t+y)%z;
        y=y*2%z;
        x/=2;
    }
    return t;
}
long long excrt()
{
    long long x=b[1],x1,x2,m=a[1];
    for (int i=2;i<=n;i++)
    {
    	long long z=((b[i]-x)%a[i]+a[i])%a[i];//为避免溢出，要对a[i]取模 
        long long gcd=exgcd(m,a[i],x1,x2);//求m*x1+a[i]*x2=gcd(m,a[i])的一组解，gcd=gcd(m,a[i]) 
		if (z%gcd!=0) return 0;//说明无解 
        x1=cheng(x1,z/gcd,a[i]);//x1=x1*z/gcd，同时调整x1为最小正数解 
        x=cheng(m,x1,m/gcd*a[i])+x;//构造新的同余方程 
        m=m/gcd*a[i];
        x=x%m;//上上行加了x之后还要取一次模
    }
    return x;//x就是最终是答案 
}
int main()
{
    cin>>n;
    for (int i=1;i<=n;i++)
    {
        scanf("%lld %lld",&a[i],&b[i]);
    }
    cout<<excrt();
}

总结

这虽然只是一个模板，但难度却也不小，我学了好久，总是因为一些取模而困惑，希望这一篇文章把细节都写到了。像这种数据规模很大的题目，每一步都要想想能不能取模，对什么取模。