[关闭]
@Junlier 2018-10-03T19:45:56.000000Z 字数 2373 阅读 3411

中国剩余定理(CRT) & 扩展中国剩余定理(ExCRT)总结

数学方法——数论
阅读体验:https://zybuluo.com/Junlier/note/1300035

前置浅讲

前置知识点:
这两个东西都是用来解同余方程组
形如

给定,求解唯一解满足上述方程组
PS:个人认为自己的讲得好一些

中国剩余定理(CRT)

心理准备:这里我觉得自己讲得不是很清晰,有点说不清的感觉
在上述方程中,存在一种特殊情况,即全部互质
有什么用呢,用处就是可以用中国剩余定理(孙子定理)
首先,古人告诉我们:
解上面那个方程相当于对于每一个

变成,其他的变成的解
然后答案就是
也就是解每一个形如下面的方程组的解在乘上

别问我为什么,我不知道
那么考虑怎么求呢?

显然解必须是的倍数
那么方程变为


求解即可,再加进答案

  1. lst CRT()
  2. {
  3. lst tot=1,Ans=0;
  4. for(int i=1;i<=n;++i)tot*=W[i];
  5. for(int i=1;i<=n;++i)
  6. {
  7. lst now=tot/W[i],x,y;
  8. Exgcd(now,W[i],x,y);//不需要我讲吧!
  9. x=(x%W[i]+W[i])%W[i];//这个取膜貌似很关键诶
  10. Ans=(Ans+Mult(Mult(x,now,tot),B[i],tot)+tot)%tot;//Mult是快速乘
  11. }return Ans>=0?Ans:Ans+tot;
  12. }

扩展中国剩余定理(ExCRT)

再把方程组放一遍:

这个时候我们的不互质了
那么我们考虑对所有的方程一个一个求解
假设我们求解到了第个方程
前面的方程组解出来答案是
那我们是不是可以把之前求解的答案看做一个这样的同余方程:
其中是前的最小公倍数,使我们想得到的新解
如果看不出请补一下同余方程。。。
很显然这些都是已知量了吧
那我们为了求出前个方程的解就相当于要解出下面这个方程组了对不对

表示Ans,表示M,表示表示表示新答案
可以化为
为了使用扩展欧几里德我们设
那么我们会解出一个,而真正的如下
我们再回代就得到了新解
有点凌乱请认真看(trust me)!
模板题:洛谷P4777 【模板】扩展中国剩余定理(EXCRT)
代码:

  1. #include<bits/stdc++.h>
  2. #define lst long long
  3. #define ldb double
  4. #define N 100050
  5. using namespace std;
  6. const int Inf=1e9;
  7. lst read()
  8. {
  9. lst s=0,m=0;char ch=getchar();
  10. while(!isdigit(ch)){if(ch=='-')m=1;ch=getchar();}
  11. while( isdigit(ch))s=(s<<3)+(s<<1)+(ch^48),ch=getchar();
  12. return m?-s:s;
  13. }
  14. lst n,Ans,x,y,M;
  15. lst A[N],B[N];
  16. lst qmul(lst x,lst y,lst p)
  17. {
  18. lst ret=0;
  19. while(y)
  20. {
  21. if(y&1)ret=(ret+x)%p;
  22. x=(x+x)%p;y>>=1;
  23. }return ret;
  24. }
  25. lst Exgcd(lst a,lst b,lst &x,lst &y)
  26. {
  27. if(!b){x=1,y=0;return a;}
  28. lst ss=Exgcd(b,a%b,x,y),t;
  29. t=x,x=y,y=t-(a/b)*y; return ss;
  30. }
  31. int main()
  32. {
  33. n=read();
  34. for(int i=1;i<=n;++i)
  35. B[i]=read(),A[i]=read();
  36. Ans=A[1],M=B[1];
  37. //根据上面的详解一步对应一行
  38. for(int i=2;i<=n;++i)
  39. {
  40. lst Get=((A[i]-Ans)%B[i]+B[i])%B[i];
  41. lst GCD=Exgcd(M,B[i],x,y);
  42. x=qmul(x,Get/GCD,B[i]);//qmul是龟速乘
  43. Ans+=M*x;
  44. M*=B[i]/GCD;//这里是更新辅助下一次计算
  45. Ans=(Ans+M)%M;
  46. }printf("%lld\n",Ans);
  47. return 0;
  48. }
添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注