解模线性方程组 非互质中国剩余定理

首先咱们得感谢KIDx大神给出这样的解法。

这里是我所学习这个算法的地方:http://972169909-qq-com.iteye.com/blog/1266328 。

我将对这个算法进行一定的总结与梳理,以及小地方的修正。

 

 

今有物不知其数,三三数之余二;五五数之余三;七七数之余二。问物几何?

 

这是经典的孙子定理。我们注意到其中的模数都是互质的,这样可以让我们进行传统孙子定理中的转化与合并。

但是如果遇到不是互质的模线性方程组我们要怎么办呢?

【主要使用手段:合并方程】

  利用一定手段,不断的合并方程,让方程组合并成最终的一个方程,从而求解。

【合并使用手段】:

  首先我们先列出一个模线性方程组:

  Ans ≡ b1 (Mod mod1)

  Ans ≡ b2 (Mod mod2)

  ......

  Ans ≡ bi (Mod modi)

 

  ……

  Ans ≡ bn (Mod modn)

 

  其中mod1,mod2....modn,表示模数,它们不一定是互质的。

  我将示范如何合并前两行,其余方程的合并和它完全相同。

  Ans ≡ b1 (Mod mod1)  ->  Ans = k1*mod1 + b1

  Ans ≡ b2 (Mod mod2)  ->  Ans = k2*mod2 + b2                 (其中k1,k2为满足条件的某整数,具体的值目前未知)

  所以我们联系上面两式得到   k1*mod1 + b1 = k2*mod2 + b2   发现这里有两个变量,可以想到利用扩展欧几里得求解,所以我们把它稍稍变化:

  k1*mod1 - k2*mod2 = b2-b1 【Wo,这不是ax+by=z的标准格式么!】这样我们就可使用扩展欧几里得算法求出一个满足条件的k1了。【请注意:这是伏笔

  但是我们还需要继续转化,仍然是这个式子 k1*mod1 + b1 = k2*mod2 + b2  

  设d=gcd(mod1,mod2)。

  观察这个等式可以发现k1*mod1可以被d整除,k2*mod2也可以被d整除,所以b2-b1也一定可以被d整除(如果不能整除就不符合k1,k2整数的条件了)。

  于是我们放心了,等式两边同时除以d,得到:

  k1*mod1/d + (b1-b2)/d = k2*mod2/d

  稍稍变一下型-> k1*mod1/d + (b1-b2)/d = (mod2/d)*k2。

  从上面的等式可以再次变成模性方程:

  k1*mod1/d ≡ (b2-b1)/d  (Mod mod2/d)   --模运算满足乘法--> 两边同时乘d,得到  k1*mod1 ≡ (b2-b1)  (Mod mod2/d)

  (原文中将这里的两边同时除以mod1,但是模运算不满足除法,所以不能进行那样的转化

  令K为k1的最小整数解【你问我为什么这样取?为了后面好算...】(如果无解....那就整个都无解了)

  则满足:k1 ≡ K  (Mod mod2/d )  ->  k1 = K + (mod2/d)*C  (C是随意的一个整数,当C=0,则k1就是它能取到的最小整数解)

  将 k1 = K + (mod2/d)*C 带入式子Ans = k1*mod1 + b1,得到 Ans =  ( K + (mod2/d)*C )*mod1 + b1 我们再展开得到:

  Ans =  K*mod1 + (mod2/d)*C*mod1 + b1

  Ans ≡  K*mod1 + b1 (Mod  mod2*mod1/d)  这就是我们所需的式子了。

  可问题还有一个:K怎么求,当然就是k1的最小正整数解了(求的方法之前埋下了伏笔)。这样我们就可使用扩展欧几里得算法求出k1的最小整数解了。

  接下来怎么办?

  将b1= K*mod1 + b1,mod1=mod2*mod1/d,然后一个个地合并就可以了。

  

最后再附上代码:

 1 /*
 2   Problem : 解模线性方程组  CRT(孙子定理) 模数不互质 
 3   Author : Robert Yuan
 4 */
 5 
 6 #include <cstdio>
 7 #include <cstring>
 8 #include <cstdlib>
 9 #include <algorithm>
10 
11 using namespace std;
12 
13 const int maxn=10010;
14 
15 int n;
16 long long x,y,GCD;
17 long long mod[maxn],b[maxn];
18 
19 void exgcd(long long a,long long b){                //扩展欧几里得求模性方程 
20     if(!b){GCD=a,x=1,y=0;return;}
21     exgcd(b,a%b);
22     long long t=x;
23     x=y;y=t-a/b*x;
24 }
25 
26 int CRT(){
27     long long mod1=mod[1],b1=b[1],mod2,b2,bb,K;        //变量名和上文介绍的意义相同 
28     for(int i=2;i<=n;i++){
29         mod2=mod[i],b2=b[i];bb=b2-b1;                 
30         exgcd(mod1,mod2);
31         if(bb%GCD)    return -1;                        //判断无解 
32         mod2/=GCD;bb/=GCD;                            //没办法,扩展欧几里得打习惯了,所以下文中的 mod2其实是 mod2/GCD 
33         if(mod2<0) mod2=-mod2;                        //扩展欧几里得中要保证最后取摸的一定是正数 
34         K=((x*bb)%mod2+mod2)%mod2;                    //这样可以求出最小正整数解...相信你一定能想出来为什么... 
35         b1+=mod1*K;                                    //更新新的 b1 & mod1 
36         mod1*=mod2;
37     }
38     return b1;
39 }
40 
41 int main(){
42 #ifndef ONLINE_JUDGE
43     freopen("x.in","r",stdin);
44 #endif
45     scanf("%d",&n);
46     for(int i=1;i<=n;i++)
47         scanf("%d%d",&mod[i],&b[i]);
48     long long flag=CRT();
49     if(flag<0)
50         printf("NO ANSWER!");
51     else
52         printf("%lld",flag);
53     return 0;
54 }
View Code

 

  

posted @ 2015-07-14 21:26  诚叙  阅读(1181)  评论(0编辑  收藏  举报