线性丢番图方程 --算法竞赛专题解析（21）：数论

本系列文章将于2021年整理出版。前驱教材：《算法竞赛入门到进阶》 清华大学出版社
网购：京东 当当   作者签名书：点我
公众号同步：算法专辑   
暑假福利：胡说三国
有建议请加QQ 群：567554289

   丢番图（Diophantus）是古希腊人，于公元前300年编写的《算术》，是最早的代数书。

1. 二元线性丢番图方程

   方程$ax+by=c$被称为二元线性丢番图方程，其中$a、b、c$是已知整数，$x、y$是变量，问是否有整数解。
   $ax+by=c$实际上是二维$x-y$平面上的一条直线，这条直线上如果有整数坐标点，方程就有解，如果没有整数坐标点，就无解。如果存在一个解，就有无穷多个解。
   下面的定理给出了有解的判断条件和通解的形式。
  定理¹：设a，b是整数且$gcd(a,b)=d$。如果$d$不能整除$c$，那么方程$ax+by=c$没有整数解，如果$d$能整除$c$，那么存在无穷多个整数解。另外，如果($x_0，y_0$)是方程的一个特解，那么所有的解（通解）可以表示为：
    $x=x_0+(b/d)n$
    $y=y_0-(a/d)n$
  其中$n$是任意整数。
  定理可以概况为：$ax+by=c$有解的充分必要条件是$d=gcd(a,b)$能整除$c$。

  例如：
  （1）方程18x + 3y = 7没有整数解，因为gcd(18, 3) = 3，3不能整除7；
  （2）方程25x + 15y = 70存在无穷个解，因为gcd(25, 15) = 5且5整除70，一个特解是$x_0$ = 4，$y_0$ = -2，通解是x = 4 + 3n，y = -2 - 5n。
  下面借助平面图解释定理。
  定理的前半部分：令$a=d{a}^{\prime }，b=d{b}^{\prime }$；有$ax+by=d(a^{\prime }x+b^{\prime }y)=c$；如果$x、y、a^{\prime }、b^{\prime }$都是整数，那么$c$必须是$d=gcd(a,b)$的倍数，才有整数解。

图1 二元线性丢番图方程

  定理的后半部分给出了通解的形式：$x$值按$b/d$递增，$y$值按 $-a/d$递增。设($x_0，y_0$)是一个格点（格点是指$x、y$坐标均为整数的点），移动到直线上另一个点($x_0+∆x，y_0+∆y$)，有$a∆x+b∆y=0$。$∆x$和$∆y$必须是整数，($x_0+∆x，y_0+∆y$)才是另一个格点。$∆x$ 最小是多少？因为$a/d$与$b/d$互素，只有$∆x=b/d，∆y=-a/d$时，$∆x$和$∆y$才是整数，并满足$a∆x+b∆y=0$。

  下面用一个例题来加强对定理的理解。

---

线段上的格点数量
题目描述：在二维平面上，给定两个格点p1 = (x1, y1)和p2 = (x2, y2)，问线段p1 p2上除了p1 、p2外还有几个格点？设x1 < x2。

---

题解：此题用暴力法逐一搜索格点复杂度太高，下面用丢番图方程的定理来求解。
思路：首先利用$p1、p2$把线段表示为方程$ax+by=c$的形式，它肯定有整数解。然后在线段范围内，根据$x$的通解的表达式$x=x_0+(b/d)n$（用$y=y_0-(a/d)n$也一样），当$x_1<x<x_2$时，求出$n$的取值情况有多少个，这就是线段内的格点数量。
  用$p_1、p_2$表示线段，经过化简，线段表示为：
    $(y_2-y_1)x+(x_1-x_2)y=y_2{x}_1-y_1{x}_2$
  对照$ax+by=c$，得$a=y_2-y_1，b=x_1-x_2，c=y_2{x}_1-y_1{x}_2，d=gcd(a,b)=gcd(|y_2-y_1|,|x_1-x_2|)$。
  对照通解公式$x=x_0+(b/d)n$，令特解是$x_1$，代入限制条件$x_1<x<x_2$，有：
    $x_1<x_1+(x_1-x_2/d)n<x_2$
  当 $-d<n<0$时满足上面的表达式，此时$n$有$d-1$种取值，即线段内有$d-1$个格点。
  下面是一个较难的题目。

