欧几里得算法功能是求两个正整数a和b的最大公因数。函数名一般设为gcd,利用的性质是gcd(a,b) = gcd(a%b, b)。
简单证明:不妨设a = kb+r, d为a,b公因数, d|a, d|b。r = a-kb,所以d|r。而r = a%b,d|a%b, d|b, 所以d = gcd(a%b,b) = gcd(a, b)
这样不断递归下去,a%b=0时,下次调用变成gcd(b,0),最大公因数是b,递归停止。
不过要保证a>b,很麻烦,只需要交换一下顺序,gcd(a,b) = gcd(b,a%b),这样即使a<b,也可以在下次递归调用交换ab顺序:
int gcd(int a, int b){ return b == 0 ? a : gcd(b, a%b); }
然而实际上a,b往往会出现负数,gcd的结果也可能出现负数。(涉及到取模取余的很烦的问题)
扩展欧几里得除了能求a,b的最大公因数之外,能求解方程ax+by=c的整数解。
方程变形为m*gcd(a,b)*x+n*gcd(a,b)=c, gcd(a,b)(m*x+n*y)=c。显然当c为gcd(a,b)的倍数时方程有整数解,否则无整数解。
考虑问题ax+by=gcd(a,b)
叫扩展欧几里得当然是因为是欧几里得的扩展。利用性质gcd(a,b) = gcd(a%b, b):
ax1+by1=gcd(a,b) ------①
b*x2+a%b*y2=gcd(b,a%b)------②
等式右侧两个数相等,所以有:ax1+by1 = b*x2+a%b*y2 = b*x2+(a-[a/b]*b)*y2 = a*y2+b*(x2-[a/b]*y2)
由恒等关系知,x1=y2, y1=x2-[a/b]*y2。就是说,①式的解x1,y1,可以由②式的解x2,y2得到,这样就形成了递归。
终止条件是b=0,此时x=1,y=0,最大公因数是a。
总结一下就是:扩展欧几里得在求a,b的最大公因数时,“顺便”算出了ax+by=gcd(a,b)的一组解。
int exgcd(int a, int b, int& x, int& y){ if(b == 0){ x = 1; y = 0; return a; } int gcd = exgcd(b, a%b, x, y); int tmp = x; x = y; y = tmp - a/b*y; return gcd; }
调用的时候先声明x,y,调用exgcd(a, b, x, y)即可把一组解写到x,y里,并且函数返回a和b的最大公因数。
现在得到了一组特解x1,y1,然而其他解呢?设另一组解x2,y2. ax1+by1=ax2+by2, a(x1-x2)=b(y2-y1), 设gcd(a,b)等于g,A等于a/g, B等于b/g, (AB互素) 等式两边同时除以g, 即A(x1-x2)=B(y2-y1), 所以x1-x2=kB, y2-y1=kA
所以通解可以表示为:x=x1-kB, y=y1+kA.
想得到最小的正数x呢?根据x通解可知,只涉及正数时最小x是x1%B。然而x1和B都可能为负数。B为负数时,B取绝对值。如果x1为负,只需x1%B+B即可得到最小x。
扩展欧几里得个很有用的用途: 求乘法逆元
ax≡1(modm) 时 x叫a的逆元。逆元作用大大地有,比如取模的时候:我们知道(a*b)%m = a%m*b%m, 但(a/b)%m不能拆开取模。例如计算组合公式C的时候,阶乘除以阶乘取模,必须一边算一边取模才能不爆范围,这时就要用到逆元:
inv表示逆元,(a/b)%m = a%m*inv(b)%m。
现在求这个逆元x,ax+my=1,这时祭出扩展欧几里得就能搞定 (ps:m经常是1e9+7,怕是为了保证a和m互素)