题目连接:
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;
}