P4777 [模板]exCRT 题解

Posted 2021-02-06 lim-817

终于过了P4777模板题...来写篇题解祝贺一下（
考虑一个式子一个式子的合并。比如说只有两个式子：
[xequiv b_1(mod a_1)]
[xequiv b_2(mod a_2) ]
那么显然有(x=a_1 imes x_1 + b_1)，(x=a_2 imes x_2 + b_2。)
所以有[a_1 imes x_1 + b_1 = a_2 imes x_2 + b_2]
移项得[a_1 imes x_1 - a_2 imes x_2 = b_2-b_1]
然后发现它很像扩展欧几里得解线性方程组的样子，于是就可以用扩展欧几里得算法解出(x_1,x_2)。
对于解一般方程的过程这里不过多叙述，不会的读者自行学习(exGCD)。
接着我们便可以成功解出一组特解(x_0).（代入原式）
不难发现(x_0)每经过(lcm(a_1,a_2))就会出现一组循环，即一组相同的解。
所以显然可以将两个式子合并成一个，为：
[xequiv x_0(mod lcm(a1,a2))]
注意：(x_0)要变成非负的
然后再把这个式子与下个式子合并，最后就可以得出最小解了。
但是还没有完。由于过程中可能出现的乘法溢出，我们结合代码来说明。

void exCRT() {
    for(int i = 2; i <= n; i++) {
        LL A = a[i - 1] , B = -a[i] , C = b[i] - b[i - 1];
        LL x , y; exGCD(A , B , x , y);
        x *= C / __gcd(A , B);
        LL X0 = a[i - 1] * x + b[i - 1]; b[i] = X0;
        a[i] = a[i] * a[i - 1] / __gcd(a[i] , a[i - 1]);
        b[i] = ((b[i] + a[i]) % a[i] + a[i]) % a[i];
    }
    cout << b[n] << endl;
}

这是我的第一版代码。虽然说都开了long long，但是还是避免不了过程中的乘法溢出。
从上往下看，发现好像每一步都不能成功取模，那么我们就来考虑从下往上看。
发现最后一步(b[i])通过取模(a[i])来得到非负整数解，也就是说过程中只要只跟(a[i])有关的变量一直都可以模(a[i])，对结果没有影响。
然后观察(a[i])的计算方式，发现这里两个long long相乘肯定会爆掉，但是(lcm)不会，要不然不会保证有非负整数解，并且(leq 10^{18}).
所以首先把(a[i])的计算提到前面，并且将其改成

LL ans = a[i] / __gcd(a[i] , a[i - 1]);
        ans *= a[i - 1]; a[i] = ans;

这样即可避免计算(lcm)中出现的溢出问题。
接着发现(X0)处计算时朴素的乘还是可能会爆long long，所以考虑采用快速乘（两个long long相乘取模一个数）
由于笔者之前写过快速乘并在某场比赛中成功丢掉50分，我很不想写递归版的快速乘（带个(log)），于是这里推荐一种黑科技：

template < typename T > inline T mul(T x , T y , T MOD) {x=x%MOD,y=y%MOD;return ((x*y-(T)(((ld)x*y+0.5)/MOD)*MOD)%MOD+MOD)%MOD;}

是(Theta(1))的，非常快。
然后要注意，当我们写一个式子的时候，比如说(xequiv 6(mod 10))，它和(xequiv 3(mod 5))是等价的。
所以一开始一定要先将(a,b)化简，即为同时除以(gcd(a,b).)
然后就可以了。
最终代码如下：

#include <bits/stdc++.h>
#define dbg(x) cerr << #x " = " << x << "
"
#define INF 0x3f3f3f3f

typedef long long LL;
typedef long double ld;
typedef unsigned long long ULL;

using namespace std;

template < typename T > inline void inp(T& t) {
    char c = getchar(); T x = 1; t = 0; while(!isdigit(c)) {if(c == ‘-‘) x = -1; c = getchar();}
    while(isdigit(c)) t = t * 10 + c - ‘0‘ , c = getchar();t *= x;
}
template < typename T , typename... Args > inline void inp(T& t , Args&... args) {inp(t); inp(args...);}
template < typename T > inline void outp(T t) {
    if(t < 0) putchar(‘-‘) , t = -t; T y = 10 , len = 1;
    while(y <= t) y *= 10 , len++; while(len--) y /= 10 , putchar(t / y + 48) , t %= y;
}
template < typename T > inline T mul(T x , T y , T MOD) {x=x%MOD,y=y%MOD;return ((x*y-(T)(((ld)x*y+0.5)/MOD)*MOD)%MOD+MOD)%MOD;}

const int MAXN = 100000 + 5;
LL a[MAXN] , b[MAXN];
int n;

void exGCD(LL a , LL b , LL& x , LL& y) {
    if(!b) {
        x = 1;
        y = 0; return ;
    }
    exGCD(b , a%b , x , y);
    LL t = x; x = y; y = t - (a / b) * y;
    return ;
}
void exCRT() {
    for(int i = 2; i <= n; i++) {
        LL A = a[i - 1] , B = -a[i] , C = b[i] - b[i - 1];
        LL g = __gcd(A , B);
        A /= g , B /= g; C /= g;
        LL x , y; exGCD(A , B , x , y);
        x = mul(x , C , a[i]);
        LL X0;
        LL ans = a[i] / __gcd(a[i] , a[i - 1]);
        ans *= a[i - 1]; a[i] = ans;
        X0 = (mul(a[i - 1] , x , a[i]) + b[i - 1] % a[i]) % a[i];
        X0 = ((X0 + a[i]) % a[i] + a[i]) % a[i];
        b[i] = X0;
    }
    cout << b[n] << endl;
}
int main() {
    inp(n);
    for(int i = 1; i <= n; i++) inp(a[i] , b[i]); 
    exCRT(); 
    return 0;
}

终于写完了...