中国剩余定理

中国剩余定理

概念

剩余定理又称孙子算经。
来源于古人的问题：今有物不知其数，三三数之余二；五五数之余三；七七数之余二。问物几何？
古人解决方案：三人同行七十稀，五树梅花甘一枝，七子团圆月正半，除百零五便得知。
正好等价于现代同余理论： 23 ≡ 2 * 70 + 3 * 21 + 2 * 15 （mod 105）
原问题为求解如下同余方程组：
x ≡ 2(mod 3）
x ≡ 3(mod 5）
x ≡ 2(mod 7)

定理介绍和证明

自然数 m₁，m₂，…m_k两两互质，并记N = m₁ * m₂* … * m_k ,则同余方程组：
x ≡ b₁ (mod m₁)
x ≡ b₂ (mod m₂)
…
…
x ≡ b_k (mod m_k)
在模N同余的意义下有唯一解，N其实是所有模数的最小公倍数LCM。
方程组的最小正整数解： x₀ = ((b₁ X₁ + b₂ X₂ +… + b_K X_K)(mod N) + N) mod N

证明：
模数两两互素，可令 x = (N/m_i） * y，方程组等价于解同余方程：
(N/m_i） * y ≡ 1（mod m_i ），由裴蜀定理，此方程一定有解，可用扩展欧几里得算法求得一个y_i。让x_i = （N/m_i） * y_i，则方程组的解 x₀ = (b₁ X₁ + b₂ X₂ +… + b_K X_K)(mod N)，在模N意义下唯一。

例题

求解模线性方程组，形如：
a ≡ B[1](mod W[1])
a ≡ B[2](mod W[2])
…
…
a ≡ B[n](mod W[n])
W,B已知，W[i]>0 且所有W互质，求a （最小正整数解）。

代码实现

#include <iostream>
using namespace std;
//中国剩余定理
//所有模数互质 
typedef long long ll;

int W[51],B[51];

ll exgcd(ll a, ll b, ll &x, ll &y)//扩欧算法 
{
	if (!b)
    {
        x = 1,y = 0;
        return a;
    }
    ll gcd = exgcd(b, a % b, x, y);
    
	ll tmp = x;
	x = y;
	y = tmp - a / b * y; 
    return gcd;
}
int China(int * w, int * b, int k)
{
	ll x, y , a = 0;
	int m;
	ll N = 1;
	for(int i = 1; i < k + 1; i ++) N *= w[i];
	
	for(int i = 1; i < k + 1; i ++)
	{
		m = N / w[i];
		exgcd(w[i],m,x,y);
		a = (a + y *m * b[i]) % N;
	}
	return (a + N) % N;
}
int main()
{
  int n;
  cin >> n;//方程数

  for(int i = 1; i <= n; i ++)
  	cin >> W[i] >> B[i];
  
  int res = China(W,B,n); 
  cout << res <<endl;
  
  return 0;
}

解一般意义同余方程通用方案（模数不一定两两互素）

x ≡ a₁ (mod m₁)
x ≡ a₂ (mod m₂)
…
…
x ≡ a_k (mod m_k)
以上模线性方程组，模数两两不一定互质

求解方案论述

模数不互质时，不能使用中国剩余定理求解。
可以用数学方法推导出求解方案。
取前两个模方程，变形为：
x mod m₁ = a₁
x mod m₂ = a₂
先求这两组方程的解x^’ 。
前两方程再变形为（其实后k - 2 个模方程都可以变形为如下等式）：
x = k₁m₁ + a₁
x = k₂m₂ + a₂
进一步推导：k₁m₁ + a₁ = k₂m₂ + a₂
k₁m₁ - k₂m₂ = a₂ - a₁ <= => k₁m₁ + k₂m₂ = a₂ - a₁
此线性同余方程有整数解的充要条件是 gcd(m1,m2) | a2 - a1 , 如无解，则整个模线性方程组无解。
用扩展欧几里得定理求得一特解k₁ 。
可得系数k1的通解为 k₁ + t * m₂ / gcd(m₁,m₂) ; t是任意整数,令 d = gcd(m₁,m₂)，c = m₂ / gcd(m₁,m₂) 显然，k1的最小正整数解为：（k1 % c + c）% c
前两个模方程的特解为： x = k₁m₁ + a₁
前两个模方程的通解为：x= （ k₁ + t * m₂ / gcd(m₁,m₂)） m₁ + a₁ = （ k₁ + t * c） m₁ + a₁
进一步整理，通解为：x = （m₁*m₂/d)t + k₁m₁+a₁ ,m₁*m₂/d是m1和m2的LCM。
如果只有两个模线性方程组，那么 (（ k₁m₁+a₁）% LCM +LCM) % LCM 即为所求最小非负整数解。
如果后续有更多的模线性方程，则可以令 m₁ = m₁*m₂/d ； a₁ = k₁m₁+a₁ ，得到新方程等式： x = m₁ t + a₁ = m₁ k₁+ a₁
新方程和第三个方程继续组成方程组：
x = k₁m₁ + a₁ 新
x = k₃m₃ + a₃
循环求解，直至最后第k号方程，循环结束后，
(（ k₁m₁+a₁）% LCM +LCM) % LCM 即为所求最小非负整数解，此刻，LCM = [m₁,m₂, … m_k]

例题

给定 2n 个整数 a1,a2,…,an 和 m1,m2,…,mn，求一个最小的非负整数 x，满足 ∀i∈[1,n],x≡mi(mod ai)。
输入格式
第 1行包含整数 n。
第 2…n+1行：每 i+1 行包含两个整数 ai 和 mi，数之间用空格隔开。
输出格式
输出最小非负整数 x，如果 x 不存在，则输出 −1。
如果存在 x，则数据保证 x 一定在 64位整数范围内。
数据范围
1≤ai≤2^31−1,0≤mi<ai 1≤n≤25
输入样例：
5
35 30
22 2
16 5
28 23
32 18
输出样例：
-1
输入样例：
10
64 52
82 80
81 20
43 33
100 32
94 58
38 30
60 32
76 68
66 50
输出样例：
365469177332
输入样例：
5
8 5
18 9
10 7
14 9
3 0
输出样例：
2277

代码实现

#include <iostream>
using namespace std;

typedef long long ll;

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;
}
int main()
{
  int n;
  cin >> n;
  
  bool flag = true;
  ll a1, m1; 
  cin >> a1 >> m1;
  
  for(int i = 0; i < n - 1; i ++)
  {
  	ll a2, m2;
  	cin >> a2 >> m2;
  	
  	ll k1, k2;
  	ll d = exgcd(a1,a2,k1,k2);
  	if((m2 - m1) % d)
  	{
  		flag = false;
  		break;
	}
	  k1 *= (m2 - m1) / d;
	  k1 = (k1 % a2 +a2) % a2;//adjust to minimal
	  
	  m1 = a1 * k1 + m1;//特解 
	  a1 = a1 / d * a2;// X = A1 K + M1 
	  
  }
  if(flag)
  {
  	cout << (m1 % a1 + a1) % a1 <<endl;
  }
  else puts("-1");
  return 0;
}