【数论算法】扩展欧几里得算法详解+编程求解不定方程

扩展欧几里得算法详解及不定方程求解模板题(洛谷P5656)解题思路

题目链接

https://www.luogu.com.cn/problem/P5656

题目描述

给定不定方程

a x + b y = c ax+by=c ax+by=c

若该方程无整数解,输出 − 1 -1 1

若该方程有整数解,且有正整数解,则输出其正整数解的数量,所有正整数解中 x x x 的最小值,所有正整数解中 y y y 的最小值,所有正整数解中 x x x 的最大值,以及所有正整数解中 y y y 的最大值。

若方程有整数解,但没有正整数解,你需要输出所有整数解中 x x x 的最小正整数值, y y y 的最小正整数值。

正整数解即为 x , y x,y x,y 均为正整数的解, 0 0 0 不是正整数。

整数解即为 x , y x,y x,y 均为整数的解。

x x x 的最小正整数值即所有 x x x 为正整数的整数解中 x x x 的最小值, y y y 同理。

数据范围

共有 T T T组数据, 1 ≤ T ≤ 2 × 1 0 5 1\le T \le 2 \times 10^5 1T2×105 1 ≤ a , b , c ≤ 1 0 9 1\le a,b,c \le 10^9 1a,b,c109

扩展欧几里得算法

本题是扩展欧几里得算法的模板题,下面先介绍扩展欧几里得算法。

代码模板

给定 a a a b b b,下面的扩展欧几里得算法可以找到一组数 x x x y y y d d d,满足 a x + b y = d ax+by=d ax+by=d 并且 d = gcd ⁡ ( a , b ) d=\gcd(a,b) d=gcd(a,b)

void ext_gcd(long long a, long long b, long long &d, long long &x, long long &y){
    if (b == 0){
        d = a;
        x = 1; 
        y = 0;
    }else{
        ext_gcd(b, a % b, d, y, x);// (*)
        y -= a / b * x;
    }
}

时间复杂度

上面算法的时间复杂度为 O ( log ⁡ ( max ⁡ ( a , b ) ) ) O(\log(\max(a,b))) O(log(max(a,b))),且最坏情况是计算斐波那契数列的相邻两项。

算法正确性的证明

用数学归纳法可以证明算法的正确性,即算法对于任意合法的自然数 a a a b b b,均可以找到一组数 x x x y y y d d d,满足 a x + b y = d ax+by=d ax+by=d 并且 d = gcd ⁡ ( a , b ) d=\gcd(a,b) d=gcd(a,b)

归纳基础

b = 0 b=0 b=0 时, gcd ⁡ ( a , 0 ) = a \gcd(a,0)=a gcd(a,0)=a,取 x = 1 , y = 0 x=1,y=0 x=1,y=0,有 a x + b y = a × 1 + 0 × 0 = a = gcd ⁡ ( a , 0 ) = gcd ⁡ ( a , b ) ax+by=a\times 1+0\times 0=a=\gcd(a,0)=\gcd(a,b) ax+by=a×1+0×0=a=gcd(a,0)=gcd(a,b) 成立。

归纳步

假设第 i i i 次递归之后,代码中的 (*) 行已经得到了正确的结果,即返回值 y ′ y' y x ′ x' x 满足:

b y ′ + ( a   m o d   b ) x ′ = gcd ⁡ ( b , a   m o d   b ) = gcd ⁡ ( a , b ) = d by'+(a \bmod b)x'=\gcd(b, a\bmod b)=\gcd(a,b)=d by+(amodb)x=gcd(b,amodb)=gcd(a,b)=d

则经过步骤 y − = ⌊ a / b ⌋ × x ′ y-= \lfloor a/b\rfloor\times x' y=a/b×x 之后,新的数对 ( x , y ) (x,y) (x,y) 满足:

a x + b y = a x ′ + b ( y ′ − ⌊ a / b ⌋ x ′ ) = a x ′ + b y ′ − b ⌊ a / b ⌋ x ′ ax+by=ax'+b(y'-\lfloor a/b\rfloor x')=ax'+by'-b\lfloor a/b\rfloor x' ax+by=ax+b(ya/bx)=ax+byba/bx
= ( a − b ⌊ a / b ⌋ ) x ′ + b y ′ = ( a   m o d   b ) x ′ + b y ′ = d = gcd ⁡ ( a , b ) =(a-b\lfloor a/b\rfloor)x'+by'=(a\bmod b)x'+by'=d=\gcd(a,b) =(aba/b)x+by=(amodb)x+by=d=gcd(a,b)

