定义
本文探讨的不定方程指的是二元一次不定方程。
exgcd,扩展欧几里得算法用于求解二元一次不定方程(Diophantine Equation)的整数解。
扩展欧几里得算法(Exgcd)
普通解
求出以下不定方程的一组解:
a
x
+
b
y
=
c
ax+by=c
ax+by=c
令
g
=
gcd
(
a
,
b
)
g=\gcd(a,b)
g=gcd(a,b):
根据裴蜀定理,方程有整数解当且仅当
g
∣
c
g|c
g∣c。
我们只需要求出
a
x
+
b
y
=
g
ax+by=g
ax+by=g的一组解
(
x
0
,
y
0
)
(x_0,y_0)
(x0,y0),令
x
0
=
c
g
x
0
,
y
0
=
c
g
y
0
x_0=\frac c g x_0,y_0=\frac c g y_0
x0=gcx0,y0=gcy0,就得到了原方程的一组解。
现在,我们求方程 a x + b y = g ax+by=g ax+by=g的一组解:
设方程
a
x
+
b
y
=
g
ax+by=g
ax+by=g的一组解用
f
(
a
,
b
)
f(a,b)
f(a,b)表示。
这里的
g
g
g表示的是参数里的
a
,
b
a,b
a,b的
gcd
(
a
,
b
)
\gcd(a,b)
gcd(a,b),不一定是一开始的
g
g
g。
一种显然的想法是, f ( a , b ) = ( x 0 , y 0 ) f(a,b)=(x_0,y_0) f(a,b)=(x0,y0)与 f ( a ′ , b ′ ) = ( x ′ , y ′ ) f(a',b')=(x',y') f(a′,b′)=(x′,y′)具有某种关系。
为了使关系好求,我们令
a
′
=
b
,
b
′
=
a
%
b
a'=b,b'=a\%b
a′=b,b′=a%b,这样的话两个方程的右边是相等的:
f
(
a
,
b
)
对应的方程
:
a
x
0
+
b
y
0
=
gcd
(
a
,
b
)
f(a,b)对应的方程:\;\;\;\;\,ax_0+by_0=\gcd(a,b)
f(a,b)对应的方程:ax0+by0=gcd(a,b)
f
(
b
,
a
%
b
)
对应的方程
:
b
x
′
+
(
a
%
b
)
y
′
=
gcd
(
b
,
a
%
b
)
=
gcd
(
a
,
b
)
f(b,a\%b)对应的方程:bx'+(a\%b)y'=\gcd(b,a\%b)=\gcd(a,b)
f(b,a%b)对应的方程:bx′+(a%b)y′=gcd(b,a%b)=gcd(a,b)
现在处理一下这个式子:
b
x
′
+
(
a
%
b
)
y
′
=
b
x
′
+
(
a
−
b
⌊
a
b
⌋
)
y
′
=
a
y
′
+
b
(
x
′
−
y
′
⌊
a
b
⌋
)
bx'+(a\%b)y'=bx'+(a-b\lfloor\frac a b\rfloor)y'=ay'+b\left(x'-y'\lfloor\frac a b\rfloor\right)
bx′+(a%b)y′=bx′+(a−b⌊ba⌋)y′=ay′+b(x′−y′⌊ba⌋)
现在对比一下系数:
f
(
b
,
a
%
b
)
:
a
y
′
+
b
(
x
′
−
y
′
⌊
a
b
⌋
)
f(b,a\%b):ay'+b\left(x'-y'\lfloor\frac a b\rfloor\right)
f(b,a%b):ay′+b(x′−y′⌊ba⌋)
f
(
a
,
b
)
:
a
x
0
+
b
y
0
f(a,b):\;\,\;\;\;ax_0+by_0
f(a,b):ax0+by0
显然有: f ( a , b ) = ( x 0 , y 0 ) = ( y ′ , x ′ − y ′ ⌊ a b ⌋ ) f(a,b)=(x_0,y_0)=\left(y',x'-y'\lfloor\frac a b\rfloor\right) f(a,b)=(x0,y0)=(y′,x′−y′⌊ba⌋)
因此可以递归计算
f
(
b
,
a
%
b
)
f(b,a\%b)
f(b,a%b)。
显然有结束条件:
f
(
b
,
0
)
=
(
1
,
0
)
f(b,0)=(1,0)
f(b,0)=(1,0)
注意结束时的 y y y理论上可以任意取值,但是可以用归纳法证明,如果这里 y y y取0,则中间过程的运算中 ∣ x ∣ , ∣ y ∣ ≤ max ( ∣ a ∣ , ∣ b ∣ ) |x|,|y|\leq \max(|a|,|b|) ∣x∣,∣y∣≤max(∣a∣,∣b∣),确保不会爆int或long long。
事实上,如果 y y y取0,可以保证在中间过程中: ∣ x ∣ < ∣ b 2 g ∣ , ∣ y ∣ < ∣ a 2 g ∣ ( a > b ) |x|<\left|\frac b{2g}\right|,|y|<\left|\frac a{2g}\right|(a>b) ∣x∣< 2gb ,∣y∣< 2ga (a>b)
通解
刚才我们已经求出了方程 a x + b y = c ax+by=c ax+by=c的一组解 ( x 0 , y 0 ) (x_0,y_0) (x0,y0),现在我们考虑一下求出它的通解 ( x , y ) (x,y) (x,y)的形式。
我们注意到:
a
x
0
+
b
y
0
=
c
ax_0+by_0=c
ax0+by0=c,把左式分为左右两部分:
左半部分
(
a
x
0
)
+
右半部分
(
b
y
0
)
=
c
左半部分(ax_0)+右半部分(by_0)=c
左半部分(ax0)+右半部分(by0)=c
如果在特解的基础上构造一组通解,要么左半部分增加,同时右半部分减少,要么左半部分减少,同时右半部分增加。
这里令左半部分增加
d
d
d:
(
a
x
0
+
d
)
+
(
b
y
0
−
d
)
=
c
(ax_0+d)+(by_0-d)=c
(ax0+d)+(by0−d)=c
因此:
a
(
x
0
+
d
a
)
+
b
(
y
0
+
d
b
)
=
c
a(x_0+\frac d a)+b(y_0+\frac d b)=c
a(x0+ad)+b(y0+bd)=c,通解
(
x
,
y
)
(x,y)
(x,y)就是
(
x
0
+
d
a
,
y
0
+
d
b
)
(x_0+\frac d a,y_0+\frac d b)
(x0+ad,y0+bd)。
这就要求
d
d
d是a、b的倍数,只有是倍数才能整除。
也就是说
d
d
d是a、b的公倍数。
同时
d
d
d应该尽可能小,如果存在更小的
d
d
d,使得有的解用现在有
d
d
d无法表示出来,就会漏解。因此
d
d
d是a、b的最小公倍数,即
l
c
m
(
a
,
b
)
lcm(a,b)
lcm(a,b)
考虑到
d
d
d可以加很多次,设
l
c
m
=
l
c
m
(
a
,
b
)
lcm=lcm(a,b)
lcm=lcm(a,b),通解就是:
(
x
0
+
l
c
m
a
⋅
k
,
y
0
−
l
c
m
b
⋅
k
)
(
k
∈
Z
)
\left(x_0+\frac {lcm} a\cdot k,y_0-\frac {lcm} b\cdot k\right)(k\in \mathbb{Z})
(x0+alcm⋅k,y0−blcm⋅k)(k∈Z)
考虑到
l
c
m
lcm
lcm往往很大,可能会爆long long,又有
gcd
(
a
,
b
)
⋅
l
c
m
(
a
,
b
)
=
a
⋅
b
\gcd(a,b)\cdot lcm(a,b)=a\cdot b
gcd(a,b)⋅lcm(a,b)=a⋅b,因此还有一种通解的形式是:
(
x
0
+
b
g
⋅
k
,
y
0
−
a
g
⋅
k
)
(
k
∈
Z
)
\left(x_0+\frac {b} g\cdot k,y_0-\frac {a} g\cdot k\right)(k\in \mathbb{Z})
(x0+gb⋅k,y0−ga⋅k)(k∈Z)
特殊解
在实际应用中,我们往往需要带有特殊性质的解。
我们可以定义一个 l , r l,r l,r,其中 l l l表示 x x x能取到正整数的最小的 k k k, r r r表示 y y y能取到正整数的最大的 k k k,这样就围成了一个区间,区间之内是 x , y x,y x,y都能取到正整数解的 k k k值,如果 l > r l>r l>r,说明方程无正整数解。
如果 l ≤ r l\leq r l≤r,则有正整数解:
- 正整数解的个数就是区间长度,即 r − l + 1 r-l+1 r−l+1。
- k = l k=l k=l时取到 x x x的最小正整数值,由于 x x x减小, y y y要增大,对应必定是 y y y的最大值(在正整数解范围内)。
- k = r k=r k=r时取到 y y y的最小正整数值,由于 y y y减小, x x x要增大,对应必定是 x x x的最大值(在正整数解范围内)。
l
,
r
l,r
l,r可以直接列不等式求出:
{
x
0
+
l
c
m
a
⋅
l
>
0
y
0
−
l
c
m
b
⋅
r
>
0
\left\{\begin{matrix} x_0+\frac {lcm}{a}\cdot l>0\\ y_0-\frac {lcm}{b}\cdot r>0 \end{matrix}\right.
{x0+alcm⋅l>0y0−blcm⋅r>0
解得:
{
l
>
−
x
0
⋅
a
l
c
m
r
<
y
0
⋅
b
l
c
m
\left\{\begin{matrix} l>-\frac{x_0\cdot a}{lcm}\\ r<\frac{y_0\cdot b}{lcm } \end{matrix}\right.
{l>−lcmx0⋅ar<lcmy0⋅b
注意到这个范围不容易取整,因此修改一下方程组:
{
x
0
+
l
c
m
a
⋅
l
≥
1
y
0
−
l
c
m
b
⋅
r
≥
1
\left\{\begin{matrix} x_0+\frac {lcm}{a}\cdot l\geq 1\\ y_0-\frac {lcm}{b}\cdot r\geq 1 \end{matrix}\right.
{x0+alcm⋅l≥1y0−blcm⋅r≥1
就有:
{
l
≥
(
1
−
x
0
)
⋅
a
l
c
m
r
≤
(
y
0
−
1
)
⋅
b
l
c
m
\left\{\begin{matrix} l\geq \frac {(1-x_0)\cdot a}{lcm}\\ r\leq \frac {(y_0-1)\cdot b}{lcm} \end{matrix}\right.
{l≥lcm(1−x0)⋅ar≤lcm(y0−1)⋅b
然后怎么取
l
,
r
l,r
l,r呢?其实
l
l
l应该取到最小整点,
r
r
r取到最大整点,因此有:
l
=
⌈
(
1
−
x
0
)
⋅
a
l
c
m
⌉
l=\lceil\frac {(1-x_0)\cdot a}{lcm}\rceil
l=⌈lcm(1−x0)⋅a⌉
r
=
⌊
(
y
0
−
1
)
⋅
b
l
c
m
⌋
r=\lfloor\frac {(y_0-1)\cdot b}{lcm}\rfloor
r=⌊lcm(y0−1)⋅b⌋
后记
于是皆大欢喜。