acm数论之欧几里得gcd-白红宇

acm数论之欧几里得gcd

阅读量：6324 次

发布时间：2019-06-22

本文共 10576 字，大约阅读时间需要 35 分钟。

1.欧几里得定理

同余定理的公式：(a+b)%mod=(a%mod+b%mod)%mod

(a*b)%mod=(a%mod*b%mod)%mod

扩展欧几里得也有自己的一个公式：a*x+b*y=gcd(a,b)

1.1定义：

gcd(a,b)=gcd(b,a%b)

1.2用法：

获得最大公约数与最小公倍数----gcd(a,b)lcm(a,b)=ab;

int gcd(int a,int b)   //最大公约数{    if(b==0)  return a;    else return gcd(b,a%b);}     int lcm(int a,int b)  //最小公倍数{    return a/gcd(a,b)*b;    //防止溢出}

2.扩展欧几里得

2.1定义：

ax+by=gcd(a,b)=d，在已知a,b的情况下求解出一组x,y

2.2推导过程：

因为a%b=a-(a/b)b, 则有 d=bx1+[a-(a/b)b]y1=bx1+ay1-(a/b)by1=ay1+b(x1-a/by1), 故 x=y1,y=x1-a/by1

2.3性质：

1.若通过扩展欧几里得求出一组特解(x0,y0)，那么有ax0+by0=d.

则方程的通解为，其中k为任意整数 x0=d/exgcd()*x，y0=d/exgcd()*y；

x=x0+k*(b/d)--------x = x0 +( b / gcd( a, b ) ) * k

y=y0-k*(a/d)--------y = y0 - ( a / gcd( a, b ) ) * k

2.已知ax+by=d的解，对于ax+by=c的解，c为任意正整数，只有当d|c时才有解

其通解为

x=(c/d)x0+k(b/d)

y=(c/d)y0-k(a/d)

//扩展欧几里得处理ax+by=gcd(a,b)=d解的问题 int exgcd(int a,int b,int &x,int &y)  //得到gcd(a,b)的值 {    if(b==0)    {        x=1,y=0;        return a;     }      int r=exgcd(b,a%b,x,y);     int t=x; x=y; y=t-a/b*y;     return r;//gcd(a,b)}

2.4常见用法：

（1）求形如ax+by=c的通解，或从中选取某些特解

（2）求乘法逆元

（3）求解线性同余方程

写的超级好的博客：

解题思路：

1.先化成不定方程ax+by=c

2.a,b参数取正数，判断a,b是否为负数，为负数就a,b,c都乘上-1，另外注意若a,b中只有一个是负数，那另一个不需要乘-1，因为a,b本身是含未知数参数的。

3.r=exgcd(a,b,x,y),然后判断是否有解采用if(c%r==0)

求特解：x0=x*c/r, y0=y*c/r;

求通解：x=x0+b/r*c , y=y0-a/r*c （其中 t 为整数）

求最小解：int s=b/r; minx= (x0%s+s)%s;//最小解

2.4.1例题实例：

问题描述：对于 ax+by=c 的不定方程求通解或特解？

设 r=gcd(a,b)，

若 c%r!=0 这此方程无整数解

若 c%r==0，特解为x1=x0*c/r , y1=y0*c/r，通解 x=x1+b/r*t , y=y1-a/r*t （其中 t 为整数）

不定方程的最小解：

r=exgcd(a,b,x,y);

x=x*c/r;

int s=b/r;

minx= (x%s+s)%s;//最小解

例题1：

题目链接：

题目大意：0<a, b<=2^31，X>=0，给定a,b求满足X*a + Y*b = 1的最小正数X，以及与之对应的Y。

解题思路：求形如ax+by=c的最小解，x=x0+k*(b/d)， y=y0-k*(a/d)，从通解可看出，其中b / d 都正数，x 和 t成正比，当x0为正数的时候 t == 0，有最小值；当x0为负数的时候，x最小就是0，与横轴相交，所以可以求出 x0 + b / t * t = 0时候的 t，因为此时 t是整数，不一定是准确相交的那一个点，所以当前的x可能是负数，如果为负数，t++，取下一个整数点。

ac代码：

