本系列文章将于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))的倍数,才有整数解。
![](https://img-blog.csdnimg.cn/20200728175824601.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3dlaXhpbl80MzkxNDU5Mw==,size_16,color_FFFFFF,t_70)
定理的后半部分给出了通解的形式:(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)
然后从最后一个方程开始,依次往前求解。特解容易求得,通解的表达很麻烦。