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

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

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

1. 二元线性丢番图方程

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

  例如:
  (1)方程18x + 3y = 7没有整数解,因为gcd(18, 3) = 3,3不能整除7;
  (2)方程25x + 15y = 70存在无穷个解,因为gcd(25, 15) = 5且5整除70,一个特解是 x 0 x_0 x0 = 4, y 0 y_0 y0 = -2,通解是x = 4 + 3n,y = -2 - 5n。
  下面借助平面图解释定理。
  定理的前半部分:令 a = d a ′ , b = d b ′ a = da',b = db' a=dab=db;有 a x + b y = d ( a ′ x + b ′ y ) = c ax + by = d(a' x + b' y) = c ax+by=d(ax+by)=c;如果 x 、 y 、 a ′ 、 b ′ x、y、a'、b' xyab都是整数,那么 c c c必须是 d = g c d ( a , b ) d = gcd(a, b) d=gcd(a,b)的倍数,才有整数解。

图1 二元线性丢番图方程

  定理的后半部分给出了通解的形式: x x x值按 b / d b/d b/d递增, y y y值按 − a / d - a/d a/d递增。设( x 0 , y 0 x_0,y_0 x0y0)是一个格点(格点是指 x 、 y x、y xy坐标均为整数的点),移动到直线上另一个点( x 0 + ∆ x , y 0 + ∆ y x_0 + ∆x,y_0 + ∆y x0+xy0+y),有 a ∆ x + b ∆ y = 0 a∆x + b∆y = 0 ax+by=0 ∆ x ∆x x ∆ y ∆y y必须是整数,( x 0 + ∆ x , y 0 + ∆ y x_0 + ∆x,y_0 + ∆y x0+xy0+y)才是另一个格点。 ∆ x ∆x x 最小是多少?因为 a / d a/d a/d b / d b/d b/d互素,只有 ∆ x = b / d , ∆ y = − a / d ∆x = b/d,∆y = - a/d x=b/dy=a/d时, ∆ x ∆x x ∆ y ∆y y才是整数,并满足 a ∆ x + b ∆ y = 0 a∆x + b∆y = 0 ax+by=0

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


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


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


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


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

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

2. 扩展欧几里得算法

  方程 a x + b y = g c d ( a , b ) ax + by = gcd(a, b) ax+by=gcd(a,b),根据前一节的定理,它有整数解。扩展欧几里得算法求一个特解 ( x 0 , y 0 ) (x_0,y_0) (x0y0)的代码如下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;
}

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

3. 二元丢番图方程 a x + b y = c ax + by = c ax+by=c的解

  用扩展欧几里德算法得到 a x + b y = g c d ( a , b ) ax + by = gcd(a, b) ax+by=gcd(a,b)的一个特解后,再利用它求方程 a x + b y = c ax + by = c ax+by=c的一个特解。步骤如下:
  (1)判断方程 a x + b y = c ax + by = c ax+by=c是否有整数解,即 g c d ( a , b ) gcd(a,b) gcd(a,b)能整除 c c c。记 d = g c d ( a , b ) d = gcd(a,b) d=gcd(a,b)
  (2)用扩展欧几里得算法求 a x + b y = d ax + by = d ax+by=d的一个特解 x 0 , y 0 x_0,y_0 x0y0
  (3)在 a x 0 + b y 0 = d ax_0 + by_0 = d ax0+by0=d两边同时乘以 c / d c/d c/d,得:
     a x 0 c / d + b y 0 c / d = c a x_0 c/d + b y_0 c/d = c ax0c/d+by0c/d=c
  (4)对照 a x + b y = c ax + by = c ax+by=c,得到它的一个解( x 0 ′ , y 0 ′ x_0', y_0' x0,y0)是:
     x 0 ′ = x 0 c / d x_0' = x_0 c / d x0=x0c/d
     y 0 ′ = y 0 c / d y_0' = y_0 c / d y0=y0c/d
  (5)方程 a x + b y = c ax + by = c ax+by=c的通解:
     x = x 0 ′ + ( b / d ) n x = x_0' + (b/d)n x=x0+(b/d)n
     y = y 0 ′ − ( a / d ) n y = y_0' - (a/d)n y=y0(a/d)n

4. 多元线性丢番图方程

  多元线性丢番图方程 a 1 x 1 + a 2 x 2 + . . . a n x n = c a_1x_1 + a_2x_2 + ... a_nx_n = c a1x1+a2x2+...anxn=c在算法竞赛中很少见。为扩展思路,这里也给出介绍。
  定理:如果 a 1 , a 2 , . . . , a n a_1,a_2,...,a_n a1a2...an,是非零整数,那么方程 a 1 x 1 + a 2 x 2 + . . . a n x n = c a_1x_1 + a_2x_2 + ... a_nx_n = c a1x1+a2x2+...anxn=c有整数解,当且仅当 d = g c d ( a 1 , a 2 , . . . , a n ) d = gcd(a_1,a_2,...,a_n) d=gcd(a1a2...an)整除 c c c。如果存在一个解,则方程有无穷多个解。
  定理可以用数学归纳法证明。
  方程的计算步骤是:
  (1)判断方程是否有解。计算
     d 2 = g c d ( a 1 , a 2 ) d_2 = gcd(a_1, a_2) d2=gcd(a1,a2)
     d 3 = g c d ( d 2 , a 3 ) d_3 = gcd(d_2, a_3) d3=gcd(d2,a3)
     d 4 = g c d ( d 3 , a 4 ) d_4 = gcd(d_3, a_4) d4=gcd(d3,a4)
     …
     d n = g c d ( d n − 1 , a n ) d_n = gcd(d_{n-1}, a_n) dn=gcd(dn1,an)
  如果 d n d_n dn能整除 c c c,方程有解。如果有解,继续以下步骤。
   (2)求解。把方程分解为 n − 1 n - 1 n1个二元方程:
     a 1 x 1 + a 2 x 2 = d 2 t 2 a_1x_1 + a_2x_2 = d_2 t_2 a1x1+a2x2=d2t2
     d 2 t 2 + a 3 x 3 = d 3 t 3 d_2t_2 + a_3x_3 = d_3 t_3 d2t2+a3x3=d3t3
    …
     d n − 1 t n − 1 + a n x n = c d_{n-1}t_{n-1} + a_nx_n = c dn1tn1+anxn=c
  然后从最后一个方程开始,依次往前求解。特解容易求得,通解的表达很麻烦。


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

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

评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

罗勇军

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值