孙子定理

孙子定理

孙子定理，又称之为中国剩余定理（Chinese Remainder Theorem, CRT）可以求解如下形式的一元线性同余方程组（其中\(n_1,n_2,\dots,n_k\)两两互质）。

中国剩余定理说明：假设整数\(m_1,m_2,\dots,m_n\)两两互质，则对任意的整数：\(a_1,a_2,\dots,a_n\)，方程组有解，并且我们可以用如下方式构造得到：

设\(M=m_1\times m_2\times\dots\times m_n=\prod_{i=1}^{n}m_i\)是整数\(m_1,m_2,\dots,m_n\)的乘积，并设\(M_i=\frac{M}{m_i},\forall i\in\{1,2,3,\dots,n\}\)是除了\(m_i\)以外的\(n-1\)个整数的乘积。设\(t_i\)是\(M_i\)模\(m_i\)下的逆元，则方程组的通解形式为：\(x=a_1t_1M_1+a_2t_2M_2+\dots+a_nt_nM_n+kM=kM+\sum_{i=1}^na_it_iM_i,k\in Z\)。
所以在模\(M\)的意义下，方程组的只有一个解：\(x=(\sum_{i=1}^na_it_iM_i)(mod\:M)\)。

算法流程：

1. 计算所有模数的积\(n\)；(有的上面说计算它们的最小公倍数，其实都一样它们两两互质)
2. 对于第\(i\)个方程：
   1. 计算\(m_i=\frac{n}{n_i}\);
   2. 计算\(m_i\)在模\(n_i\)意义下的逆元\(m_i^{-1}\)；
   3. 计算\(C_i=m_i\times m_i^{-1}\)（不要对\(n_i\)取模）；
3. 方程组的唯一解为：\(a=\sum_{i=1}^ka_i\times c_i(mod\:n)\)。

下面我们来通过一道题，加深对孙子定理的运用

猜数字

现有两组数字，每组\(k\)个，第一组中的数字分别用\(a_1,a_2,\dots,a_k\)表示，第二组中的数字分别用\(b_1,b_2,\dots,b_k\)表示。其中第二组的数字是两两互素的。求最小的\(n\in N\)，满足对于\(\forall i\in[i,k]\)，有\(b_i|(n-a_i)\)。
输入格式
第一行一个整数\(k\)。
第二行\(k\)个整数，表示：\(a_1,a_2,\dots,a_k\)。
第三行\(k\)个整数，表示：\(b_1,b_2,\dots,b_k\)。
输出格式
输出一行一个整数，为所求的答案\(n\)。\(1≤k≤10，∣a_i∣≤10^9，1≤b_i≤6×10^3\prod_{i=1}^k b_i\le 10^{18}\)。

这是孙子定理的典型运用，同样也是一道模板题，在这我我们观察到数据的范围很大，所以在计算的最后一步我们不能直接使用乘法运算，而是使用快速乘。

#include <bits/stdc++.h>
#include <map>
using namespace std;
typedef long long ll;
#define endl '\n'
#define IOS ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
#define _for(i, a, b) for (int i=(a); i<=(b); i++)
const int INF = 0x7fffffff;
const int MAXN = 1e5 + 10;

ll m[MAXN], a[MAXN], n;

void exgcd(ll a, ll b, ll &x, ll &y){
	if(b == 0){
		x = 1;
		y = 0;
		return ;
	}
	exgcd(b,a%b,x,y);
	ll temp = x;
	x = y;
	y = temp - a / b * y;
	return ;
}

ll qmul(ll n, ll a, ll mod){
	ll res = 0;
	a = a % mod;
	while(n){
		if(n & 1) res = (res + a) % mod;
		n >>= 1;
		a = (a + a) % mod;
	}
	return res % mod;
}

ll CHT(){
	ll res = 0, M = 1;
	for(int i = 1; i <= n; i++){
		M *= m[i];
	}
	for(int i = 1; i <= n; i++){
		ll x, y;
		ll Mi = M / m[i];
		exgcd(Mi,m[i],x,y);
		x = (x % m[i] + m[i]) % m[i];//防止x为负数
		res = (res + qmul(qmul(Mi,x,M),a[i],M)) % M;//使用快速乘不然有的会爆long long
	}
	if(res < 0) res += M;
	return res;
}