1 #include
      
        2 using namespace std; 3 //扩展欧几里得处理ax+by=gcd(a,b)=d解的问题  4 int exgcd(int a,int b,int &x,int &y)  //得到gcd(a,b)的值  5 { 6     if(b==0) 7     { 8         x=1,y=0; 9         return a;10      } 11      int r=exgcd(b,a%b,x,y);12      int t=x; x=y; y=t-a/b*y;13      return r;//gcd(a,b)14 }15 int main()16 {17     int a,b;18     while(~scanf("%d%d",&a,&b))19     {20         int x,y;21         int d=exgcd(a,b,x,y);22         if(1%d)//有余数 23         {24             cout<<"sorry"<
       
        0)//由于a,b为正数，所以k=0取最小 32             cout<
        
         <<" "<
         
          <

View Code

例题2：

题目链接：

题意：一个人要从A走到B 只能走a布、b步、(a+b)步，可以往左或右走，问最少要走几步，不能走到的输出-1

题目思路：可以设走x步a米的，走y部b米的，得到式子：ax + by = B - A;

通过扩展欧几里德可得x和y

x、y的解集：X = x + b/gcd(a,b) * t; Y = y - a /gcd(a,b) * t;

x与y最接近时为答案，这样的话c可以取的更多，那么步数会最小

它们

（1）X、Y同号时，即方向相同，可以走a+b来减少步数，取max(X,Y)（这里可以画X、Y的坐标图来帮助理解）

（2）X、Y异号时，即方向相反，步数绝对值相加，abs(X) + abs(Y)

1 #include
      
        2 #include
       
         3 #include
        
          4 using namespace std; 5 const long long maxn = 0x3f3f3f3f; 6 //#define maxn 0x3f3f3f3f 7 long long exgcd(long long a,long long b,long long &x,long long &y) 8 { 9     if(b==0)10     {11         x=1,y=0;12         return a;13     }14     long long r=exgcd(b,a%b,x,y);15     long long t=x;x=y;y=t-a/b*y;16     return r;17 }18 long long Abs(long long a)19 {20     if(a<0)21     a=-1*a;22     return a;23 }24 int main()25 {26     long long a,b,A,B;27     long long ans;28     int T;29     cin>>T;30     while(T--)31     {32         cin>>A>>B>>a>>b;33         long long x,y;34         long long h=exgcd(a,b,x,y);35         if((B-A)%h!=0)36             printf("-1\n");37         else{38             x=(B-A)/h*x;//注释修改39             y=(B-A)/h*y;//注释修改40             long long minx,miny;41             long long k;42             b=b/h;43             a=a/h;44             k=(y-x)/(b+a);//注释修改45             ans=maxn*maxn;   //注释修改  这里定义成maxn数据定义小了46             for(int i=k-1;i<=k+1;i++)47             {48                 minx=x+b*i,miny=y-a*i;49                 long long m;50                 if(Abs(minx+miny)==Abs(minx)+Abs(miny))//说明x与y同号51                     m=max(Abs(minx),Abs(miny));52                 else53                     m=Abs(minx)+Abs(miny);54                 if(m
         
          0&&miny>0)||(minx<0&&miny<0))//x与y同号57 //                {58 //                m=Abs(min(minx,miny))+Abs(minx-miny);//这句有问题的  测试这个存在数据Abs(min(minx,miny))+Abs(minx-miny)！=max(Abs(minx),Abs(miny));59 //                }60 //                 else61 //                 m=Abs(minx)+Abs(miny);62 //                 if(m

ac代码

wa的心得或者可以说是经验：

1.maxn定义不同而wa:

const long long maxn = 0x3f3f3f3f; //正确

//#define maxn 0x3f3f3f3f 错误

2.除法运算要考虑是否倍数除，可以倍数除最好倍数除

x=(B-A)/h*x; y=(B-A)/h*y;//正确

x=x/h*(B-A); y=y/h*(B-A); //错误

3.ans=maxn*maxn; //注释修改这里定义成maxn数据定义小了

3.前一个为正确写法，后一个为错误写法

分析：两者的不同在于在for循环里面minx,miny，后一种不能保证i=k-1~k+1这里面取得答案，因为多了一个参数。

1 b=b/h;2 a=a/h;3 k=(y-x)/(b+a);//注释修改4 for(int i=k-1;i<=k+1;i++)5 minx=x+b*i,miny=y-a*i;

1 b=b;2 a=a;3 k=(y-x)/(b+a)*h;4 for(int i=k-1;i<=k+1;i++)5 minx=x+b*i/h,miny=y-a*i/h;

2.4.2例题实例：求乘法逆元

