算法竞赛专题解析(21):数论--线性丢番图方程

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

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

1. 二元线性丢番图方程

   方程(ax + by = c)被称为二元线性丢番图方程,其中(a、b、c)是已知整数,(x、y)是变量,问是否有整数解。
   (ax + by = c)实际上是二维(x-y)平面上的一条直线,这条直线上如果有整数坐标点,方程就有解,如果没有整数坐标点,就无解。如果存在一个解,就有无穷多个解。
   下面的定理给出了有解的判断条件和通解的形式。
  定理[1]:设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 = da',b = db');有(ax + by = d(a' x + b' y) = c);如果(x、y、a'、b')都是整数,那么(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))的代码如下[2]

//返回 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)在(ax_0 + by_0 = d)两边同时乘以(c/d),得:
    (a x_0 c/d + b y_0 c/d = c)
  (4)对照(ax + by = c),得到它的一个解((x_0', y_0'))是:
    (x_0' = x_0 c / d)
    (y_0' = y_0 c / d)
  (5)方程(ax + by = c)的通解:
     (x = x_0' + (b/d)n)
     (y = y_0' - (a/d)n)

4. 多元线性丢番图方程

  多元线性丢番图方程(a_1x_1 + a_2x_2 + ... a_nx_n = c)在算法竞赛中很少见。为扩展思路,这里也给出介绍。
  定理:如果(a_1,a_2,...,a_n),是非零整数,那么方程(a_1x_1 + a_2x_2 + ... a_nx_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_1x_1 + a_2x_2 = d_2 t_2)
    (d_2t_2 + a_3x_3 = d_3 t_3)
    ...
    (d_{n-1}t_{n-1} + a_nx_n = c)
  然后从最后一个方程开始,依次往前求解。特解容易求得,通解的表达很麻烦。


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

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

原文地址:https://www.cnblogs.com/luoyj/p/13394962.html