int main(){
	cin >> n;
	for(ll i = 1; i <= n; i++){
		cin >> a[i];
	}
	for(int i= 1; i <= n; i++){
		cin >> m[i];
	}
	for(int i = 1; i <= n; i++){
		a[i] = (a[i] % m[i] + m[i]) % m[i];
	}
	cout << CHT();
	return 0;
}

扩展孙子定理——求解模数不互质情况下的线性方程组

将\(1,2\)两式合并我们可以得到：\(m_1\times x_1+a_1=m_2\times x_2+a_2\),即\(m_1\times x_1+(-m_2)\times x_2=a_2-a_1\)。运用扩展欧几里得算法，我们可以算出\(x_1,x_2\)进而我们可以知道\(x_1,x_2\)的解集\(x_1=x_1+k\times\frac{m_2}{d},x_2=x_2+k\times\frac{m_1}{d}\)(\(k\)为任意常数)。其中\(gcd(m_1,m_2)=d\)，并且我们知道\(gcd(m_1,m_2)|(a_2-a_1)\)，否则无解。这里对于\(x_1,x_2\)解集的证明可以自行证明。

接下来，将\(x_1\)的解集单独带入\((1)\)式中：\(x=m_1\times x_1+a_1+k\times\frac{m_1\times m_2}{d}\)。这里我们知道\(m_1,x_1,a_1，[m_1,m_2]\)的值是已经确认的，只有\(k\)的值是不定的，所以我们可以将式子看成\(x=x_0+k\times a\)（即\(x\equiv x_0(mod\:a)\)）。所以我们不断重复\(n-1\)次，即可将\(n\)个式子合并，所以最后的结果为\(x_0\:mod\:a\)。

同样的我们通过一道题来更加详细的了解它。

扩展中国剩余定理

给定\(n\)组非负整数\(a_i,b_i\)，求解关于\(x\)的方程组的最小非负整数解。

输入格式
输入一行包含整数\(n\)。
接下来\(n\)行，每行两个非负整数\(a_i,b_i\)。
输出格式
输出一行，为满足条件的最小非负整数\(x\)。
数据范围
\(1\leq n\leq10^5,1\leq b_i,a_i\leq10^{12}\)，保证所有\(a_1\)的最小公倍数不超过\(10^{18}\)。

这里也是直接运用扩展孙子定理，但是我们可以看到数据的范围很大，所以在编写代码时，我们不能直接使用乘法运算，而是使用快速乘。但是对于最后求最小公倍数我们必须使用乘法运算，怎么办呐(╯▔皿▔)╯，这里我们可以先除再乘。

代码如下：

#include <bits/stdc++.h>
#include <map>
using namespace std;
typedef long long ll;
#define endl '\n'
#define IOS ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
#define _for(i, a, b) for (int i=(a); i<=(b); i++)
const int INF = 0x7fffffff;
const int MAXN = 1e5 + 10;

ll n, m[MAXN], a[MAXN];

ll exgcd(ll a, ll b, ll &x, ll &y){
    if(b == 0){
        x = 1;
        y = 0;
        return a;
    }
    ll gcd = exgcd(b,a%b,x,y);
    ll temp = x;
    x = y;
    y = temp - a / b * y;
    return gcd;
}

ll qmul(ll n, ll b, ll mod){
    ll res = 0;
    while(b > 0){//这里使用b是因为防止n为负数
        if(b & 1) res = (res + n) % mod;
        n = (n + n) % mod;
        b >>= 1;
    }
    return res;
}

int main(){
    IOS
    cin >> n;
    for(int i = 1; i <= n; i++){
        cin >> m[i] >> a[i];
    }
    ll m1 = m[1], a1 = a[1];    
    ll x, y;
    int ok  = 1;
    for(int i = 2; i <= n; i++){
        ll a2 = a[i], m2 = m[i];
        ll c = ((a2 - a1) % m2 + m2) % m2;//防止c为负数
        ll gcd = exgcd(m1,m2,x,y);
        if(c % gcd){ ok = 0; break;}
        x = qmul(x,c/gcd,m2);//防止爆long long
        a1 = a1 + x * m1;
        m1 = m2 / gcd * m1;//先除后成，防止爆long long
        a1 = (a1 + m1) % m1;
    }
    if(ok == 0) cout << -1;
    else    cout << a1;
    return 0;
}