由数学归纳法得知,任意步骤之后,上面的算法都会返回正确的结果,即都会找到一组数 x x x y y y d d d,满足 a x + b y = d ax+by=d ax+by=d 并且 d = gcd ⁡ ( a , b ) d=\gcd(a,b) d=gcd(a,b)

算法时间复杂度的证明

我们分几种情况考虑算法的执行。

b = 0 b=0 b=0

到达递归出口,执行三条赋值语句,这三条语句的时间复杂度为 O ( 1 ) O(1) O(1)

b ≠ 0 b \neq 0 b=0 a < b a < b a<b

这种情况需要递归。传入 (*) 行递归调用的前两个参数 a ′ = b a'=b a=b b ′ = a   m o d   b = a b'=a \bmod b=a b=amodb=a, 这时有 a ′ = b > b ′ = a a'=b > b'=a a=b>b=a, 递归后 a ′ a' a b ′ b' b 的关系将变为下一种 a ′ ≥ b ′ a' \ge b' ab 的情况。递归调用之外的时间开销为 O ( 1 ) O(1) O(1)

b ≠ 0 b \neq 0 b=0 a ≥ b a \ge b ab

这种情况也需要递归。传入 (*) 行递归调用的前两个参数 a ′ = b a'=b a=b b ′ = a   m o d   b b'=a \bmod b b=amodb。由于 b ′ = a   m o d   b < b = a ′ b'=a \bmod b < b=a' b=amodb<b=a,因此递归后 a ′ a' a b ′ b' b 的关系依然为 a ′ ≥ b ′ a' \ge b' ab 的情况。注意到,当 a ≥ b a \ge b ab 时, a   m o d   b ≤ ⌊ a 2 ⌋ a \bmod b \le \lfloor \frac{a}{2}\rfloor amodb2a (证明在下面)。 因此我们有下一次递归调用的前两个参数 a ′ ′ a'' a b ′ ′ b'' b 满足 a ′ ′ = b ′ = a   m o d   b ≤ ⌊ a 2 ⌋ a'' = b' = a \bmod b \le \lfloor \frac{a}{2}\rfloor a=b=amodb2a, 进而有
max ⁡ ( a ′ ′ , b ′ ′ ) = a ′ ′ ≤ ⌊ a 2 ⌋ = ⌊ max ⁡ ( a , b ) 2 ⌋ \max(a'', b'')=a''\le \lfloor \frac{a}{2}\rfloor = \lfloor \frac{\max(a,b)}{2}\rfloor max(a,b)=a2a=2max(a,b)

因此在 a ≥ b a \ge b ab 的情况下,每两次递归调用就可以让递归函数前两个参数 a a a b b b 的最大值减少至少一半。而递归调用之外的时间开销为 O ( 1 ) O(1) O(1), 因此算法总的时间复杂度为 O ( log ⁡ ( max ⁡ ( a , b ) ) ) O(\log(\max(a,b))) O(log(max(a,b)))

引理: ∀ a , b ∈ Z + \forall a,b \in \mathbb{Z}^{+} a,bZ+,若 a ≥ b a \ge b ab,则有 a   m o d   b ≤ ⌊ a 2 ⌋ a \bmod b \le \lfloor \frac{a}{2}\rfloor amodb2a

b ≤ ⌊ a 2 ⌋ + 1 b \le \lfloor\frac{a}{2}\rfloor+1 b2a+1, 则 a   m o d   b ≤ b − 1 = ⌊ a 2 ⌋ a \bmod b \le b - 1=\lfloor \frac{a}{2}\rfloor amodbb1=2a,引理成立;

b > ⌊ a 2 ⌋ + 1 b > \lfloor\frac{a}{2}\rfloor+1 b>2a+1, 由于 a ≥ b a \ge b ab,因而有 ⌊ a b ⌋ ≥ 1 \lfloor\frac{a}{b}\rfloor \ge 1 ba1,因此

a   m o d   b = a − ⌊ a b ⌋ b ≤ a − b < a − ⌊ a 2 ⌋ − 1 a \bmod b = a - \lfloor\frac{a}{b}\rfloor b \le a - b < a -\lfloor\frac{a}{2}\rfloor-1 amodb=ababab<a2a1
= 2 ⌊ a 2 ⌋ + ( a   m o d   2 ) − ⌊ a 2 ⌋ − 1 = ⌊ a 2 ⌋ + ( a   m o d   2 ) − 1 ≤ ⌊ a 2 ⌋ =2\lfloor\frac{a}{2}\rfloor+(a\bmod 2)-\lfloor\frac{a}{2}\rfloor-1=\lfloor\frac{a}{2}\rfloor+(a\bmod 2)-1\le \lfloor\frac{a}{2}\rfloor =22a+(amod2)2a1=2a+(amod2)12a

