• 各种友(e)善(xin)数论总集,从入门到绝望6---同余系列!(ex)GCD,(ex)CRT


    参考文献

    SCY资料、手推

    博客:https://www.cnblogs.com/MashiroSky/p/5918158.html很好的阐述了中国剩余定理的公式等等等等,反正我看了这个博客很快就懂了。

    GCD的家族

    GCD

    如果我们要求两个数字的最大公约数怎么求?

    如果两个数字(a,b)存在最大公约数(k),那么把(a)写成(xk),把(b)写成(yk)

    那么我们知道(b\%a=(y\%x)k),但是(y)又与(x)互质,所以仅当(x=1)时,(y)为最大公约数,此时(a=k)也就是最大公约数,而且我们也知道(y\%x)也是与x互质的。

    所以我们可以设新的(a')与新的(b')(a'=b\%a),(b'=a),当(a'=0)时,(b)就是最大公约数。

    而且一开始(a>b)也不怕,做了一次交换后就换过来了((b\%a=b))

    int  gcd(int  x,int  y)
    {
    	if(x==0)return  y;
    	return  gcd(y%x,x);
    }
    

    扩展:更相减损法

    另外补充一种GCD的形式。

    二进制的GCD用的是更相减损法。

    首先,我们有两个数字(x,y)

    1. (x\%2==0,y\%2==0,gcd(x,y)=gcd(x>>1,y>>1)*2)
    2. (x\%2==0,y\%2==1,gcd(x,y)=gcd(x>>1,y))
    3. (x\%2==1,y\%2==0,gcd(x,y)=gcd(x,y>>1))
    4. (x\%2==1,y\%2==1,x>=y,gcd(x,y)=gcd((x-y)>>1,y)),其实这里的(x-y)不用除(2)也可以,只是这样写更快,反正都是偶数。
    inline  LL  gcd(LL  x,LL  y)
    {
    	int  ans=0;
    	while(x  &&  y)
    	{
    		if(x&1  &&  y&1)
    		{
    			y>x?x^=y^=x^=y:0;
    			x=(x-y)>>1;
    		}
    		else  if(x&1)y>>=1;
    		else  if(y&1)x>>=1;
    		else  x>>=1,y>>=1,ans++;
    	}
    	return  (x+y)<<ans;
    }
    

    常数好像比GCD更小,复杂度都是(log),但是貌似跑起来差不多,写写吧,反正都差不多了,还更稳一点,你说是吧。

    扩展欧几里得(EXGCD)

    如何求(ax+by=c)的一个整数解?

    尽看扩展欧几里得。

    首先,我们设(ax+by=gcd(a,b))

    那我们就想了,如果有个方程为(0x+gcd(a,b)y=gcd(a,b)),不就很妙了吗?
    (y=1,x)随便等于什么都不影响最终结果,只不过一般取0,怕爆long long。

    然后我们再证一个:

    当我们解出了((b\%a)x+ay=gcd(a,b)),怎么推到(ax+by=gcd(a,b))

    我们可以推导一下((⌊x⌋)为向下取整):

    [ax+by=(b-⌊b/a⌋*a)x'+ay' ]

    [ax+by=-a*⌊b/a⌋x'+ay'+bx' ]

    [ax+by=a(y'-⌊b/a⌋x')+bx' ]

    (⌊b/a⌋)就是C++的(b/a)

    所以,我们发现(b\%a)(a)就是(gcd)过程,所以最后也会得到(0)(gcd(a,b))

    然后再看:如果(gcd(a,b)|c)的话,设(d=c/gcd(a,b)),然后把(x,y)直接乘以(d)就没什么毛病了。

    但是如果不整除就无解,我们又可以证一证:


    首先(gcd(a,b)|c)是肯定有解的,我们可以知道,而且(ax+by=c)可以类似于(ax≡cmod b),我们可以根据同余的性质得出我们可以把(a,b,c)的同一个因数(gcd(a,b))提出,得到(a',b',c'),化成(a'x≡c'mod b'),而且a'与b'不互质,那么我们可以得出(欧拉定理证明过程可以说明),同余方程(ax≡cmod b)(a,b)互质时有解,或者说当(c)整除(gcd(a,b))时有解。(好像说偏了

    然后,如果(gcd(a,b))与c不整除,我们再设(z=gcd(a,b),c=kz+d)

    那么我们把原式化为:(a'zx+b'zy=kz+d)同除(z)可得

    (a'x+b'y=k+ frac{d}{z})

    但是(d<z)所以(frac{d}{z})是个分数,但是(a'x+b'y)都是整数,所以无解

    当然这个同时也有个推论。------https://blog.csdn.net/bcr_233/article/details/87897021
    (ax+by+cz+⋯+nm=f)(其中 (a,b,c,…,n,fa,b,c,)为整数),
    它有整数解的充要条件是 (f)(gcd(a,b,c,…,n))的整数倍。

    必要性就不证明了,充分性说一下,我们采用数学归纳法,对于有(n)个未知数:(a_1x_1+a_2x_2+...+a_nx_n)而言,我们已经证明了(a_1x_1+a_2x_2+...+a_{n-1}x_{n-1}=gcd(a_1,a_2,...))一定有解,那么我们就可以把前面的统一换成:(gcd(a_1,a_2,a_3...)x_{n+1}+a_nx_n),然后就变成了二个元,所以(gcd(a_1,a_2,a_3...,a_{n-1})x_{n+1}+a_nx_n=gcd(gcd(a_1,a_2,a_3...,a_{n-1}),a_n)=gcd(a_1,a_2,a_3...,a_n))必定有解,得证。


    至此,结束了(感觉又过了一个世纪)!

    #include<cstdio>
    #include<cstring>
    using  namespace  std;
    typedef  long  long  LL;
    LL  exgcd(LL  a,LL  b,LL  &x,LL  &y/*引用*/)
    {
    	if(a==0)
    	{
    		x=0;y=1;//0x+gcd(a,b)y=gcd(a,b) 
    		return  b;
    	}
    	else
    	{
    		LL  tx,ty;
    		LL  d=exgcd(b%a,a,tx,ty);
    		x=ty-(b/a)*tx;//公式 
    		y=tx;//公式 
    		return  d;
    	}
    }
    int  main()
    {
    	LL  a,b,c,x,y;scanf("%lld%lld%lld",&a,&b,&c);
    	LL  d=exgcd(a,b,x,y);
    	if(c%d!=0)printf("no solution!
    ");
    	else
    	{
    		c/=d;
    		printf("%lld %lld
    ",x*c,y*c);
    	}
    	return  0;
    }
    

    其实也没什么难度。

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

    同余方程(基本的EXGCD应用)

    题目描述
    
    【题意】
    已知a,b,m,求x的最小正整数解,使得ax≡b(mod m)
    【输入格式】
    一行三个整数 a,b,m。 1 ≤ a,b,m ≤ 10^9 
    【输出格式】
    一行一个整数x,无解输出"no solution!"
    【样例输入】
    2 5 7
    【样例输出】
    6
    

    我们会发现(ax≡b)((mod) (m))其实就是(ax-my=b)也就是(ax+my=b)(m)的正负无多大必要,反正只是求(x),而且(m>0)方便的地方在于(gcd)为正数)

    然后我们运用扩展欧几里德求出(x),但是如何求出最小的(x)

    学过不定式的人都知道,当(ax)([a,b])时,(my)也可以相应的加([a,b])时(([a,b])就是(a,b)的最小公倍数),(x,y)能得出另外一组整数解,也就是(x-[a,b]/a,y+[a,b]/b),又因为([a,b]=a*b/(a,b)),所以就是(x-b/(a,b),y-a/(a,b))((a,b))(gcd(a,b)))。

    而且([a,b])是能同时被(a,b)整除的最小的数字。

    所以,我们可以用(d=gcd(a,b),x=(x\%(b/d)+b/d)\%(b/d))得出(x)最小的正整数解

    #include<cstdio>
    #include<cstring>
    using  namespace  std;
    typedef  long  long  LL;
    LL  exgcd(LL  a,LL  b,LL  &x,LL  &y)
    {
    	if(a==0)
    	{
    		x=0;y=1;
    		return  b;
    	}
    	else
    	{
    		LL  tx,ty;
    		LL  d=exgcd(b%a,a,tx,ty);
    		x=ty-(b/a)*tx;
    		y=tx;
    		return  d;
    	}
    }
    LL  zbs(LL  x){return  x<0?-x:x;}
    int  main()
    {
    	LL  a,b,c,x,y;scanf("%lld%lld%lld",&a,&c,&b);
    	LL  d=exgcd(a,b,x,y);//求扩展欧几里得 
    	if(c%d!=0)printf("no solution!
    ");
    	else
    	{
    		c/=d;
    		x*=c;
    		LL  k=b/d;
    		x=((x%k)+k)%k;//最小的正整数解 
    		printf("%lld
    ",x);
    	}
    	return  0;
    }
    

    中国剩余定理CRT

    题目描述
    
    【题意】 
    同余方程是这样的:已知a,b,n,求x的最小正整数解,使得ax=b(mod m)
    同余方程组是这样:也是求x的最小正整数解,但已知b数组和m数组的情况下,
    x=b[1](mod m[1]),
    x=b[2](mod m[2]),
    x=b[3](mod m[3]),
    ~~~~~~
    x=b[n](mod m[n])
    m之间互质。
    【输入格式】
    一行一个整数 n(1<=n<=?)
    下来n行每行两个整数b[i],m[i]。 1 ≤b[i],m[i] ≤ 10^9 
    【输出格式】
    一行一个整数x,无解输出"no solution!"
    【样例输入】
    3
    3 5
    2 3
    2 7
    
    【样例输出】
    23
    

    我们考虑一下答案有可能是什么?

    (M)为所有(m)的总乘,可以发现答案其实是可以加减(M)的,所以我们就是把所有的方程组合并到一个模(M)的方程里面。

    那么我们只要找到一个一定是解的公式模(M)就行了

    就比如数据:

    3
    3 5
    2 3
    2 7
    

    转换一下思路,我们需要找到一个数字(a),是(3,7)的倍数,但是模(5)(3),一个数字(b)(5,7)的倍数,模(3)(2),一个数字(c)(5,3)的倍数,但是模(7)(2),那么(a+b+c)仔细想想就发现一定是其中的一组解了。

    但是怎么找到(a,b,c)呢,是(3,7)的倍数,模(5)(3)的一个数字?其实很简单的一个做法就是先找到是(3,7)的倍数,模(5)(1)的数字,然后乘以(3),对于前面的数字不就是在模(5)意义下,(lcm(3,7))(因为互质肯定为(3*7))的逆元。

    所以推到一下就会发现公式了,对于这道题目,我们设(M_i^{-1})表示模(m_i)意义下(M_i)的逆元,(M_i=frac{M}{m_{i}})那么公式就是:

    (x=(b_1M_1M_1^{-1}+...+b_nM_nM_n^{-1})\% M)

    这个东西最有用的地方在于它可以用于很多模数不为质数的算法中,可以把模数拆成各个不互质且更好计算的数字,最后合并。

    没有代码QMQ。

    扩展中国剩余定理exCRT

    题目描述
    
    【题意】 
    同余方程是这样的:已知a,b,n,求x的最小正整数解,使得ax=b(mod m)
    同余方程组是这样:也是求x的最小正整数解,但已知b数组和m数组的情况下,
    x=b[1](mod m[1]),
    x=b[2](mod m[2]),
    x=b[3](mod m[3]),
    ~~~~~~
    x=b[n](mod m[n])
    【输入格式】
    一行一个整数 n(1<=n<=?)
    下来n行每行两个整数b[i],m[i]。 1 ≤b[i],m[i] ≤ 10^9 
    【输出格式】
    一行一个整数x,无解输出"no solution!"
    【样例输入】
    3
    3 5
    2 3
    2 7
    
    【样例输出】
    23
    

    这道题就十分的友(e)善(xin)了

    现在我们可以列出(x=b_{1}+m_{1}y_{1},x=b_{2}+m_{2}y_{2}...)

    我们拿上面减下面的出:(m_{1}y_{1}-m_{2}y_{2}=b_{2}-b_{1})也就是(ax+by=c)的形式。

    那么我们写成(ax'+by'=c,a=m_{1},b=m_{2})(这里是(-m_{2})也无所谓,但是是正数后面算最小正整数比较方便。)(,c=b_{2}-b_{1},x'=y_{1},y'=y_{2})

    然后解出(x'),将(x')带入(x=b_{1}+ax')得出了(x)的其中一个解,而且我们又知道(x)要变只能加上或减去(b/(a,b)),所以(x)就只能加减((b/(a,b)*a=[a,b])),所以我们又可以列出一个新的方程:(x≡b_{1}+ax')((mod) ([a,b]))(其中的(b_{1}+ax')表示的是x的其中一组解),为了后面的方便,我们把这个方程代替为第二个方程,后面就是个二三做一次,三四做一次...

    然后我们就可以拿这个方程与第三个方程做这种事,重复如此,最后的(b_{?})就是答案,但是,别忘了是求最小整数解。


    求最小整数解有两种方法:

    我常用的方法:

    开局直接把每个(b_{i})(不包括求出来的(b_{?}))模(a_{i}),并且把每次求出的(x)都做一次最小整数解,那么得出来的新的方程的(b_{?})一定小于新的(m_{?}),为什么,(m_{?})([a,b]),而(b_{?})(b_{i}+ax)得出的。(这里的(a)表示(m_{i})(b)表示(m_{i+1})),我们知道(x≤b/(a,b)-1),那么(ax≤[a,b]-a)(b_{i}<a),那么(ax+b_{i}<[a,b]),所以(b_{?})小于(m_{?}),也就是不用模了。

    #include<cstdio>
    #include<cstring>
    using  namespace  std;
    typedef  long  long  LL;
    LL  exgcd(LL  a,LL  b,LL  &x,LL  &y)
    {
    	if(a==0)
    	{
    		x=0;y=1;
    		return  b;
    	}
    	else
    	{
    		LL  tx,ty;
    		LL  d=exgcd(b%a,a,tx,ty);
    		x=ty-(b/a)*tx;
    		y=tx;
    		return  d;
    	}
    }
    LL  a1,b1;
    bool  solve(LL  a2,LL  b2)
    {
    	LL  a=a1,b=a2,c=b2-b1,x,y;
    	LL  d=exgcd(a,b,x,y);
    	if(c%d!=0)return  false;
    	else
    	{
    		c/=d;
    		x*=c;
    		LL  lca=b/d;
    		x=((x%lca)+lca)%lca;//最小正整数解 
    		b1=a1*x+b1;a1=a*b/d;//不用模 
    		return  true;
    	}
    }
    int  main()
    {
    	int  n;scanf("%d",&n);
    	scanf("%lld%lld",&b1,&a1);b1%=a1;//开局模一模 
    	for(int  i=2;i<=n;i++)
    	{
    		LL  a,b;scanf("%lld%lld",&b,&a);b%=a;//开局模一模 
    		if(!solve(a,b))
    		{
    			printf("no solution!
    ");
    			return  0;
    		}
    	}
    	printf("%lld
    ",b1/*不用模*/);
    	return  0;
    }
    

    但是后面,我又发现了一个不是很麻烦的方法,就是结尾模一下就行了。。。

    //少了去正过程,中间会有许多数字变成负数
    //而且更容易爆long long
    #include<cstdio>
    #include<cstring>
    using  namespace  std;
    typedef  long  long  LL;
    LL  exgcd(LL  a,LL  b,LL  &x,LL  &y)
    {
    	if(a==0)
    	{
    		x=0;y=1;
    		return  b;
    	}
    	else
    	{
    		LL  tx,ty;
    		LL  d=exgcd(b%a,a,tx,ty);
    		x=ty-(b/a)*tx;
    		y=tx;
    		return  d;
    	}
    }
    LL  zbs(LL  x){return  x<0?-x:x;}
    LL  a1,b1;
    bool  solve(LL  a2,LL  b2)
    {
    	LL  a=a1,b=a2,c=b2-b1,x,y;
    	LL  d=exgcd(a,b,x,y);
    	if(c%d!=0)return  false;
    	else
    	{
    		c/=d;
    		x*=c;//可能为负数 
    		b1=a1*x+b1;a1=a*b/d;
    		return  true;
    	}
    }
    int  main()
    {
    	int  n;scanf("%d",&n);
    	scanf("%lld%lld",&b1,&a1);
    	for(int  i=2;i<=n;i++)
    	{
    		LL  a,b;scanf("%lld%lld",&b,&a);
    		if(!solve(a,b))
    		{
    			printf("no solution!
    ");
    			return  0;
    		}
    	}
    	LL  zhl=a1<0?-a1:a1;//绝对值 
    	printf("%lld
    ",(b1%zhl+zhl)%zhl/*尾巴模一模*/);
    	return  0;
    }
    

    以上便是作者的心得。

    一点想法

    当然,如果每个(x)前面也有系数(a_i)怎么办?我们求出的值可能不能整除于(x)的系数?

    在一定有解的情况下,我们可以拆开成两个方程:一个是(x≡bmod m)

    另外一个是(x≡0mod a)

    当然,这只是本蒟蒻的一点点假设,并没有实践过。

    要注意的地方

    这个EXCRT一定要特别注意正负问题,尤其是当数字特别大的时候,尽量做到每个地方都打模,才会比较保险。

  • 相关阅读:
    Kafka Replication: The case for MirrorMaker 2.0
    CentOs7.3 搭建 SolrCloud 集群服务
    Redis 集群:CLUSTERDOWN The cluster is down
    Redis Cluster: (error) MOVED
    Kafka重启出错:Corrupt index found
    redis cluster slots数量 为何是16384(2的14次方)
    redis监控工具汇总
    Redis运维利器 -- RedisManager
    redis三种集群策略
    Deploy custom service on non hadoop node with Apache Ambari
  • 原文地址:https://www.cnblogs.com/zhangjianjunab/p/12179941.html
Copyright © 2020-2023  润新知