zoukankan      html  css  js  c++  java
  • 扩展欧几里得(exgcd)与同余详解

    exgcd入门以及同余基础

    gcd,欧几里得的智慧结晶,信息竞赛的重要算法,数论的...(编不下去了

    讲exgcd之前,我们先普及一下同余的性质:

    1. left ( a+b 
ight )mod p=left (amod p+ bmod p
ight )mod p
    2. left ( a-b 
ight )mod p=left (amod p- bmod p
ight )mod p
    3. left ( a	imes b 
ight )mod p=left (amod p)	imes (bmod p
ight )mod p
    4. aequiv b(mod p),那么a^{n}equiv b^{n}(mod p)
    5. amod p1=xamod p2=x,且p1,p2互质,amod (p	imes q)=x

     有了这三个式子,就不用怕在计算时溢出了。

    下面我会用gcdleft ( a,b 
ight )lcmleft ( a,b 
ight )分别表示a与b的最大公约数与最小公倍数。

    首先会来学扩欧的同学肯定都会欧几里得算法(即辗转相除法)了吧

    而通过观察发现:lcmleft ( a,b 
ight )=frac{a	imes b}{gcdleft ( a,b 
ight )}=frac{a}{gcdleft ( a,b 
ight )}	imes b,先除后乘防溢出。

    所以gcdleft ( a,b 
ight )lcmleft ( a,b 
ight )的代码如下:

    1 inline int gcd(int a,int b)
    2 {return (b==0)?a:gcd(b,a%b);}
    3 inline int lcm(int a,int b)
    4 {return a/gcd(a,b)*b;

    讲exgcd之前先引入一种方程——不定方程

     所谓不定方程,是指未知数的个数多于方程个数,且未知数受到某些限制(如要求是有理数、整数或正整数等等)的方程或方程组。

                                                                                             ————百度百科 

    就是形如ax+by=c的方程,其中a,b,c已知。

    1.判断是否有解

    如果cmod gcdleft ( a,b 
ight )
eq 0,那么方程无解。

    2.转化

    方程可转化为{a}'x+{b}'y={c}'

    其中{a}'=frac{a}{gcdleft (a,b 
ight )}{b}'=frac{b}{gcdleft (a,b 
ight )}{c}'=frac{c}{gcdleft (a,b 
ight )}

    3.求一组特解

    接着就用到了exgcd。

    我们知道gcd有一个性质gcd(a,b) = gcd(b, amod b)

    如果,一直循环下去,b将等于0,那么x将等于c/a,y=0。

    1 inline void exgcd(int a,int b,int c)
    2 {
    3     if(!b)
    4     {x=c/a;y=0;return;}
    5     exgcd(b,a%b,c);
    6     x=y;
    7     y=(c-a*x)/b;
    8     return;
    9 }

    这就求出了一组特解。

    exgcd的模板我也在这摆出来

    1 inline void exgcd(int a,int b)
    2 {
    3     if(!b)
    4     {x=1;y=0;return;}
    5     exgcd(b,a%b);
    6     k=x;x=y;
    7     y=k-a/b*y;
    8     return;
    9 }

    这是求ax+by 
