spoj 3179 DPEQN 题解

题目连接：

http://www.spoj.pl/problems/DPEQN/

本题是探求多元一次同余方程的解法，虽然本题只要求一个解，但是我们完全可以求出所有解。

我们先来复习一下一元一次同余方程的解法（算法导论中有详细介绍，可以参看）。

一元一次同余方程的形式如下：

我们知道，本方程有解的充分必要条件是 gcd(a,n) | b 。如果有解，解的个数（当然是在 [0,n)范围内 ），为d ，其中d = gcd(a,n)

设 d = gcd(a,n),假定对整数x1 和 y1,有 d = ax1 + ny1(其中x1可以由欧几里德的扩展算法求出)。如果d|b，则方程有一个解

的值为x0,满足： 

x0 = x1(b/d)mod n

那么其所有解怎么求呢，我把算法导论上的伪代码搬来：

MODULAR-LINEAR-EQUATION-SOLVER(a,b,n) (d,x1,y1) <- EXTENDED-EUCLID(a,n) if d|b    then x0 <- x1(b/d) mod n        for i <- 0 to d-1            do print (x0 + i(n/d)) mod n    else print "no solution"

其中，EXTENDED-EUCLID(a,n)是欧几里德的扩展算法，用于求gcd(a,n)，并把其中的线性关系求出，伪代码描述：

EXTENDED-EUCLID(a,n)    if(b=0)        then return (a,1,0)    (d1,x1,y1) <- (d1,y1,x1-a/b*y )    return (d,x,y)

求一元一次的同余方程很简单，只要利用这两个函数就能搞定。

那么如何求多元同余方程的解呢？

形式如下：

可以转换为一元一次的同余方程来做。具体做法如下：

1.首先我们要判断是否有解，类似于一元一次同余方程，多元的有解条件也是：gcd(a1,a2,a3...an,m) | b

2.求解。求解的方法是一次求出gcd。即：d1 = gcd(a1,m), d2 = gcd(a1,a2,m),d3 = gcd(a1,a2,a3,m)...,dn = gcd(a1,a2,a3,...an,m)，这些其实可以作为判断是否有解的中间结果。

3.接着，我们如是来求解，从右至左：

先求解：。解出的xn肯定是所有解中的一个。（读者可以自己证明）

然后令 b1 = an * xn ,其中xn已经求出。b1 就是常数，现在的方程消掉一元，变为了：

我们再求解。同理，解出的xn-1肯定是所有解中的一个.

令b2 = an-1 * xn-1  ,其中xn-1已经求出。b2 就是常数，现在的方程消掉二元，变为了：

。依次类推求解。

最后一次我们只剩下：

解这个一元一次的同余方程，求出x1即可。。综上，就是一组解。如果要求所有的解，只要利用求一元的所有解的解法即可。

可以证明，只要多元同余方程有解。解的个数为：

。。。

下面贴出本题的代码，作为模板用：

#include <iostream>#include <stdio.h>using namespace std;#define MAX 102int a[MAX];int d[MAX];int X[MAX];//产生 a mod b == x mod b 的最小非负数 xint procMod(int a,int b){    if(a % b >=0)    {        return a % b;    }    else    {        return a % b + b;    }}//欧几里德算法推广void exGcd(int a,int b,int &d,int &x,int &y){    if(b == 0)    {        x = 1;        y = 0;        d = a;        return ;    }    exGcd(b,a%b,d,x,y);    int temp = x;    x = y;    y = temp - a/b*y;    return ;}int judge(int n,int b,int m){    int x,y;    for(int i=0; i<n; i++)    {        if(i == 0)        {            exGcd(a[0],m,d[0],x,y);        }        else        {            exGcd(d[i-1],a[i],d[i],x,y);        }    }    if(b%d[n-1] == 0)    {        return 1;    }    else    {        return 0;    }}int main(){    int T;    int n,b,m;    int b2;    #ifndef ONLINE_JUDGE        freopen("in.txt","r",stdin);    #endif    scanf("%d",&T);    while(T--)    {        scanf("%d",&n);        for(int i=0; i<n; i++)        {            scanf("%d",&a[i]);        }        scanf("%d %d",&b,&m);        b = b%m;        for(int i=0; i<n; i++)        {            a[i] = a[i]%m;        }        if(judge(n,b,m) == 1)        {            b2 = b;            for(int i=n; i>0; i--)            {                int x,y;                int temp_d;                int current;                if(i!=n)                {                    current = int((long long)a[i] * X[i] % m);                }                else                {                    current = 0;                }                b = procMod(-current + b ,m);                if(i>1)                {                    exGcd(a[i-1],d[i-2],temp_d,x,y);                    b2 =  b % d[i-2];                    X[i-1] = int((long long)x * (b2/temp_d) % d[i-2]);                    if(X[i-1]<0)                    {                        X[i-1] = X[i-1] + d[i-2];                    }                }                else                {                    exGcd(a[0],m,temp_d,x,y);                    b2 =  b % m;                    X[0] = int((long long)x * (b2/temp_d) % m);                    if(X[0]<0)                    {                        X[0] = X[0] + m;                    }                }            }            for(int i=0;i<n-1;i++)            {                printf("%d ",X[i]);            }            printf("%d\n",X[n-1]);        }        else        {            printf("NO\n");        }    }    return 0;}