因而引理成立。

解题步骤

运用扩展欧几里得算法

对于题目输入的三个数 a , b , c a,b,c a,b,c,首先可以通过扩展欧几里得算法求出一组数 x 0 , y 0 , d x_0,y_0,d x0,y0,d,满足 a x 0 + b y 0 = d a x_0 + b y_0 =d ax0+by0=d d = gcd ⁡ ( a , b ) d=\gcd(a,b) d=gcd(a,b)

原方程是否有整数解?

首先检查 c c c 是不是 d d d 的倍数,如果不是,则说明无解(因为方程左边 a x + b y ax+by ax+by 一定是 d d d 的倍数,但右边 c c c 不是 d d d 的倍数)。

如何求出原方程的一组特解?

接着将方程 x 0 a + y 0 b = d x_0 a + y_0 b=d x0a+y0b=d 两边同乘以 c d \frac{c}{d} dc,得到 a x 0 c d + b y 0 c d = d c d = c ax_0 \frac{c}{d}+by_0 \frac{c}{d}=d\frac{c}{d}=c ax0dc+by0dc=ddc=c,所以 ( x 0 c d , y 0 c d ) (x_0\frac{c}{d},y_0\frac{c}{d}) (x0dc,y0dc) 是不定方程 a x + b y = c ax+by=c ax+by=c 的一组解。

如何通过特解找到其他的解?

假设有两组解 ( x 1 , y 1 ) (x_1,y_1) (x1,y1) ( x 2 , y 2 ) (x_2,y_2) (x2,y2) 使得原不定方程成立,那么有 a x 1 + b y 1 = a x 2 + b y 2 = c ax_1+by_1=ax_2+by_2=c ax1+by1=ax2+by2=c,因此 a ( x 1 − x 2 ) = b ( y 2 − y 1 ) a(x_1-x_2)=b(y_2-y_1) a(x1x2)=b(y2y1),该等式两边同时除以 d = gcd ⁡ ( a , b ) d=\gcd(a,b) d=gcd(a,b),得到:

a d ( x 1 − x 2 ) = b d ( y 2 − y 1 ) \frac{a}{d}(x_1-x_2)=\frac{b}{d}(y_2-y_1) da(x1x2)=db(y2y1)

此时 gcd ⁡ ( a d , b d ) = 1 \gcd(\frac{a}{d},\frac{b}{d})=1 gcd(da,db)=1,因此如果要让等式成立, x 1 − x 2 x_1-x_2 x1x2 必须要为 b d \frac{b}{d} db 的整数倍。不妨设 x 1 − x 2 = t b d x_1-x_2=t\frac{b}{d} x1x2=tdb,那么有 y 2 − y 1 = t a d y_2-y_1=t\frac{a}{d} y2y1=tda

因此给定一组特解 ( x , y ) (x,y) (x,y),我们可以将不定方程所有可能的解表示为 ( x + t b d , y − t a d ) (x+t\frac{b}{d},y-t\frac{a}{d}) (x+tdb,ytda),其中 t t t 为任意整数。结合上一小节的分析,我们可以将这一题中的所有可能的解表示为 ( x 0 c d + t b d , y 0 c d − t a d ) (x_0\frac{c}{d}+t\frac{b}{d},y_0\frac{c}{d}-t\frac{a}{d}) (x0dc+tdb,y0dctda),其中 t t t 为任意整数。

如何找到 x x x 为尽可能小的正整数的解?

我们已经有了一组初始解 ( x 0 c d , y 0 c d ) (x_0\frac{c}{d},y_0\frac{c}{d}) (x0dc,y0dc)。要找到 x x x 为尽可能小的正整数的解,实际上就是要找到 t x m i n t_{x_{min}} txmin,使得 x 0 c d + t x m i n b d x_0\frac{c}{d}+t_{x_{min}}\frac{b}{d} x0dc+txmindb 为尽可能小的正整数 x m i n x_{min} xmin,即要让 x m i n = x 0 c d + t x m i n b d x_{min}=x_0\frac{c}{d}+t_{x_{min}}\frac{b}{d} xmin=x0dc+txmindb 为尽可能小的正整数。

