线性丢番图方程Linear Diophantine equation,名字听得非常高大上的。但是实际上很简单的,线性丢番图就是下面这种形式的方程:
a
x
+
b
y
=
c
ax+by=c
ax+by=c
公式中a,b,c都是已知数,x和y是未知数,然后求x和y的整数解。这种有两个未知数,却只有一个方程的情况,解不止一个。这种解不止一个的情况,一般是用字母t或k作为可变的整数来表示通解的形式。首先我们不考虑a和b为零的场景,存在0的场景下,这个方程就没什么意义了。
举个例子:
9
x
+
100
y
=
1
9x+100y=1
9x+100y=1
其中一个解为
x
=
−
11
y
=
1
9
×
−
11
+
100
×
1
=
1
x=-11\\ y=1\\ 9\times-11+100\times 1=1
x=−11y=19×−11+100×1=1
线性丢番图方程和前面我写的扩展欧几里得法的方程形式很类似。扩展欧几里得法方程是以下形式:
a
x
+
b
y
=
g
ax+by=g
ax+by=g
扩展欧几里得法方程可以说是线性丢番图方程的特殊形式,因为g是a和b的最大公约数。而线性丢番图方程的c是一个更一般的值(后面会讲实际上是g的整数倍)。
一 与同余的联系(觉得难可以略过)
首先这个方程组可以转变为同余方程组:
a
x
≡
c
(
m
o
d
b
)
b
x
≡
c
(
m
o
d
a
)
ax\equiv c(mod \ b)\\ bx\equiv c(mod \ a)\\
ax≡c(mod b)bx≡c(mod a)
这个步骤我不能跳过,所以推导一下:
∵
a
x
=
c
−
b
y
∴
c
=
a
x
−
b
y
(
c
m
o
d
b
)
=
(
a
x
−
b
y
)
(
m
o
d
b
)
=
(
a
x
m
o
d
b
)
∴
c
≡
a
x
(
m
o
d
b
)
\because ax = c-by\\ \therefore c = ax - by\\ (c \ mod \ b) = (ax - by) (mod \ b)= (ax \ mod \ b)\\ \therefore c \equiv ax (mod \ b)
∵ax=c−by∴c=ax−by(c mod b)=(ax−by)(mod b)=(ax mod b)∴c≡ax(mod b)
当a、b互质的时候,有下列公式:
x
≡
c
a
−
1
(
m
o
d
b
)
x \equiv ca^{-1}(mod \ b)
x≡ca−1(mod b)
a
−
1
a^{-1}
a−1是
a
a
a模b的模乘逆元,也就是
a
⋅
a
−
1
≡
1
(
m
o
d
b
)
a\cdot a^{-1}\equiv 1(mod\ b)
a⋅a−1≡1(mod b)。上式的证明很简单,两边同时乘以a,有:
a
x
≡
a
c
a
−
1
(
m
o
d
b
)
≡
a
a
−
1
c
(
m
o
d
b
)
≡
1
⋅
c
(
m
o
d
b
)
ax \equiv aca^{-1}(mod \ b)\equiv aa^{-1}c(mod \ b)\equiv 1\cdot c(mod \ b)
ax≡aca−1(mod b)≡aa−1c(mod b)≡1⋅c(mod b)
当a、b不互质的时候,只有c能被a与b的最大公约数整除时才有解,有:
g
=
g
c
d
(
a
,
b
)
x
≡
c
g
(
a
g
)
−
1
(
m
o
d
b
g
)
y
=
c
−
a
x
b
g=gcd(a,b)\\ x\equiv\frac c g (\frac a g)^{-1}(mod \ \frac b g) \\ y=\frac{c-ax}{b}
g=gcd(a,b)x≡gc(ga)−1(mod gb)y=bc−ax
所以说线性丢番图方程可以转换为线性同余方程组,用中国剩余定理去解,但是我们不需要这么做。
二 计算方法
线性丢番图方程
a
x
+
b
y
=
c
ax+by=c
ax+by=c,有c能被a与b的最大公约数整除时才有解。这个的证明超级简单啊。设g是a和b的最大公约数,那么ax是g的倍数,by也是g的倍数,那么c一定要是g的倍数,否则没有整数解。而且有多个解,但是通解是特解的同余类。这个有点像微分方程,通解等于特解加上一个可变的部分。
举个例子,
6
x
+
9
y
=
21
6x+9y=21
6x+9y=21,其中一个特解是
x
=
2
,
y
=
1
x=2,y=1
x=2,y=1,通解是
x
=
3
t
+
2
,
y
=
−
2
t
+
1
x=3t+2,y=-2t+1
x=3t+2,y=−2t+1,可以看出通解中包含了特解。
在前面,我们学过扩展欧几里得法可以计算如下的方程:
a
x
+
b
y
=
g
ax+by=g
ax+by=g
扩展欧几里得法(如果没学过可以看我的文章《3 模乘逆元-扩展欧几里得法》)要求g必须是a和b的最大公约数,否则不能计算。而对于方程
a
x
+
b
y
=
c
ax+by=c
ax+by=c,c不一定是a或b的最大公约数g,但一定是g的整数倍。所以我们可以这样计算:
第一步,先用扩展欧几里得法解出一组数
x
g
x_g
xg和
y
g
y_g
yg,也就是解下面的方程:
a
x
g
+
b
y
g
=
g
ax_g+by_g=g
axg+byg=g
第二步,因为我们知道c是g的整数倍,上面的方程两边同时乘以
c
g
\frac c g
gc:
a
x
g
c
g
+
b
y
g
c
g
=
g
c
g
=
c
ax_g\frac c g+by_g\frac c g=g\frac c g=c
axggc+byggc=ggc=c
所以我们就找到了一个特解:
x
0
=
x
g
c
g
y
0
=
y
g
c
g
x_0=x_g\frac c g\\ y_0=y_g\frac c g
x0=xggcy0=yggc
那特解呢?我们把通解代入方程然后再加上俩个特殊的数字
a
(
x
0
+
k
b
g
)
+
b
(
y
0
−
k
a
g
)
=
a
x
0
+
b
y
0
+
k
a
b
g
−
k
a
b
g
a(x_0+\frac{kb}{g})+b(y_0-\frac{ka}{g})=ax_0+by_0+\frac{kab}g-\frac{kab}g
a(x0+gkb)+b(y0−gka)=ax0+by0+gkab−gkab
所以特解就是(k为任意整数):
x
=
x
0
+
k
b
g
y
=
y
0
−
k
a
g
x=x_0+\frac{kb}{g}\\ y=y_0-\frac{ka}{g}
x=x0+gkby=y0−gka
读者会诧异为什么要除以g,虽然不除于g也是解,但是除以g是为了找出所有的解,只有包含了所有的解才能叫通解all solutions嘛。
三 代码实现
我用python实现的,但是我把以前的扩展欧几里得法代码修改了下,增加了一个返回值最大公约数g。因为无论是校验方程是否可解,还是求特解和表示通解,这个g都特别重要。
def extend_euclid(a: int, b: int):
if a < b:
x, y, g = extend_euclid(b, a)
return y, x, g
# 保存a/b的系数
coefficients = []
while b != 0:
coefficients.append(a // b)
a, b = b, a % b
# 循环代入
x, y = 1, 0
coefficients.reverse()
for i in coefficients:
x, y = y, x - y * i
return x, y, a
def diophantine(a, b, c):
x_g, y_g, g = extend_euclid(a, b)
if c % g != 0:
raise 'No solutions'
return x_g * c // g, y_g * c // g, b // g, -a // g
if __name__ == '__main__':
x0, y0, xk, yk = diophantine(6, 9, 27)
print(f'solution is x ={x0} {"+" if xk > 0 else ""}{xk}t, y = {y0} {"+" if yk > 0 else ""}{yk}t')
测试结果:
solution is x = -9 +3t, y = 9 -2t
测试结果完全正确,哦,耶,挑战下一数论题目。