使用扩展欧几里得算法求解整数二元一次不定方程
问题:形如ax+by=c(a,b均不为0)的方程,a,b,c都是整数,求(x,y)整数解。
整数二元一次方程有解的充要条件是gcd(a,b)|c。如果不能整除则无解。如果a,b其中一个是质数的话,那么gcd(a,b)=1肯定存在解。
简单的欧几里得算法如下,也称辗转相除法:
def gcd(a, b): while b != 0: a, b = b, a % b return a
扩展欧几里得算法可以求解ax+by=gcd(a,b). 求解到x,y之后,那么r/gcd(a,b) * x, r/gcd(a,b) * y就是方程ax+by=c的解。 扩展的推导过程是这样的:
- ax1+by1=gcd(a,b)=gcd(b,a%b)=bx2+(a%b)y2
- 假设a%b=r, a/b=k, 那么a%b=(a-bk). 带回到上面式子中
- ax1+by1=bx2+(a-bk)y2=ay2+b(x2-ky2)
- x1=y2, y1=(x2-ky2)
- 当a%b=0时,gcd(a,b)=b,x2=1,y2=0(其实可以为任意值?)
我们求解到x2,y2的基本情况之后,就可以倒退回x1,y1的情况
def ext_gcd(a, b, init_y2 = 0): if b == 0: return a, 1, init_y2 d, x2, y2 = ext_gcd(b, a % b, init_y2) x1, y1 = y2, x2 - (a // b) * y2 return d, x1, y1
测试了一下没有什么问题, y2的确可以为任意初始值
def test_ext_gcd(a, b, init_y2=0): d, x, y = ext_gcd(a, b, init_y2) assert a*x+b*y == d return d, x, y print(test_ext_gcd(10, 12)) print(test_ext_gcd(10, 12, 10)) print(test_ext_gcd(17, 31)) print(test_ext_gcd(17, 31, 20))