关于约数,余数的算法
欧几里得算法-gcd
gcd(a,b)=gcd(b,a%b),其中a为整数,b是不为0的整数.
证明:
如果a<b,a%b=a,明显有,gcd(a,b)=gcd(b,a).
如果a>=b,设a=k*b+r.设d为a,b的任意公因数,因为d|a,d|b,所以d|(a-k*b),因此d也是b和a%b的公因数,反过来亦成立,因此a和b公因数集合和b和a%b的公因数集合相同,a和b的最大公因数即为b和a%b的最大公因数.
作用:求最大公约数
实现:
1 int gcd(int a,int b) 2 { 3 return b?gcd(b,a%b):a; 4 }
复杂度:O(log(a+b))
简单说明一下递归次数.
a=k(0)b+r(0)
b=k(1)r(0)+r(2)
...
其中r(0)>r(1)>r(2)...,k(i)>=0
这个过程中有r(i+2)<r(i)/2.证明如下:
若r(i+1)<r(i)/2,等式显然成立.
若r(i+1)>=r(i)/2,因为r(i)=k(i+1)r(i+1)+r(i+2),所以r(i+2)<r(i)/2.
即递归次数n=2*log2b=log2b=log(a+b)
扩展欧几里得算法-exgcd
作用:求整数解x和y使得ax+by=gcd(a,b)
原理:
给出方程ax+by=gcd(a,b)(这里假设b不为0),由公式可知有b*x1+a%b*y1=gcd(b,a%b).a%b可表示为a-a/b*b(这里的除是整除后向下取整),因此有:
1:ax+by=gcd(a,b)
2:b*x1+(a-a/b*b)*y1=gcd(b,a%b)
由gcd性质可知gcd(a,b)=gcd(b,a%b),因此有:
ax+by=b*x1+(a-a/b*b)*y1
把等式右边化简一下有:
ax+by=a*y1+b*(x1-a/b*y1)
因此有:
x=y1
y=x1-a/b*y1
这样不断递归求解下去,当递归到b为0时,x=1,y=0,再回溯一遍就能得到原方程的解了.当然扩展欧几里得在实现的过程中也能求出gcd(a,b).
1 int exgcd(int a,int b,int &x,int &y) 2 { 3 if(!b) 4 { 5 x=1,y=0; 6 return a; 7 } 8 int g=exgcd(b,a%b,y,x); //简化写法 9 y-=a/b*x; //y从下一个方程回溯回来时是下一个方程的x 10 return g; 11 }
复杂度:O(log(a+b))
扩展:事实上ax+by=gcd(a,b)的解有穷多个,而exgcd求出来的解只是其中一个特解.假设exgcd求出x0和y0,则方程解集为x=x0+k(b/gcd(a,b)),y=y0-k(a/gcd(a,b))---其中k为整数.把解集代入方程就能知道理由了.若要求ax+by=c的解x,y,先看gcd(a,b)能否整除c,如果不能整除的话此方程无解.如果能整除的话只需要把x0,y0扩大c/gcd(a,b)倍即可得到一个特解,然后把这个特解套入上面公式即可得出解集.
中国剩余定理
设m1,m2...mn是两两互质的整数,设m=∏(n,i=1)mi,Mi=m/mi,ti是方程Miti≡1(modmi)的一个解.对于任意的n个整数a1,a2...an,方程组
x≡a1(modm1)
x≡a2(modm2)
...
x≡an(modmn)
有整数解x=∑(n,i=1)aiMiti
证明:
对于第i个方程x≡ai(modmi)来说,ajMjtj(j!=i)%mi=0,因为Mj中包含mi.因此把x代入到这个方程中就只剩下aiMiti了,而ti又是Miti≡1(modmi)的解,因此对于第i个方程x=aiMiti≡ai(modmi)是成立的.以此类推可知x为该方程组的解.
中国剩余定理求出的只是该方程组的一个特解x0,该方程组的通解为x=x0+k*m---其中k为整数.
实现(省略exgcd的代码,变量参考上文):
1 LL solve(LL n) //n代表n个方程 2 { 3 LL M=1,ans=0,x,y; //这里的M表示上文的m 4 for(int i=0;i<n;i++) M*=m[i],a[i]%=m[i]; 5 for(int i=0;i<n;i++) 6 { 7 LL t=M/m[i]; 8 exgcd(t,m[i],x,y); //求上述Miti=1(modmi)的解ti 9 ans+=a[i]*t*x; 10 } 11 ans=((ans%M)+M)%M; //转化为最小整数解 12 return ans; 13 }
扩展:事实上中国剩余定理的应用条件有点苛刻,要求模数两两互质.因此这里讲一下更一般的同余方程组通用解法.
假设已经求出了前i个方程的特解ans(i),这里说明一下lcm(k)=lcm(m1,m2...mk),则通解x(i)=ans(i)+k*lcm(i),把x(i)代入到第i+1个方程可得ans(i)+k*lcm(i)≡a(i+1)(mod m(i+1)),化简一下可得k*lcm(i)≡a(i+1)-ans(i)(mod m(i+1)),这里可以运用扩展欧几里得判一下有无解,有解的话可以求出特解k(特)=k0,通解为k(通)=k0+k*m(i+1)/gcd(lcm(i),m(i+1))---在扩展欧几里得篇已提及.
把k(特)代入x(i)=ans(i)+k*lcm(i)得到ans(i+1)=ans(i)+k0*lcm(i)
把k(通)代入x(i)=ans(i)+k*lcm(i)得到x(i+1)=ans(i+1)+k*lcm(i+1)
这样递推下去即可求出整个方程组的通解.
实现(省略exgcd的代码,变量参考上文):
1 LL solve(LL n) //方程组由n个方程 2 { 3 LL ans=0,lcm=1,x,y; 4 for(int i=0;i<n;i++) 5 { 6 if(!i) ans=a[i],lcm*=m[i]; //只有一个方程时的解 7 else 8 { 9 LL g=exgcd(lcm,m[i],x,y),c=a[i]-ans; 10 if(c%g) //判定有无解 11 { 12 ans=-1; 13 break; 14 } 15 x*=c/g; //别忘了扩大解,因为扩展欧几里得求得的是c为g时的解 16 LL t=m[i]/g; 17 x=(x%t+t)%t; //保证解x为最小非负整数解,防爆ans 18 ans+=x*lcm; 19 lcm*=m[i]/g; //求新的lcm 20 ans=((ans%lcm)+lcm)%lcm; //与上面x同理,防爆ans 21 } 22 } 23 return ans; 24 }
