中国剩余定理

来源：互联网发布：购买域名的网站编辑：程序博客网时间：2024/06/03 09:19

参考自：http://blog.csdn.net/acdreamers/article/details/8050018

中国剩余定理（CRT）的表述如下

设正整数两两互素，则同余方程组

有整数解。并且在模下的解是唯一的，解为

其中，而为模的逆元。

代码：

int CRT(int a[],int m[],int n)  {      int M = 1;      int ans = 0;      for(int i=1; i<=n; i++)          M *= m[i];      for(int i=1; i<=n; i++)      {          int x, y;          int Mi = M / m[i];          extend_Euclid(Mi, m[i], x, y);          ans = (ans + Mi * x * a[i]) % M;      }      if(ans < 0) ans += M;      return ans;  }

习题：poj1006
代码：

#include <iostream>  #include <string.h>  #include <stdio.h>    using namespace std;    int a[4], m[4];    void extend_Euclid(int a, int b, int &x, int &y)  {      if(b == 0)      {          x = 1;          y = 0;          return;      }      extend_Euclid(b, a % b, x, y);      int tmp = x;      x = y;      y = tmp - (a / b) * y;  }    int CRT(int a[],int m[],int n)  {      int M = 1;      int ans = 0;      for(int i=1; i<=n; i++)          M *= m[i];      for(int i=1; i<=n; i++)      {          int x, y;          int Mi = M / m[i];          extend_Euclid(Mi, m[i], x, y);          ans = (ans + Mi * x * a[i]) % M;      }      if(ans < 0) ans += M;      return ans;  }    int main()  {      int p, e, i, d, t = 1;      while(cin>>p>>e>>i>>d)      {          if(p == -1 && e == -1 && i == -1 && d == -1)              break;          a[1] = p;          a[2] = e;          a[3] = i;          m[1] = 23;          m[2] = 28;          m[3] = 33;          int ans = CRT(a, m, 3);          if(ans <= d)              ans += 21252;          cout<<"Case "<<t++<<": the next triple peak occurs in "<<ans - d<<" days."<<endl;      }      return 0;  }

普通的中国剩余定理要求所有的互素，那么如果不互素呢，怎么求解同余方程组？

这种情况就采用两两合并的思想，假设要合并如下两个方程

那么得到

在利用扩展欧几里得算法解出的最小正整数解，再带入

得到后合并为一个方程的结果为

即新方程为X = A mod M,其中A就是上式的x，M为lcm(m1,m2)。

这样一直合并下去，最终可以求得同余方程组的解。

题目：poj2891

#include <iostream>  #include <string.h>  #include <stdio.h>    using namespace std;  typedef long long LL;  const int N = 1005;    LL a[N], m[N];    LL gcd(LL a,LL b)  {      return b? gcd(b, a % b) : a;  }    void extend_Euclid(LL a, LL b, LL &x, LL &y)  {      if(b == 0)      {          x = 1;          y = 0;          return;      }      extend_Euclid(b, a % b, x, y);      LL tmp = x;      x = y;      y = tmp - (a / b) * y;  }    LL Inv(LL a, LL b)  {      LL d = gcd(a, b);      if(d != 1) return -1;      LL x, y;      extend_Euclid(a, b, x, y);      return (x % b + b) % b;  }    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];      LL m1 = m[1];      for(int i=2; i<=n; i++)      {          LL a2 = a[i];          LL m2 = m[i];          LL m3, a3;          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)!=EOF)      {          for(int i=1; i<=n; i++)              scanf("%I64d%I64d",&m[i], &a[i]);          LL ans = CRT(a, m, n);          printf("%I64d\n",ans);      }      return 0;  }

阅读全文

0 0