zoukankan      html  css  js  c++  java
  • 解模线性方程组 非互质中国剩余定理

    首先咱们得感谢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

      

  • 相关阅读:
    使用电脑模拟微信内置浏览器
    手机QQ浏览器属于代理服务器吗?
    艾伟:[你必须知道的.NET]第三十二回,深入.NET 4.0之,Tuple一二 狼人:
    艾伟:Silverlight 2.0在IE6 SP2上的虚线边框问题 狼人:
    艾伟:ASP.NET安全问题--Forms验证(后篇)--实战篇 狼人:
    艾伟:基于.NET平台的Windows编程实战(四)—— 数据库操作类的编写 狼人:
    艾伟:WMGPS开发 狼人:
    艾伟:基于.NET平台的Windows编程实战(一)——前言 狼人:
    艾伟:.NET,你忘记了么?(八) 从dynamic到特性误用 狼人:
    艾伟:小巧优美的ORM框架doodads入门指南[转载] 狼人:
  • 原文地址:https://www.cnblogs.com/Robert-Yuan/p/4646632.html
Copyright © 2011-2022 走看看