算法总结之求解模线性方程组

算法总结之求解模线性方程组

1）求解模线性方程 ax = b(mod n)

　　方程ax = b(mod n) -> ax = b + ny ->ax - ny = b

　　-> ax + n (-y) =b 其中a,n,b已知。 可用扩展欧几里得来求解该方程的一组特解。

　　这里给出下列几个定理用来求解方程：

　　1.当且仅当d|b时，方程ax = b(mod n)有解。d=gcd(a,n)

　　2.ax = b(mod n) 或者有d个不同解，或者无解。

　　3.令d=gcd(a,n) 假定对整数x', y', 有d = ax' + ny'， 如果d | b, 则方程ax = b(mod n)有一个解的值为x0, 满足：

　　　　x0=x‘(b/d)(mod n)

　　4.假设方程ax = b(mod n)有解， x0是方程的任意一个解， 则方程对模n恰有d个不同的解， 分别为:

　　　　xi = x0 + i * (n / d), 其中 i = 1,2,3......d - 1

　　根据这4个定理，运用扩展欧几里得算法就能轻易的求出模线性方程的所有解了。

伪代码如下:

MODULAR_LINEAR_EQUATION_SOLVER(a,b,n)
(d,x',y')=EXTENDED_EUCLID(a,n)
if (d|b)
    x0=x'(b/d) mod n
    for i=0 to d-1
    print (x0+i(n/d)) mod n
else
    print "no solutions"

2)求解模线性方程组

　　x = a1(mod m1)

　　x = a2(mod m2)

　　x = a3(mod m3)

　　

　　先求解方程组前两项。 x=m1*k1+a1=m2*k2+a2

　-> m1*k1+m2*(-k2)=a2-a1

　　这个方程可以通过欧几里得求解出最小正整数的k1 则x=m1*k1+a1 显然x为两个方程的最小正整数解。

　　则这两个方程的通解为 X=x+k*LCM(m1,m2) -> X=x(mod LCM(m1,m2)) 就转换成了一个形式相同方程了

　　在通过这个方程和后面的其他方程求解。最终的结果就出来了。

　　以POJ2891为例 贴上代码：

Code：

/*************************************************************************
    > File Name: poj2891.cpp
    > Author: Enumz
    > Mail: 369372123@qq.com
    > Created Time: 2014年10月28日 星期二 02时50分07秒
 ************************************************************************/

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<string>
#include<cstring>
#include<list>
#include<queue>
#include<stack>
#include<map>
#include<set>
#include<algorithm>
#include<cmath>
#include<bitset>
#include<climits>
#define MAXN 100000
#define LL long long
using namespace std;
LL extended_gcd(LL a,LL b,LL &x,LL &y) //返回值为gcd(a,b)
{
    LL ret,tmp;
    if (b==0)
    {
        x=1,y=0;
        return a;
    }
    ret=extended_gcd(b,a%b,x,y);
    tmp=x;
    x=y;
    y=tmp-a/b*y;
    return ret;
}
int main()
{
    LL N;
    while (cin>>N)
    {
        long long a1,m1;
        long long a2,m2;
        cin>>a1>>m1;
        if (N==1)
            printf("%lld
",m1);
        else
        {
            bool flag=0;
            for (int i=2;i<=N;i++)
            {
                cin>>a2>>m2;
                if (flag==1) continue;
                long long x,y;
                LL ret=extended_gcd(a1,a2,x,y);
                if ((m2-m1)%ret!=0)
                    flag=1;
                else
                {
                    long long ans1=(m2-m1)/ret*x;
                    ans1=ans1%(a2/ret);
                    if (ans1<0) ans1+=(a2/ret);
                    m1=ans1*a1+m1;
                    a1=a1*a2/ret;
                }
            }
            if (!flag)
                cout<<m1<<endl;
            else
                cout<<-1<<endl;
        }
    }
    return 0;
}