zoukankan      html  css  js  c++  java
  • 扩展欧几里得&&中国剩余定理(学习笔记)

    原本是想把CRT、扩展CRT、欧几里得、扩展欧几里得都写在这,但由于博主太菜,刚刚才会EXCRT qwq

    在退组边缘徘徊的我果然还是菜了一点啊!!!

    布吉岛为什么但就是想奶一口gql省队稳了2333

    不闲扯了,进入正题!


    欧几里得(gcd)&&扩展欧几里得(exgcd):

    先来一个众人皆知的欧几里得算法(gcd ( a,b )=gcd ( b,a mod b ))

    证明过程自己完成,这里不多加叙述

    再来普及一个简单的裴蜀定理:若a,b是整数,且 (gcd(a,b)=d),那么一定存在整数 (x,y),使得(ax+by=d) 成立。

    裴蜀定理证明请看这里(主要是因为我之前打了一段后发现没证完然后就懒得打了qwq)

    那么 (fHuge ext{扩展欧几里得}​) 来了!!!

    扩展欧几里得算法:可以用来求解 (ax+by=gcd(a,b)) ,也就是求同余方程 (axequiv gcd(a,b)(mod b)) 的解

    算法过程及解如下:

    我们要求解 (ax+by=gcd(a,b)) 其实主要是想算出 (x) 的值,(y) 是一个辅助解

    然后我们构造这么一个式子:(bx_1+(a mod b)y_1=gcd(b,(a mod b))) ,根据欧几里得算法,我们可得出 (bx_1+(a mod b)y_1=ax+by) ,不妨把 (a mod b) 变成 (a-b*(a/b)) ,然后带入原式:


    (bx_1+(a-b*(a/b))y_1=ax+by)

    (bx_1+ay_1-b*(a/b)y_1=ax+by)

    (ay_1+b(x_1-(a/b)y1)=ax+by)


    那么我们求出了一组解:(x=y_1,y=x_1-(a/b)y_2)

    然后将我们求解的两个式子对比一下:


    (ax+by=gcd(a,b))

    (bx_1+(a mod b)y_1=gcd(b,(a mod b)))


    是不是发现了什么?是不是有点像欧几里得算法?

    然后我们按这种方式递推下去,直到 (b=0) ,那么:


    (ax+by=gcd(a,b))

    (ax=gcd(a,0))


    显然 (gcd(a,0)=a,x=1) ,此时已经和 (y) 的值没有关系了,即此时 (y) 的值是任意数,但是这里建议把 (y) 赋成 (0) ,可以避免在返回中爆 (long long)

    之后我们已经求出来了关于同余方程的一个解 (x) ,虽然 (x) 不一定是最小的,但显然 (x) 加上或减去 (b) 是没有任何影响的,所以用一个 (x = (x \% b + b) \% b) 就可以了求出满足同余方程的最小正整数解了

    上代码:

    void exgcd(int a,int b,int &x,int &y)
    {//x,y都是要返回的值
        if(b==0)
        {
            x=1,y=0;
            return;
        }
        exgcd(b,(a%b),x,y);//递归
        int xx=y;//记录一下解
        y=x-(a/b)*xx,x=xx;
        return;
    }
    

    至此,扩展欧几里得讲解完毕qwq


    中国剩余定理(CRT):

    中国剩余定理其实真的不难qwq

    一般情况下我们是要求解如下式子:

    (egin{cases} xequiv a_1(mod b_1)quad \ xequiv a_2(mod b_2)quad \ ...quad \ xequiv a_n(mod b_n)quad \ end{cases})

    其中所有的 (b_1,b_2...b_n​) 都互质,这里的 (a_i​) 都小于 (b_i​)

    首先我们构造一个数 (B=b_1*b_2*...*b_n​) ,那么显然当我们求出 (x​) 之后加上或减去 (B​) 都是成立的。

    接下来考虑对于每一个同余方程的处理:

    对于同余式 (xequiv a_k(mod b_k)) ,我们可以构造一个 (B_k=B/b_k) ,显然 (gcd(B_k,b_k)=1) ,根据裴蜀定理可得:存在整数 (i,j) 使得 (iB_k+jb_k=1) ,即 (iB_kequiv 1(mod b_k)) 。因为 (iB_k mod b_k=1) ,那么有 (a_i*iB_k mod b_k=a_i)(i) 可以通过扩展欧几里得求得,所以该同余方程的一个解是 (x=a_i*i*B_k)

    再回到整个方程组,我们继续构造一个数 (x=a_1i_1B_1+a_2i_2B_2+...+a_ni_nB_n) ,那么这个数就是同余方程组的一个解。因为对于第 (k) 个同余方程,(B_{1...n}) 中除 (B_k) 之外所有数都是 (b_k) 的倍数,所以同余方程是成立的。

    放一下CRT的代码:

    void exgcd(int a,int b,int &x,int &y){
        if(b==0)
        {
            x=1,y=0;
            return;
        }
        exgcd(b,(a%b),x,y);
        int xx=y;
        y=x-(a/b)*xx,x=xx;
        return;
    }//扩展欧几里得求解同余方程
    
    lt china()
    {
        for(int i=1;i<=n;i++) N*=B[i];
        for(int i=1;i<=n;i++)
        {
            int a=N/B[i],b=B[i],x=0,y=0;
            exgcd(a,b,x,y);//此处的x还要变成最小正整数解
            ans+=a*((x%b+b)%b)*A[i];
        }
        return (ans%N+N)%N;
    }
    

    扩展中国剩余定理(EXCRT)

    好了,同样是形如下面的式子:

    (egin{cases} xequiv a_1(mod b_1)quad \ xequiv a_2(mod b_2)quad \ ...quad \ xequiv a_n(mod b_n)quad \ end{cases}​)

    只不过这次不保证所有的 (b_1,b_2,...,b_n) 互质

    因为不互质,所以我们无法像CRT一样使用扩欧。

    考虑你已经求出来了前 (k-1​) 个同余方程的解,得到的解为 (ans​) ,设 (lcm​) 为前 (k-1​) 个方程中所有的 (b_i​) 的最小公倍数,则前 (k-1​) 个方程的解为 (x=ans+i*lcm​) ,而我们只需要确定一个 (i​) 使得 (ans+i*lcmequiv a_k(mod b_k)​) ,然后更新一下 (lcm​) 就好了。

    又是一波转化:


    (ans+i*lcmequiv a_k(mod b_k))

    (i*lcmequiv a_k-ans(mod b_k))


    注意:此处的 (lcm)(b_k) 不一定互质,需要稍加处理,设 (gcd=gcd(lcm,b_k),c=(a_k-ans) mod b_k) ,由上式可得:

    (i*lcm+h*b_k=c)


    由裴蜀定理可得此方程有解的必要条件是 (gcd|c) ,于是继续变形:


    (i/gcd*lcm+h*b_k/gcd=c/gcd)

    (i/gcd*lcmequiv c/gcd(mod b_k/gcd))


    此时 (gcd(lcm,b_k/gcd)=1) ,对于上面这个方程,我们依然可以先求出一个 (j) 满足 (j*lcmequiv 1(mod b_k/gcd)) ,然后再乘以 (c/gcd) 倍就好了。

    具体代码中有解释:

    int exgcd(int a,int b,int &x,int &y)
    {//扩展欧几里得,一并求出gcd
        if(b==0)
        {
            x=1,y=0;
            return a;
        }
        int res=exgcd(b,(a%b),x,y);
        int xx=x;
        x=y,y=xx-(a/b)*y;
        return res;
    }
    
    int mul(int s,int p,int mod)
    {//这是一个快速乘,防止s*p爆long long
        int res=0;
        while(p)
        {
            if(p%2) res=(res+s)%mod;
            s=(s+s)%mod;
            p/=2;
        }
        return res%mod;
    }
    
    int EXCRT()
    {
        m=b[1],ans=a[1];//第一个方程要特殊处理,直接赋值就好了
        for(int i=2;i<=n;i++)
        {
            int x,y,c=(a[i]-ans%b[i]+b[i])%b[i];
            //c就是前面的ak-ans
            int gcd=exgcd(m,b[i],x,y);
            x=mul(x,c/gcd,b[i]/gcd);
            //这一段我也弄了好久,最后终于搞懂了
            //大致把x为什么要乘c/gcd,和%(b/gcd)的原因写在上面了
            //不懂欢迎提问
            ans+=x*m;//更新ans
            m*=(b[i]/gcd);//更新lcm
            ans=(ans%m+m)%m;
        }
        return (ans%m+m)%m);
    }
    

    终于更完了qwq心累

  • 相关阅读:
    android高级页面效果集锦
    2018年Android的保活方案效果统计
    程序员如何预估自己的项目开发时间?
    Google开发者大会:你不得不知的Tensorflow小技巧
    练就Java24章真经—你所不知道的工厂方法
    一个完整Java Web项目背后的密码
    怎么捕获和记录SQL Server中发生的死锁
    使用跟踪标志位分析死锁
    通过SQL Server Profiler来监视分析死锁
    SQL Server中的死锁
  • 原文地址:https://www.cnblogs.com/ajy-shi-cj-zui-cai/p/10420706.html
Copyright © 2011-2022 走看看