问题描述：求解一个数 a 对于另一个数 m 的乘法逆元：a*x=1(mod m) 也即：a*x%m=1求x 也即a*x+m*y=1

当gcd(a , m) != 1 ,a*x + m*y = 1没有解

当1 % gcd(a , m) == 0，a*x + m*y = 1有解

x0=x*1/exgcd();

最小的解:x0 % m

x 的通解: x0 + m*t

#include
     
      #include
      
       using namespace std;//扩展欧几里得处理ax+by=gcd(a,b)=d解的问题 int exgcd(int a,int b,int &x,int &y)  //得到gcd(a,b)的值 {    if(b==0)    {        x=1,y=0;        return a;     }      int r=exgcd(b,a%b,x,y);     int t=x; x=y; y=t-a/b*y;     return r;//gcd(a,b)}//求解一个数 a 对于另一个数 m 的乘法逆元：a*x=1(mod m) 也即：a*x%m=1求x   也即a*x+m*y=1 int inversenum(int a,int m){    int x,y,ans,gcd;    gcd=exgcd(a,m,x,y);    if(1%gcd!=0)//无解     return -1;    x=x*1/gcd;    m=abs(m);    ans=x%m;    if(ans<=0)    ans+=m;    return ans;}int main(){    int a,m;    cin>>a>>m;    cout<
       
        <

例题：ZOJ3609 Modular Inverse

1 #include
      
        2 #include
       
         3 using namespace std; 4 //扩展欧几里得处理ax+by=gcd(a,b)=d解的问题  5 int exgcd(int a,int b,int &x,int &y)  //得到gcd(a,b)的值  6 { 7     if(b==0) 8     { 9         x=1,y=0;10         return a;11      } 12      int r=exgcd(b,a%b,x,y);13      int t=x; x=y; y=t-a/b*y;14      return r;//gcd(a,b)15 }16 //求解一个数 a 对于另一个数 m 的乘法逆元：a*x=1(mod m) 也即：a*x%m=1求x   也即a*x+m*y=1 17 int inversenum(int a,int m)18 {19     int x,y,ans,gcd;20     gcd=exgcd(a,m,x,y);21     if(1%gcd!=0)//无解 22     return -1;23     x=x*1/gcd;24     m=abs(m);25     ans=x%m;26     if(ans<=0)27     ans+=m;28     return ans;29 }30 int main()31 {32     int T;33     cin>>T;34     while(T--){35     int a,m;36     cin>>a>>m;37     if(inversenum(a,m)==-1)38     cout<<"Not Exist"<

View Code

2.4.2例题实例：求线性同余方程

问题描述：关于 x 的模方程 ax%b=c (或者是(x+mt)%L=(y+nt)%L)的解，方程转换为 ax+by=c 其中 y 一般为非正整数

解得 x1=x0*c/r，通解为 x=x1+b/r*t

设 s=b/r （已证明 b/r 为通解的最小间隔），则 x 的最小正整数解为 (x1%s+s)%s

例题：poj 1061

根据题意：可写出(x+mt)%L=(y+nt)%L，求可取的最小t

(m-n)*t=(y-x)%L; (m-n)*t+L*b=y-x; (m-n)*a+L*b=y-x; 看成a*m1+b*n1=t;(m1=m-n,n1=L,t=y-x)

注意：数据类型应该为long long

1 #include
      
        2 #include
       
         3 using namespace std; 4 //扩展欧几里得处理ax+by=gcd(a,b)=d解的问题  5 long long t; 6 long long exgcd(long long a,long long b,long long &x,long long &y)  //得到gcd(a,b)的值  7 { 8     if(b==0) 9     {10         x=1,y=0;11         return a;12      } 13      long long r=exgcd(b,a%b,x,y);14      long long t=x; x=y; y=t-a/b*y;15      return r;//gcd(a,b)16 }17 long long solve(long long m,long long n)18 {19     long long x,y;20     long long r=exgcd(m,n,x,y);21     if(t%r)22     {23         return -1;24     }25     x=x*t/r;26     long long s=n/r;27     return (x%s+s)%s;//不定方程的最小解x 28 }29 int main()30 {31     long long x,y,m,n,L;32     cin>>x>>y>>m>>n>>L;33     //不定式a*(m-n)+b*L=y-x;   看成a*m+b*n=t; 34     t=y-x;35     m=m-n;36     n=L;37     if(m<0)//整个不定式a,b取正数，由于b>0，b不需要变(b前面有系数),所以t需要改变 38     {39         m=-1*m;40         t=-1*t;41     }42     if(solve(m,n)==-1)43     cout<<"Impossible"<

