1. 程式人生 > >【POJ2891】【一般模線性方程 模板題】Strange Way to Express Integers

【POJ2891】【一般模線性方程 模板題】Strange Way to Express Integers

題意:

給的模數不滿足互質,求同餘方程組的解

思路:

同樣是求這個東西。。
X mod m1=r1
X mod m2=r2
...
...
...
X mod mn=rn

首先,我們看兩個式子的情況
X mod m1=r1……………………………………………………………(1)
X mod m2=r2……………………………………………………………(2)
則有 
X=m1*k1+r1………………………………………………………………(*)
X=m2*k2+r2
那麼 m1*k1+r1=m2*k2+r2
整理,得
m1*k1-m2*k2=r2-r1
令(a,b,x,y,m)=(m1,m2,k1,k2,r2-r1),原式變成
ax+by=m
熟悉吧?

此時,因為GCD(a,b)=1不一定成立,GCD(a,b) | m 也就不一定成立。所以應該先判 若 GCD(a,b) | m 不成立,則!!!方程無解!!!。
否則,繼續往下。

解出(x,y),將k1=x反代回(*),得到X。
於是X就是這兩個方程的一個特解,通解就是 X'=X+k*LCM(m1,m2)
這個式子再一變形,得 X' mod LCM(m1,m2)=X
這個方程一出來,說明我們實現了(1)(2)兩個方程的合併。
令 M=LCM(m1,m2),R=r2-r1
就可將合併後的方程記為 X mod M = R。

然後,擴充套件到n個方程。
用合併後的方程再來和其他的方程按這樣的方式進行合併,最後就能只剩下一個方程 X mod M=R,其中 M=LCM(m1,m2,...,mn)。
那麼,X便是原模線性方程組的一個特解,通解為 X'=X+k*M。

如果,要得到X的最小正整數解,就還是原來那個方法:

X%=M;
if (X<0) X+=M;

程式碼一:

//不滿足模數兩兩互質
#include <iostream>  
#include <string.h>  
#include <stdio.h>  
#define ll __int64
using namespace std;

const int inf = 0x3f3f3f3f;
const int maxn = 1e5 + 5;
int n;
ll m[maxn], r[maxn];

void extgcd(ll a, ll b, ll &d, ll &x, ll &y){
  if (!b) {d = a, x = 1, y = 0;}
  else{
    extgcd(b, a % b, d, y, x);
    y -= x * (a / b);
  }
}

ll extcrt(ll *m, ll *r, int n){
  ll M = m[1], R = r[1], x, y, d;
  for (int i = 2; i <= n; ++i){
    extgcd(M, m[i], d, x, y);
    if ((r[i] - R) % d) return -1;
    x = (r[i] - R) / d * x % (m[i] / d);
    R += x * M;
    M = M / d * m[i];
    R %= M;
  }
  return R > 0 ? R : R + M;
}

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


普通的中國剩餘定理要求所有的互素,那麼如果不互素呢,怎麼求解同餘方程組?

這種情況就採用兩兩合併的思想,假設要合併如下兩個方程

      

那麼得到

       

在利用擴充套件歐幾里得演算法解出的最小正整數解,再帶入

       

得到後合併為一個方程的結果為

       

這樣一直合併下去,最終可以求得同餘方程組的解。

程式碼二:
#include <iostream>  
#include <string.h>  
#include <stdio.h>  
#define ll __int64
#define rep(i,k,n) for(int i=k;i<=n;i++)
using  namespace  std;

const int N=1010;
ll a[N],m[N];

ll gcd(ll a,ll b){
  return b ? gcd(b, a%b) : a;
}

void extgcd(ll a, ll b, ll &x, ll &y)  {  
  if(b == 0){  
    x = 1; y = 0;  
    return;  
  }  
  extgcd(b, a % b, x, y);  
  ll tmp = x;  
  x = y;  
  y = tmp - (a / b) * y;  
}  
  
ll inv(ll a,ll m){//a模m的逆元
  ll x, y;  
  extgcd(a, m, x, y);//a*x+m*(-k)=1
  return (x % m + m) % m;
}  

bool merge(ll a1, ll m1, ll a2, ll m2, ll &a3, ll &m3){
  ll d = gcd(m1, m2);  
  ll c = a2 - a1;  
  if(c % d) return false;  
  c = (c % m2 + m2) % m2;  
  m1 /= d;  
  m2 /= d;  
  c /= d;  
  c *= inv(m1, m2);  
  c %= m2;  
  c *= m1 * d;  
  c += a1;  
  m3 = m1 * m2 * d;  
  a3 = (c % m3 + m3) % m3;  
  return true;  
}
  
ll CRT(ll a[],ll m[],int n){
  ll a1=a[1],m1=m[1];
  rep(i, 2, n){  
    ll a2=a[i], m2=m[i];  
    ll a3, m3;
    if(!merge(a1, m1, a2, m2, a3, m3))
      return -1;
    a1 = a3, m1 = m3;
  }  
  return (a1 % m1 + m1) % m1;  
} 

int main(){
  int n;
	while(~scanf("%d",&n)){
		rep(i, 1, n)scanf("%I64d%I64d",&m[i], &a[i]);
    ll ans=CRT(a, m, n);
    printf("%I64d\n", ans);
	}
	return 0;
} 

轉載自http://blog.csdn.net/acdreamers/article/details/8050018