ight =gcdleft ( a,b 
ight )时用的,后面讲同余方程会讲。

    4.构造通解

    我们假设x1,y1是我们求出的一组特解,那么

    left{egin{matrix} & &x=ccdot x1+ccdot bcdot k \ & & y=ccdot y1-ccdot acdot k end{matrix}
ight.     left ( kepsilon z 
ight )


    同余类问题

    1.单个同余方程

    求的是axequiv b( mod p)关于x的解

    转化一下,就成了ax+py=b,就可以直接套exgcd模板。

    2.同余方程组

    left{egin{matrix} & &xequiv b1left (mod p1 
ight )\ & & xequiv b2left (mod p2 
ight ) end{matrix}
ight.

    1.有解的充要条件left ( b1-b2 
ight )mod gcd(p1,p2)=0

    2.left{egin{matrix} & &x-p1	imes y1=b1\ & & x-p2	imes y2=b2end{matrix}
ight.

    下式减上式得p1	imes y1-p2	imes y2=b2-b1

    再用exgcd求出y1和y2,{x}'=b1+p1	imes y1=b2+p2	imes y2

    3.关于通解

    所有的x mod lcm(p1,p2)有唯一解,这样就可以通过特解,求通解了。

    4.至于式子更多的同余方程组,就先联立两个,就可以得出新的方程xequiv {x}'left ( mod lcm(p1,p2) 
ight )

    再联立下一个。


    exgcd用处及题目讲解

    1.求同余方程的解

    例如这道题P1082

    这是一道裸的扩欧模板题,变形之后就是求ax+by=1

    套模板即可。

     1 inline void exgcd(int a,int b)
     2 {
     3     if(b==0)
     4     {x=1;y=0;return;}
     5     exgcd(b,a%b);
     6     k=x;x=y;
     7     y=k-a/b*y;
     8     return;
     9 }
    10 int main()
    11 {
    12     int n,m;
    13     read(n),read(m);
    14     exgcd(n,m);
    15     printf("%d",(x+m)%m);
    16 }

    还有一道模板P1516

    仔细观察,推一下后我们发现,这在就是在求x+k	imes mequiv y+k	imes n(mod l)的解。

    进而可以推出x+k	imes m-y-k	imes n=l	imes p

    合并同类项后(x-y)+k	imes(m-n)-l	imes p=0

    把一些东西移到左边来后(x-y)=k	imes(n-m)+l	imes p

    把(x-y),(n-m)各看成一个整体后,问题就成了解一个不定方程。

     1 inline int exgcd(long long a,long long b)
     2 {
     3     if(b==0)
     4     {x=1;y=0;return a;}
     5     ans=exgcd(b,a%b);
     6     k=x;x=y;
     7     y=k-a/b*y;
     8     return ans;
     9 }
    10 int main()
    11 {
    12     long long x1,y1,m,n,l;
    13     read(x1),read(y1),read(m),read(n),read(l);
    14     if(n-m<0)swap(x1,y1);
    15     exgcd(std::abs(n-m),l);
    16     if((x1-y1)%ans!=0)
    17     printf("Impossible");
    18     else
    19     printf("%lld",((x*((x1-y1)/ans))%(l/ans)+(l/ans))%(l/ans));
    20 }

    还有一道也是模板P4777,涉及同余方程组求解,上面已详细的讲了,近期我也会发一篇中国剩余定理的博客

     1 inline long long mul(long long a,long long b,long long mod)
     2 {
     3     long long res=0;
     4     while(b>0)
     5     {
     6         if(b&1) res=(res+a)%mod;
     7         a=(a+a)%mod;
     8         b>>=1;
     9     }
    10     return res;
    11 }
    12 long long exgcd(long long a,long long b,long long &x,long long &y)
    13 {
    14     if(!b)
    15     {x=1;y=0;return a;}
    16     long long gcd=exgcd(b,a%b,x,y);
    17     k1=x;x=y;
    18     y=k1-a/b*y;
    19     return gcd;
    20 }
    21 int main()
    22 {
    23     io::begin();
    24     io::read(n);
    25     for(register int i=1;i<=n;i++)
    26     io::read(b1[i]),io::read(a1[i]);
    27     long long x,y,k;
    28     long long m=b1[1],ans=a1[1];
    29     for(int i=2;i<=n;i++)
    30     {
    31         long long a=m,b=b1[i],c=(a1[i]-ans%b+b)%b;
    32         long long gcd=exgcd(a,b,x,y);
    33         long long p=b/gcd;
    34         x=mul(x,c/gcd,p);
    35         ans+=x*m;
    36         m*=p;
    37         ans=(ans%m+m)%m;
    38     }
    39     printf("%lld",(ans%m+m)%m);
    40 }

    2.扩欧求逆元

    这是一种很重要的算法,至于逆元怎么跟扩欧扯上关系,大家可以点这里乘法逆元及两道模板题详解

    这里就不多赘述了,大家可以用扩欧a一下P3811,P2613。


    我要讲的讲完了,如果觉得讲的还好,请关注我的blog,谢谢

  • 相关阅读:
    Atitit 趋势管理之道 attilax著
    Atitit 循环处理的新特性 for...else...
    Atitit 2017年的技术趋势与未来的大技术趋势
    atitit 用什么样的维度看问题.docx 如何了解 看待xxx
    atitit prj mnrs 项目中的几种经理角色.docx
    Atitit IT办公场所以及度假村以及网点以及租房点建设之道 attilax总结
    Atitit 工具选型的因素与方法 attilax总结
    Atitit.团队文化建设影响组织的的一些原理 法则 定理 效应 p826.v4
    Atiitt 管理方面的误区总结 attilax总结
    Atitit 未来趋势把控的书籍 attilax总结 v3
  • 原文地址:https://www.cnblogs.com/mashiro-/p/9526323.html
Copyright © 2011-2022 走看看