View Code

例题：p2421荒岛野人

题目链接：

题目思路：暴力+exgcd

首先判断两个野人是否会在同一个洞穴：(Ci+Pi*t)%M==(Cj+Pj*t)%M (p[i]-p[j])*t=(C[j]-C[i])%M (P[i]-P[j])*a+M*b=C[j]-C[i]

暴力求M,从1递增，满足条件就停下来，check(M)

check(M)的作用是判断n个野人之间没有住在同一个洞穴，即余数不同或者是相同余数数字大于两野人最小寿命。

题目代码：

1 #include
      
        2 #include
       
          3 using namespace std; 4 int c[200],p[200],l[200],n; 5 int gcd(int a,int b){
    return a%b==0?b:gcd(b,a%b);} 6 //扩展欧几里得处理ax+by=gcd(a,b)=d解的问题  7 int exgcd(int a,int b,int &x,int &y)  //得到gcd(a,b)的值  8 { 9     if(b==0)10     {11         x=1,y=0;12         return a;13      } 14      int r=exgcd(b,a%b,x,y);15      int t=x; x=y; y=t-a/b*y;16      return r;//gcd(a,b)17 }18 bool go(int i,int j,int b){19       int a=p[i]-p[j],f=c[j]-c[i],g,x,y;20       if(a<0)21       {22           a=a*-1,f=-1*f;23       }24       g=exgcd(a,b,x,y);25       if(f%g==0){26           x=x*f/g;27           b=b/g;28           x=(x%b+b)%b;//最小解29           if(x<=min(l[i],l[j]))30           return 0;31       }32       return 1;33 }34 bool check(int m){35       int i,j,k;36       for(i=1;i

View Code

题目集应用：

第二题：hdu 1576

题目思路：求（A/B)%9973，已知条件有：n=A%9973,gcd(B,9973)=1,A%B=0;输入n,B，求（A/B)%9973 的值？

设（A/B)%9973=x,即（A/B)=9973h+x; A=9973*h*B+x*B ; n=A%9973即n=(9973*h*B+x*B)%9973，即n=(x*B)%9973，即x*B-9973y=n;求最小的正数x

1 #include
      
        2 using namespace std; 3 //扩展欧几里得处理ax+by=gcd(a,b)=d解的问题  4 int exgcd(int a,int b,int &x,int &y)  //得到gcd(a,b)的值  5 { 6     if(b==0) 7     { 8         x=1,y=0; 9         return a;10      } 11      int r=exgcd(b,a%b,x,y);12      int t=x; x=y; y=t-a/b*y;13      return r;//gcd(a,b)14 }15 int main()16 {17     int T,n,B;18     scanf("%d",&T);19     while(T--)20     {21         scanf("%d%d",&n,&B);22         int x,y;23         int d=exgcd(B,9973,x,y);24         x=x*n/d;//最小解 25         int s=9973/d;26         x=(x%s+s)%s;27         cout<
       
        <

View Code

第三题：UVA12169

题目大意：有3个整数 x[1], a, b 满足递推式x[i]=(a*x[i-1]+b)mod 10001。由这个递推式计算出了长度为2T的数列，现在要求输入x[1],x[3],......x[2T-1], 输出x[2],x[4]......x[2T]. T<=100,0<=x<=10000. 如果有多种可能的输出，任意输出一个结果即可。

题目思路：第一种直接暴力方法，枚举a,b，找到合适的a,b（即满足式子，又因为我们只知道x[1],x[3]..，所以只要它满足所有的x[i]=(a*((a*x[i-2]+b)mod 10001)+b)mod 10001）然后输出相对于应的x[2],x[4]...

ac代码：

1 #include
      
        2 using namespace std; 3 #define maxn 100005 4 int a[maxn]; 5 int main() 6 { 7     int n; 8     scanf("%d",&n); 9     for(int i=0;i

暴力法

第二种方法，比第一种方法快，我们枚举a,利用简化式子(a+1)*b mod 10001 = (x[3]-a*a*x[1]) mod 10001。这里就变成了同模方程，扩展欧几里得即可解答,求b。然后由a,b求相对于应的x[2],x[4]...

相应博客：