---

Area poj 1265
题目描述：给出二维平面上的一个封闭的多边形，多边形的顶点都是格点。请计算多边形边界上格点数量$j$，内部格点数量$k$，以及多边形的面积$s$。

---

题解：边界上的格点$j$用前面的方法计算，面积$s$用几何叉积计算，最后内部的格点$k$通过Pick定理计算。
  Pick定理：在一个平面直角坐标系内，如果一个多边形的顶点都是格点，多边形的面积等于边界上格点数的一半加上内部格点数再减一，即$s=j/2+k-1$。

  求解方程$ax+by=c$的关键是找到一个特解。根据定理的描述，解和求GCD有关，所以求特解用到了欧几里得求GCD的思路，称为扩展欧几里得算法。

2. 扩展欧几里得算法

  方程$ax+by=gcd(a,b)$，根据前一节的定理，它有整数解。扩展欧几里得算法求一个特解$(x_0，y_0)$的代码如下²：

//返回 d = gcd(a,b); 并返回 ax + by = d的特解x,y
typedef long long ll;
ll extend_gcd(ll a,ll b,ll &x,ll &y){ 
    if(b == 0){ x=1; y=0; return a;}
    ll d = extend_gcd(b,a%b,y,x);
    y -= a/b * x;
    return d;
}

  有时候为了简化描述，把$ax+by=gcd(a,b)$两边除以$gcd(a,b)$，得到$cx+dy=1$，其中$c=a/gcd(a,b),d=b/gcd(a,b)$。$c，d$是互素的。$cx+dy=1$的通解是：$x=x_0+dn，y=y_0-cn$。

3. 二元丢番图方程$ax+by=c$的解

  用扩展欧几里德算法得到$ax+by=gcd(a,b)$的一个特解后，再利用它求方程$ax+by=c$的一个特解。步骤如下：
  （1）判断方程$ax+by=c$是否有整数解，即$gcd(a,b)$能整除$c$。记 $d=gcd(a,b)$。
  （2）用扩展欧几里得算法求$ax+by=d$的一个特解$x_0，y_0$。
  （3）在$a{x}_0+b{y}_0=d$两边同时乘以$c/d$，得：
    $a{x}_{0}c/d+b{y}_{0}c/d=c$
  （4）对照$ax+by=c$，得到它的一个解($x_0^{\prime },y_0^{\prime }$)是：
    $x_0^{\prime }=x_{0}c/d$
    $y_0^{\prime }=y_{0}c/d$
  （5）方程$ax+by=c$的通解：
     $x=x_0^{\prime }+(b/d)n$
     $y=y_0^{\prime }-(a/d)n$

4. 多元线性丢番图方程

  多元线性丢番图方程$a_1{x}_1+a_2{x}_2+...a_n{x}_n=c$在算法竞赛中很少见。为扩展思路，这里也给出介绍。
  定理：如果$a_1，a_2，...，a_n$，是非零整数，那么方程$a_1{x}_1+a_2{x}_2+...a_n{x}_n=c$有整数解，当且仅当$d=gcd(a_1，a_2，...，a_n)$整除$c$。如果存在一个解，则方程有无穷多个解。
  定理可以用数学归纳法证明。
  方程的计算步骤是：
  （1）判断方程是否有解。计算
     $d_2=gcd(a_1,a_2)$
     $d_3=gcd(d_2,a_3)$
     $d_4=gcd(d_3,a_4)$
     …
     $d_n=gcd(d_{n-1},a_n)$
  如果$d_n$能整除$c$，方程有解。如果有解，继续以下步骤。
   （2）求解。把方程分解为$n-1$个二元方程：
    $a_1{x}_1+a_2{x}_2=d_2{t}_2$
    $d_2{t}_2+a_3{x}_3=d_3{t}_3$
    …
    $d_{n-1}{t}_{n-1}+a_n{x}_n=c$
  然后从最后一个方程开始，依次往前求解。特解容易求得，通解的表达很麻烦。

---

1. 详细证明参考：《初等数论及其应用》102页。 ↩︎

2. 程序的执行过程参考：《算法导论》，Thomas H. Cormen等，机械工业出版社，31.2节。 ↩︎