模線性方程組

描述：

給定了n組除數m[i]和餘數r[i]，通過這n組(m[i],r[i])求解一個x，使得x mod m[i] = r[i]。

演算法思想：

一開始就直接求解多個方程不是太容易，我們從n=2開始遞推：

已知：

x mod m[1] = r[1]
x mod m[2] = r[2]

根據這兩個式子，我們存在兩個整數k[1],k[2]：

x = m[1] * k[1] + r[1]
x = m[2] * k[2] + r[2]

由於兩個值相等，因此我們有：

	m[1] * k[1] + r[1] = m[2] * k[2] + r[2]
=>	m[1] * k[1] - m[2] * k[2] = r[2] - r[1]

由於m[1],m[2],r[1],r[2]都是常數，若令A=m[1],B=m[2],C=r[2]-r[1],x=k[1],y=k[2]，則上式變為：Ax + By = C。

是不是覺得特別眼熟。

沒錯，這就是我們之前講過的擴充套件歐幾里德。

我們可以先通過gcd(m[1], m[2])能否整除r[2]-r[1]來判定是否存在解。

假設存在解，則我們通過擴充套件歐幾里德求解出k[1],k[2]。

再把k[1]代入x = m[1] * k[1] + r[1]，就可以求解出x。

同時我們將這個x作為特解，可以擴展出一個解系：

X = x + k*lcm(m[1], m[2]) k為整數

lcm(a,b)表示a和b的最小公倍數。其求解公式為lcm(a,b)=a*b/gcd(a,b)。

將其改變形式為：

X mod lcm(m[1], m[2]) = x。

令M = lcm(m[1], m[2]), R = x，則有新的模方程X mod M = R。

此時，可以發現我們將x mod m[1] = r[1]，x mod m[2] = r[2]合併為了一個式子X mod lcm(m[1], m[2]) = x。滿足後者的X一定滿足前兩個式子。

每兩個式子都可以通過該方法化簡為一個式子。那麼我們只要重複進行這個操作，就可以將n個方程組化簡為一個方程，並且求出一個最後的解了。

將其寫做虛擬碼為：

M = m[1], R = r[1]
For i = 2 .. N 
	d = gcd(M, m[i])
	c = r[i] - R
	If (c mod d) Then	// 無解的情況
		Return -1
	End If
	(k1, k2) = extend_gcd(M / d, m[i] / d)	// 擴充套件歐幾里德計算k1,k2
	k1 = (c / d * k1) mod (m[i] / d)	// 擴充套件解系
	R = R + k1 * M		// 計算x = m[1] * k[1] + r[1]
	M = M / d * m[i] 	// 求解lcm(M, m[i])
	R %= M 			// 求解合併後的新R，同時讓R最小
End For		
If (R < 0) Then 
	R = R + M
End If
Return R

AC程式碼：

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;

typedef long long ll;
int n;
ll m[1005],r[1005];

ll gcd(ll a, ll b){
    if(b == 0)
        return a;
    return gcd(b,a%b);
}

ll extend_gcd(ll a, ll b, ll &x, ll &y){
    if(a==0 && b==0) return -1;// 無最大公約數
    if(b == 0){x = 1; y = 0; return a;}
    ll d = extend_gcd(b,a%b,y,x);
    y -= a/b*x;
    return d;
}

ll solve(){
    ll M = m[0],R = r[0];
    for(int i = 1; i < n; ++i){
        ll d = gcd(M,m[i]);
        ll c = r[i]-R;
        if(c % d)
            return -1;// 無解的情況
        ll k1,k2;
        extend_gcd(M/d,m[i]/d,k1,k2);// 擴充套件歐幾里德計算k1,k2
        k1 = (c/d*k1)%(m[i]/d);// 擴充套件解系
        R += k1*M;// 計算x = m[1] * k[1] + r[1]
        M = M/d*m[i];// 求解lcm(M, m[i])
        R %= M;// 求解合併後的新R，同時讓R最小
    }
    if(R < 0)
        R += M;
    return R;
}

int main(){
    while(~scanf("%d",&n)){
        for(int i = 0; i < n; ++i){
            scanf("%d%d",&m[i],&r[i]);
        }
        printf("%lld\n",solve());
    }
    return 0;
}