2.4 线性丢番图方程

  线性丢番图方程Linear Diophantine equation，名字听得非常高大上的。但是实际上很简单的，线性丢番图就是下面这种形式的方程：
$$ax+by=c$$
  公式中a，b，c都是已知数，x和y是未知数，然后求x和y的整数解。这种有两个未知数，却只有一个方程的情况，解不止一个。这种解不止一个的情况，一般是用字母t或k作为可变的整数来表示通解的形式。首先我们不考虑a和b为零的场景，存在0的场景下，这个方程就没什么意义了。
  举个例子：
$$9x+100y=1$$
  其中一个解为
$$\begin{gathered} x=-11 \\ y=1 \\ 9\times -11+100\times 1=1 \end{gathered}$$
  线性丢番图方程和前面我写的扩展欧几里得法的方程形式很类似。扩展欧几里得法方程是以下形式：
$$ax+by=g$$
  扩展欧几里得法方程可以说是线性丢番图方程的特殊形式，因为g是a和b的最大公约数。而线性丢番图方程的c是一个更一般的值（后面会讲实际上是g的整数倍）。

一 与同余的联系（觉得难可以略过）

  首先这个方程组可以转变为同余方程组：
$$\begin{gathered} ax\equiv c(modb) \\ bx\equiv c(moda) \end{gathered}$$
  这个步骤我不能跳过，所以推导一下：
$$\begin{gathered} \because ax=c-by \\ \therefore c=ax-by \\ (cmodb)=(ax-by)(modb)=(axmodb) \\ \therefore c\equiv ax(modb) \end{gathered}$$
  当a、b互质的时候，有下列公式：
$$x\equiv c{a}^{-1}(modb)$$
  $a^{-1}$是$a$模b的模乘逆元，也就是$a\cdot a^{-1}\equiv 1(modb)$。上式的证明很简单，两边同时乘以a，有：
$$ax\equiv ac{a}^{-1}(modb)\equiv a{a}^{-1}c(modb)\equiv 1\cdot c(modb)$$
  当a、b不互质的时候，只有c能被a与b的最大公约数整除时才有解，有：
$$\begin{gathered} g=gcd(a,b) \\ x\equiv \frac{c}{g}(\frac{a}{g}{)}^{-1}(mod\frac{b}{g}) \\ y=\frac{c-ax}{b} \end{gathered}$$
  所以说线性丢番图方程可以转换为线性同余方程组，用中国剩余定理去解，但是我们不需要这么做。

二 计算方法

  线性丢番图方程$ax+by=c$，有c能被a与b的最大公约数整除时才有解。这个的证明超级简单啊。设g是a和b的最大公约数，那么ax是g的倍数，by也是g的倍数，那么c一定要是g的倍数，否则没有整数解。而且有多个解，但是通解是特解的同余类。这个有点像微分方程，通解等于特解加上一个可变的部分。
  举个例子，$6x+9y=21$，其中一个特解是$x=2,y=1$，通解是$x=3t+2,y=-2t+1$，可以看出通解中包含了特解。
  在前面，我们学过扩展欧几里得法可以计算如下的方程：
$$ax+by=g$$
  扩展欧几里得法（如果没学过可以看我的文章《3 模乘逆元-扩展欧几里得法》）要求g必须是a和b的最大公约数，否则不能计算。而对于方程$ax+by=c$，c不一定是a或b的最大公约数g，但一定是g的整数倍。所以我们可以这样计算：
  第一步，先用扩展欧几里得法解出一组数$x_g$和$y_g$，也就是解下面的方程：
$$a{x}_g+b{y}_g=g$$
  第二步，因为我们知道c是g的整数倍，上面的方程两边同时乘以$\frac{c}{g}$：
$$a{x}_g\frac{c}{g}+b{y}_g\frac{c}{g}=g\frac{c}{g}=c$$
  所以我们就找到了一个特解：
$$\begin{gathered} x_0=x_g\frac{c}{g} \\ {y}_0=y_g\frac{c}{g} \end{gathered}$$
  那特解呢？我们把通解代入方程然后再加上俩个特殊的数字
$$a(x_0+\frac{kb}{g})+b(y_0-\frac{ka}{g})=a{x}_0+b{y}_0+\frac{kab}{g}-\frac{kab}{g}$$
  所以特解就是（k为任意整数）：
$$\begin{gathered} x=x_0+\frac{kb}{g} \\ y=y_0-\frac{ka}{g} \end{gathered}$$
  读者会诧异为什么要除以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

  测试结果完全正确，哦，耶，挑战下一数论题目。