zoukankan      html  css  js  c++  java
  • 拓展中国剩余定理(exCRT)摘要

    清除一个误区

    虽然中国剩余定理和拓展中国剩余定理只差两个字,但他俩的解法相差十万八千里,所以会不会CRT无所谓

    用途

    求类似$$egin{cases}x equiv b_{1}pmod{a_{1}} \x equiv b_{2}pmod{a_{2}} \...\x equiv b_{n}pmod{a_{n}} \ end{cases}$$的线性同余方程组的解

    具体过程

    假设现在我们只有两个同余方程$$x equiv b_{1}pmod{a_{1}}$$ $$x equiv b_{2}pmod{a_{2}}$$

    改写它们得$$x=a_{1}k_{1}+b_{1}                  ①$$ $$x=a_{2}k_{2}+b_{2}$$

    联立$$a_{1}k_{1}+b_{1}=a_{2}k_{2}+b_{2}$$

    所以$$a_{1}k_{1}=b_{2}-b_{1}+a_{2}k_{2}$$

    记$$c=gcd(a_{1},a_{2})$$ $$d_{1}=frac{a_{1}}{c}$$ $$d_{2}=frac{a_{2}}{c}$$

    两边同除以$c$得$$d_{1}k_{1}=frac{b_{2}-b_{1}}{c}+d_{2}k_{2}$$

    即$$d_{1}k_{1} equiv frac{b_{2}-b_{1}}{c}\%d_{2}pmod{d_{2}}$$

    显然,方程组有解的条件是$c|(b_{2}-b_{1})$

    为什么要在这个时候取模$\%d_{2}$?因为之后取模的时候模数和被模数已经不互质了

    记$d_{1}$在模$d_{2}$意义下的逆元为$p$,两边同乘$p$得$$k_{1} equiv frac{p(b_{2}-b_{1})}{c}\%d_{2} pmod{d_{2}}$$

    即$$k_{1} = frac{p(b_{2}-b_{1})}{c}\%d_{2}+d_{2}*y$$

    带入①式得$$x=frac{p(b_{2}-b_{1})}{c}\%d_{2}*a_{1}+d_{2}a_{1}y+b_{1}$$

    即$$x equiv frac{p(b_{2}-b_{1})}{c}\%d_{2}*a_{1}+b_{1}pmod{a_{1}d_{2}}$$

    注意,上式中的$a_{1}$万万不可乘到$\%d_{2}$前面

    所以我们就弄出了一个新的同余方程。这个方程也可以写成$$x equiv b pmod a$$其中$$b=(inv(frac{a_{1}}{(a_{1},a_{2})},frac{a_{2}}{(a_{1},a_{2})})*frac{b_{2}-b_{1}}{(a_{1},a_{2})}\%frac{a_{2}}{(a_{1},a_{2})}*a_{1}+b_{1})$$ $$a=frac{a_{1}a_{2}}{(a_{1},a_{2})}$$

    其中$inv(m,n)$表示$m$在模$n$意义下的逆元,$(m,n)$表示$m$和$n$的最大公约数

    然后把新求出的方程和后面的方程联立,迭代求解,直到只剩一个方程

    精度问题

    毫无疑问的是,在求解方程组的时候很容易爆$long long$,那该怎么办呢?

    __int128 快(龟)速(速)乘(乘)是个不错的主意。

    何为快速乘?

    就是一种和快速幂类似的东西,利用二进制在log时间内求出两个数相乘对另一个数取模的结果并且不会爆long long

    但要注意的是,千万不要让负数出现在快速乘里面,尤其是被拆分的那个乘数,会死循环的

    1 int mul(int x,int k,int mod) {
    2     ll ans=0;
    3     while(k) {
    4         if(k&1)(ans+=x)%=mod;
    5         k>>=1;
    6         (x<<=1)%=mod;
    7     }
    8     return ans;
    9 }
    快速乘

    题目

    洛谷【模板】扩展中国剩余定理(EXCRT)

     1 #include<bits/stdc++.h>
     2 using namespace std;
     3 #define INF 0x7fffffff
     4 #define ME 0x7f
     5 #define FO(s) freopen(s".in","r",stdin);freopen(s".out","w",stdout)
     6 #define fui(i,a,b,c) for(int i=(a);i<=(b);i+=(c))
     7 #define fdi(i,a,b,c) for(int i=(a);i>=(b);i-=(c))
     8 #define fel(i,a) for(register int i=hd[a];i;i=dg[i].nxt)
     9 #define ll long long
    10 #define int long long
    11 #define MEM(a,b) memset(a,b,sizeof(a))
    12 #define maxn 100010
    13 ll n;
    14 ll a[maxn],b[maxn],a0,b0;
    15 template<class T>
    16 inline T read(T &n){
    17     n=0;int t=1;double x=10;char ch;
    18     for(ch=getchar();!isdigit(ch)&&ch!='-';ch=getchar());(ch=='-')?t=-1:n=ch-'0';
    19     for(ch=getchar();isdigit(ch);ch=getchar()) n=n*10+ch-'0';
    20     if(ch=='.') for(ch=getchar();isdigit(ch);ch=getchar()) n+=(ch-'0')/x,x*=10;
    21     return (n*=t);
    22 }
    23 template<class T>
    24 T write(T n){
    25     if(n<0) putchar('-'),n=-n;
    26     if(n>=10) write(n/10);putchar(n%10+'0');return n;
    27 }
    28 template<class T>
    29 T writeln(T n){write(n);putchar('
    ');return n;}
    30 int exgcd(int a,int b,int &x,int &y){
    31     if(!b)return a+((x=1)*(y=0));
    32     int xx,yx,z;z=exgcd(b,a%b,xx,yx);
    33     x=yx,y=xx-a/b*yx;return z;
    34 }
    35 int gcd(int a,int b){if(!b)return a;return gcd(b,a%b);}
    36 int mul(int x,int k,int mod){
    37     ll ans=0;
    38     while(k){
    39         if(k&1)(ans+=x)%=mod;
    40         k>>=1;
    41         (x<<=1)%=mod;
    42     }
    43     return ans;
    44 }
    45 signed main(){
    46     read(n);
    47     fui(i,1,n,1)read(a[i]),read(b[i]);
    48     a0=a[1],b0=b[1];
    49     fui(i,2,n,1){
    50         int c=gcd(a0,a[i]),d1=a0/c,d2=a[i]/c,p,x,y;
    51         int a1=a0,a2=a[i],b1=b0,b2=b[i];
    52         if(((b2-b1)/c*c)^(b2-b1))return 0*writeln(-1);
    53         exgcd(d1,d2,p,y);((p%=d2)+=d2)%=d2;
    54         x=mul(p,(((b2-b1)/c)%d2+d2)%d2,d2)*a1+b1;
    55         y=a1*d2;b0=x,a0=y;((b0%=a0)+=a0)%=a0;
    56     }
    57     int x,y;exgcd(1,a0,x,y);
    58     ((x%=a0)+=a0)%=a0;
    59     x=mul(x,b0,a0);
    60     writeln(x);
    61     return 0;
    62 }
    AC代码

     好像还有NOI2018屠龙勇士?23333先等我去A掉它

    作者:A星际穿越
    我的博客写得这么烂,应该不会有人想转载的吧
    如果要转载的话,请在文章显眼处标明作者和出处谢谢
  • 相关阅读:
    1024X768大图 (Wallpaper)
    (Mike Lynch)Application of linear weight neural networks to recognition of hand print characters
    瞬间模糊搜索1000万基本句型的语言算法
    单核与双核的竞争 INTEL P4 670对抗820
    FlashFTP工具的自动缓存服务器目录的功能
    LDAP over SSL (LDAPS) Certificate
    Restart the domain controller in Directory Services Restore Mode Remotely
    How do I install Active Directory on my Windows Server 2003 server?
    指针与指针变量(转)
    How to enable LDAP over SSL with a thirdparty certification authority
  • 原文地址:https://www.cnblogs.com/Axjcy/p/9801119.html
Copyright © 2011-2022 走看看