注意到 1 ≤ x m i n ≤ b d 1 \le x_{min} \le \frac{b}{d} 1xmindb (若 x m i n < 1 x_{min}<1 xmin<1,则 x m i n x_{min} xmin 不为正数;若 x m i n > b d x_{min}> \frac{b}{d} xmin>db,则可以让 x m i n x_{min} xmin 减小为一个更小的解。因此让 x 0 c d + t x m i n b d = x m i n x_0\frac{c}{d}+t_{x_{min}}\frac{b}{d}=x_{min} x0dc+txmindb=xmin 两边同时整除以 b d \frac{b}{d} db,可以得到 ⌊ x 0 c / d b / d ⌋ + t x m i n = 0 \lfloor\frac{x_0c/d}{b/d}\rfloor+t_{x_{min}}=0 b/dx0c/d+txmin=0 或者 ⌊ x 0 c / d b / d ⌋ + t x m i n = 1 \lfloor\frac{x_0c/d}{b/d}\rfloor+t_{x_{min}}=1 b/dx0c/d+txmin=1,即 t x m i n = − ⌊ x 0 c / d b / d ⌋ t_{x_{min}}=-\lfloor\frac{x_0c/d}{b/d}\rfloor txmin=b/dx0c/d 或者 t x m i n = 1 − ⌊ x 0 c / d b / d ⌋ t_{x_{min}}=1-\lfloor\frac{x_0c/d}{b/d}\rfloor txmin=1b/dx0c/d

因此,我们计算出 x 0 c d + ( − ⌊ x 0 c / d b / d ⌋ ) b d x_0\frac{c}{d}+(-\lfloor\frac{x_0c/d}{b/d}\rfloor)\frac{b}{d} x0dc+(b/dx0c/d)db 以及 x 0 c d + ( 1 − ⌊ x 0 c / d b / d ⌋ ) b d x_0\frac{c}{d}+(1-\lfloor\frac{x_0c/d}{b/d}\rfloor)\frac{b}{d} x0dc+(1b/dx0c/d)db,看一下哪一个是 x x x 最小的正整数值。

求出 x x x 最小的正整数值 x m i n x_{min} xmin 之后,我们可以进一步求出当前的 y = c − a x b y=\frac{c-ax}{b} y=bcax,将其记为 y m a x y_{max} ymax

情况一:原方程不存在正整数解

如果 y m a x y_{max} ymax 小于等于 0 0 0,则说明原不定方程不存在 x x x y y y 均为正整数的解(因为要使得方程 a x + b y = c ax+by=c ax+by=c 成立, y y y 变大时需要让 x x x 变小,但 x m i n x_{min} xmin 再变小就不为正整数了,因此 x m i n x_{min} xmin 无法变小, y m a x y_{max} ymax 也无法变大)。此时可以用类似于上一小节的方法求出 y m i n y_{min} ymin,然后可以输出 x m i n x_{min} xmin y m i n y_{min} ymin

情况二:原方程存在正整数解

如果 y m a x y_{max} ymax 大于 0 0 0,则 ( x m i n , y m a x ) (x_{min},y_{max}) (xmin,ymax) 为原方程的一组正整数解,且 x m i n x_{min} xmin 为可以达到的最小整数 x x x 值, y m a x y_{max} ymax 为可以达到的最大整数 y y y 值。用类似的方法可以求出 ( x m a x , y m i n ) (x_{max},y_{min}) (xmax,ymin)。进而我们可以求出解的数目为 n = x m a x − x m i n b / d + 1 n=\dfrac{x_{max}-x_{min}}{b/d}+1 n=b/dxmaxxmin+1,即 n = t x m a x − t x m i n + 1 n=t_{x_{max}}-t_{x_{min}}+1 n=txmaxtxmin+1。然后将 n , x m i n , y m i n , x m a x , y m a x n,x_{min},y_{min},x_{max},y_{max} n,xmin,ymin,xmax,ymax 均输出即可。

代码

模板题代码就不贴了,理解了上一节解题步骤后按照解题步骤调用模板并实现算术运算与条件判断等语句即可。

相关思考题

  1. 如果要解含有三个未知数的不定方程 a x + b y + c z = e ax+by+cz=e ax+by+cz=e,该怎样处理呢?其中 x , y , z x,y,z x,y,z 为未知数, a , b , c , e a,b,c,e a,b,c,e 为输入的常数。
    • 基本思路是先引入一个新的未知数 w w w,解未知数为 x , w x,w x,w 的不定方程 a x + g c d ( b , c ) w = e ax + gcd(b, c)w = e ax+gcd(b,c)w=e,然后再解不定方程 b y + c z = g c d ( b , c ) by + cz = gcd(b, c) by+cz=gcd(b,c)
  2. 如果要解含有 n n n 个未知数的不定方程,该怎样处理呢?
    • 基本思路类似上一题,引入一个新的未知数转化为解含有 n − 1 n-1 n1 个未知数的不定方程,然后可以递归处理。
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值