标签:
首先咱们得感谢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 }
标签:
原文地址:http://www.cnblogs.com/Robert-Yuan/p/4646632.html