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 102

int a[MAX];
int d[MAX];

int X[MAX];

//产生 a mod b == x mod b 的最小非负数 x

int 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;
}