求解线性同余方程组

求解线性同余方程组

线性同余方程组就是形如这样的方程组：

$$\{\begin{array}{l}{a}_{1}x\equiv b_1(mod{m}_1)\\ a_{2}x\equiv b_2(mod{m}_2)\\ ⋮\\ a_{n}x\equiv b_n(mod{m}_n)\end{array}$$

与扩展中国剩余定理的证明类似，首先找出一个形如这样的方程：$x\equiv b_0(mod{m}_0)$

可以化为：$x=b_0+m_0{k}_0$

第一步可以假定 $b_1=0,m_1=1$ ，代表解是所有的正整数。

之后与方程组的第一个方程联立：

$$\{\begin{array}{l}x\equiv b_0(mod{m}_0)\\ a_{1}x\equiv b_1(mod{m}_1)\end{array}$$

把一式带入二式得：$$a_1(b_0+m_0{k}_0)\equiv b_1(mod{m}_1)$$

$$a_1{b}_0+a_1{m}_0{k}_0=b_1+m_1{k}_1$$

$$a_1{m}_0{k}_0=b_1-a_1{b}_0+m_1{k}_1$$

此时令$(a,b)$表示$gcd(a,b)$，则有：

$$\frac{a_1{m}_1}{(a_1{m}_0,m_1)}{k}_0=\frac{b_1-a_1{b}_0}{(a_1{m}_0,m_1)}+\frac{m_1}{(a_1{m}_0,m_1)}{k}_1$$

（此时$(a_1{m}_0,m_1)$一定要能够整除$b_1-a_1{b}_0$，否则方程无解）

所以有：$$\frac{a_1{m}_0}{(a_1{m}_0,m_1)}{k}_0\equiv \frac{b_1-a_1{b}_0}{(a_1{m}_0,m_1)}(mod\frac{m_1}{(a_1{m}_0,m_1)})$$

$$k_0\equiv \frac{b_1-a_1{b}_0}{(a_1{m}_0,m_1)}inv(\frac{a_1{m}_0}{(a_1{m}_0,m_1)},\frac{m_1}{(a_1{m}_0,m_1)})(mod\frac{m_1}{(a_1{m}_0,m_1)})$$

$$k_0=\frac{b_1-a_1{b}_0}{(a_1{m}_0,m_1)}inv(\frac{a_1{m}_0}{(a_1{m}_0,m_1)},\frac{m_1}{(a_1{m}_0,m_1)})+y\frac{m_1}{(a_1{m}_0,m_1)}$$

代入高亮的上式：

$$x=\frac{b_1-a_1{b}_0}{(a_1{m}_0,m_1)}inv(\frac{a_1{m}_0}{(a_1{m}_0,m_1)},\frac{m_1}{(a_1{m}_0,m_1)})m_0+y\frac{m_0{m}_1}{(a_1{m}_0,m_1)}+b_0$$

$$x\equiv \frac{b_1-a_1{b}_0}{(a_1{m}_0,m_1)}inv(\frac{a_1{m}_0}{(a_1{m}_0,m_1)},\frac{m_1}{(a_1{m}_0,m_1)})m_0+b_0(mod\frac{m_0{m}_1}{(a_1{m}_0,m_1)})$$

此时可以将方程看作：$$x\equiv b(modm)$$

其中$$b=(\frac{b_1-a_1{b}_0}{(a_1{m}_0,m_1)}inv(\frac{a_1{m}_0}{(a_1{m}_0,m_1)},\frac{m_1}{(a_1{m}_0,m_1)}))\%\frac{m_1}{(a_1{m}_0,m_1)}{m}_0+b_0$$

$$m=\frac{m_0{m}_1}{(a_1{m}_0,m_1)}$$

到这里，式子里所有的值我们都知道了，并且得到了一个新的同余式。再把新的和其他的联立，依此循环，最终求出通解。

代码

//返回一个（b，m）的数对
pair<int, int>linear_congruence(const vector<int>&A, const vector<int>&B, const vector<int>&M) {
	//由于最开始没有任何限制，所以先把解设为表示所有整教的 x三0（mod 1）
	int b = 0, m = 1;
	for (int i = 0; i < A.size(); i++) {
		int a = A[i] * m, b = B[i] - A[i] * b, d = gcd(M[i], a);
		if (b%d != 0) return make_pair(0, -1);//无解
		int t = b / d * mod_inverse(a / d, M[i] / d) % (M[i] / d);
		b = b + m * t;
		m *= M[i] / d;
	}
	return make_pair(x % m, m);
}