• CRT&EXCRT 中国剩余定理及其扩展


    前言:

    中国剩余定理又名孙子定理。因孙子二字歧义,常以段子形式广泛流传。

     

    中国剩余定理并不是很好理解,我也理解了很多次。

    CRT 中国剩余定理

    中国剩余定理,就是一个解同余方程组的算法。

    求满足n个条件的最小的x。

    看起来很麻烦。

    先找一个特殊情况:$m_1,m_2,...m_n$两两互质。

    这个时候,构造$M=m_1*m_2*...m_n$;

    令$M_i=M/m_i$;

    所以,构造$n$个数,其中第$i$个数是除$i$之外的其他所有数的倍数,并且第$i$个数$mod m_i =1$

    即:$M_i x = 1 ( mod m_i ) $求出这样一个x,就求出了 这个数。

    因为$m$之间两两互质,所以对于$n$个这样的方程,$x$本质上就是$M_i$在$m_i$意义下的乘法逆元。

    (不会$exgcd$?左转:EXGCD 扩展欧几里得

    因为互质,一定有解的。

    用扩展欧几里得算就可以。

    同理,构造$n$个数。$b_1,b_2....b_n$

    其中,$b_i=M_i imes x_i$

    那么,因为$b_i  = 1 (mod m_i)$,所以$ b_i * a_i = a_i (mod m_i)$

    那么,原题目中的这个x就是:$x=(a_1 imes b_1+a_2 imes b_2+...+a_n imes b_n) $验证一下,是不是?

    总得来说,

    对于$mi$互质的情况,$x=sum_1^n M_i imes a_i imes inv_i$

    其中,$inv_i$表示,$M_i$在$modspace m_i$意义下的逆元。

    当然,为了保证$x$最小,要让$x$和$lcm$做一些处理。当然,因为互质,所以$lcm$就是$M$了

    x=(x%M+M)%M

    例题:poj1006 生理周期 Biorhythms

    在CRT上的小小变形。注意开始计算的时间d就好了。上述最小的$x$就是$n+d$

    答案$n=(a_1 imes b_1+a_2 imes b_2+a_3 imes b_3-d)%lcm$

    #include<cstdio>
    #include<iostream>
    #include<algorithm>
    #include<cstdlib>
    using namespace std;
    int p,e,l;
    int p1,e1,l1;
    int tt,d;
    void exgcd(int a,int b,int &x,int &y){
        if(b==0){
            x=1,y=0;return;
        }
        exgcd(b,a%b,y,x);
        y-=(a/b)*x;
    }
    int main(){
        exgcd(23*28,33,l1,tt);
        l1=(l1%33+33)%33;
        l1*=23*28;
        
        exgcd(28*33,23,p1,tt);
        p1=(p1%23+23)%23;
        p1*=28*33;
        
        exgcd(23*33,28,e1,tt);
        e1=(e1%28+28)%28;
        e1*=23*33;
        
        int lcm=23*28*33;
        int cnt=0;
        while(1){
            ++cnt;
            scanf("%d%d%d%d",&p,&e,&l,&d);
            if(p==-1) break;
            int op=(p1*p+e1*e+l1*l-d+lcm)%lcm;
            if(op==0) op=lcm;
            printf("Case %d: the next triple peak occurs in %d days.
    ",cnt,op);
        }
        return 0;
    }
    poj1006

    EXCRT 扩展中国剩余定理

    但是,不是所有的方程,$m$都是互质的。

    $m$不是互质的时候,我们的$M_i x = 1(mod m_i) $可能就没有解了。所以挂掉。

    孙子解决不了了。

    但是现代人不是孙子也解决了。

    (你是孙子吗)

    现在,$m_1,m_2,...m_n$之间没有任何关系。

    只考虑两个怎么处理?

    可以得到:$x=a_1+k_1*m_1 ; x=a_2+k_2*m_2$

    所以 $a_1+k_1*m_1 = a_2+k_2*m_2$

    $k_2*m_2-k_1*m_1=a_1-a_2$;

    很像:$a imes x +b imes y=c$

    设$m_1,m_2$ 的$gcd$ 为 $g$

    设$a_1-a_2=c$

    当$c$不是$g$的倍数的时候,那就完了。($exgcd$无解情况)

    如果是,就用$exgcd$求出$k_2 imes m_2+k_1 imes m_1=gcd(m_1,m_2)$的$k1$

    因为$c$是$g$的倍数,所以,两边同时乘上$c/g$,即$k_1$乘上$c/g$

    就得到了$k_2 imes m_2+k_1 imes m_1=c$的解$k_1$。

    当然,最好$k_1$ 再模一下$m_2$ ($k_1$,$k_2$做出调整),防止爆$long long$

    然后可以反推x,

    但是注意,我们列的原方程是:$k_2 imes m_2-k_1 imes m_1=c$

    差一个符号,所以$k_1$ 实际上是 $-k_1$

    $x=-k_1 imes m_1+a_1$;

    这样就求出了$x$。(可以把这个$x0$转化成最小的非负数解)

    这个$x$符合第一个方程,也符合第二个方程。设这个$x$为$x_0$

    所以,可以得到通解是:$x=x_0+k imes lcm(m_1,m_2)$

    满足这个条件的x就满足第一第二两个方程。满足第一第二两个方程的所有的解也都是这个方程的解。

    所以第一第二个方程和这个方程是等价的。

    将这个方程转化一下,可以得到新的同余方程:$x=x_0 (mod lcm(m_1,m_2))$

    这样,我们成功的把两个方程转化成了一个方程,以此类推。

    最后留下的这个方程,它的x_0的最小非负数解,就是我们要的最终答案!!!!!!!!!!!!!!

    例题:poj2891

    这是一个模板题,直接上代码:

    #include<cstdio>
    #include<algorithm>
    #include<iostream>
    using namespace std;
    typedef long long ll;
    const int N=100000+10;
    int n;
    ll a[N],r[N];
    ll exgcd(ll a,ll b,ll &x,ll &y){
        if(b==0){
            x=1,y=0;return a;
        }    
        ll ret=exgcd(b,a%b,y,x);
        y-=(a/b)*x;
        return ret;
    }
    ll excrt(){
        ll M=a[1],R=r[1],x,y,d;
        for(int i=2;i<=n;i++){
            d=exgcd(M,a[i],x,y);
            if((R-r[i])%d) return -1;
            x=(R-r[i])/d*x%a[i];
            R-=M*x;
            M=M/d*a[i];
            R%=M;
        }
        return (R%M+M)%M;
    }
    int main()
    {
        while(scanf("%d",&n)!=EOF){
            for(int i=1;i<=n;i++){
                scanf("%lld%lld",&a[i],&r[i]);
            }
            printf("%lld
    ",excrt());
        }
        return 0;
    }
    poj 2891

    例题:NOI2018屠龙勇士

    1.exgcd转化成同余方程的一般形式

    2.求解同余方程

    细节:

    1.Pi=1的情况特殊处理

    2.开始的ai、ki可以不用mod p[i]以免出现0的情况较为麻烦

    代码:

    #include<bits/stdc++.h>
    #define reg register int
    #define il inline
    #define fi first
    #define se second
    #define mk(a,b) make_pair(a,b)
    #define numb (ch^'0')
    #define int long long
    using namespace std;
    typedef long long ll;
    template<class T>il void rd(T &x){
        char ch;x=0;bool fl=false;
        while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true);
        for(x=numb;isdigit(ch=getchar());x=x*10+numb);
        (fl==true)&&(x=-x);
    }
    template<class T>il void output(T x){if(x/10)output(x/10);putchar(x%10+'0');}
    template<class T>il void ot(T x){if(x<0) putchar('-'),x=-x;output(x);putchar(' ');}
    template<class T>il void prt(T a[],int st,int nd){for(reg i=st;i<=nd;++i) ot(a[i]);putchar('
    ');}
    
    namespace Miracle{
    const int N=1e5+5;
    int n,m;
    ll a[N],p[N],st[N],k[N];
    multiset<int>s;
    multiset<int>::iterator it;
    ll ad(ll x,ll y,ll mod){
        return x+y>=mod?x+y-mod:x+y;
    }
    ll qk(ll x,ll y,ll mod){
        ll ret=0;
        while(y){
            if(y&1) ret=ad(ret,x,mod);
            x=ad(x,x,mod);
            y>>=1;
        }
        return ret;
    }
    void clear(){
        s.clear();
    }
    ll exgcd(ll a,ll b,ll &x,ll &y){
        if(!b){
            x=1;y=0;return a;
        }
        ll ret=exgcd(b,a%b,y,x);
        y-=(a/b)*x;
        return ret;
    }
    ll wrk(){//warning!!!  x!=0
    //    cout<<" kk "<<endl;
    //    prt(k,1,n);
        for(reg i=1;i<=n;++i){//warning!!!
    //        a[i]%=p[i];
    //        k[i]%=p[i];
    //        if(!k[i]) return -1;
            ll x=0,y=0;
            ll g=exgcd(k[i],p[i],x,y);
            x=(x%(p[i]/g)+(p[i]/g))%(p[i]/g);
            if(a[i]%g) return -1;
            x=qk(x,a[i]/g,p[i]/g);
            p[i]=p[i]/g;
            a[i]=x;
        }
    //    cout<<" after 1 "<<endl;
    //    prt(a,1,n);
    //    prt(p,1,n);
        for(reg i=1;i<n;++i){
    //        cout<<" turn "<<i<<endl;
            ll p1=p[i],p2=p[i+1];
    //        cout<<" p1 "<<p1<<" p2 "<<p2<<" a1 "<<a[i]<<" a2 "<<a[i+1]<<endl;
            ll x=0,y=0;
            ll g=exgcd(p1,p2,x,y);
            ll tmp=a[i+1]-a[i];
            ll mo=p1/g*p2;//lcm
            tmp=(tmp%mo+mo)%mo;
    //        cout<<" mo "<<mo<<" tmp "<<tmp<<" x "<<x<<" y "<<y<<endl;
            if(tmp%g) return -1;
            x=(x%(mo/p1)+(mo/p1))%(mo/p1);
            x=qk(x,tmp/g,mo/p1);
            ll nw=(a[i]+qk(x,p1,mo))%mo;
            a[i+1]=nw;
            p[i+1]=mo;
        }
        if(!a[n]) a[n]+=p[n];
        return a[n];
    }//warning!!! x!=0
    int main(){
        int t;
        rd(t);
        while(t--){
            clear();
            rd(n);rd(m);
            bool fl=true;
            for(reg i=1;i<=n;++i) rd(a[i]);
            for(reg i=1;i<=n;++i) {
                rd(p[i]);
                if(p[i]!=1) fl=false;
            }
            for(reg i=1;i<=n;++i) rd(st[i]);
            ll x;
            for(reg i=1;i<=m;++i){
                rd(x);s.insert(x);
            }
            for(reg i=1;i<=n;++i){
                it=s.upper_bound(a[i]);
                if(it!=s.begin()){
                    --it;
                    k[i]=*it;
                    s.erase(it);
                }else{
                    it=s.begin();
                    k[i]=*it;
                    s.erase(it);
                }
                s.insert(st[i]);
            }
            if(fl){
                ll ans=0;
                for(reg i=1;i<=n;++i){
                    ans=max(ans,(a[i]+k[i]-1)/k[i]);
                }
                printf("%lld
    ",ans);
            }else{
                printf("%lld
    ",wrk());
            }
        }
        return 0;
    }
    
    }
    signed main(){
        Miracle::main();
        return 0;
    }
    
    /*
       Author: *Miracle*
       Date: 2019/4/2 19:09:28
    */
    View Code
  • 相关阅读:
    62. Unique Paths
    24. Swap Nodes in Pairs
    83. Remove Duplicates from Sorted List
    21. Merge Two Sorted Lists
    141. Linked List Cycle
    268. Missing Number
    191. Number of 1 Bits
    231. Power of Two
    9. Palindrome Number
    88. Merge Sorted Array
  • 原文地址:https://www.cnblogs.com/Miracevin/p/9254795.html
Copyright © 2020-2023  润新知