浅谈扩展更相减损法在解二元一次不定方程领域的应用

watermouthhang · 2025-07-12 18:37:35 · 算法·理论

算法引入

在学习数论时，我们会学习一种算法，叫做扩展欧几里得算法，可以用于解二元一次不定方程，时间复杂度为 \operatorname O(\log(\max(a,b)) 。

不过，我们还可以找到另一种算法来解二元一次不定方程。
是谁呢？ 更相减损法（Alternating Subtraction Algorithm，ASA）。

更相减损法

本体

式子如下（想看证明的看这）：

所以，114 和 514 的最大公因数是 2 。

到省略号部分时，我们可以发现，当 a ， b 差距过大时，该算法效率会大幅降低。如果构造 \{n,1\} 这样的数据，复杂度会退化回 \operatorname O(n) 。

所以， \text{\color{red}需要优化} 。

优化·式子

显然， 当 a , b 都是 2 的倍数时

否则，当 a 是 2 的倍数时

否则，当 b 是 2 的倍数时

我们，就可以就此，利用移位运算（不会见此）除以 2 优化。

则最终式子为（其中 x \mid y 表示 y 是 x 的倍数，x \nmid y 表示 y 不是 x 的倍数）：

优化·时间复杂度

先说结论:稳定为 \operatorname O(\log(\max(a,b))) 。
为什么呢？

证明如下：

命题 1

内容：当递归不在 (2) 式时，以后也不可能进入 (2)式。

证明：对于 (1)(5) 式，显然。
对于 (3) 式，已经 2 \nmid b ，在运行完之后依旧如此，不会进入 (2) 式。
对于 (4) 式，同上。
对于 (6) 式，因为 2 \nmid a ， 2 \nmid b 。
即

所以

即 2 \mid (a-b) 。
因为 2 \nmid b , 2 \mid (a-b) ，
所以，在 (6) 式之后的下一个式子，一定是 (4)(5) 式。

综上，命题成立。

命题 2

内容： (1) 式只能出现一次。

证明：易证，略。

命题 3

内容：(5)(6) 式不能连续 4 次出现。

证明：对于 (5) 式，显然。
对于 (6) 式，见命题 1(6) ，成立。
对于 (5)(6) 式一起，由命题 1(6) ，（6）式以后为 (4)(5) 式。
在此时 (5) 式以后，只能为 (4)$式。

综上，命题成立。

证明

现在开始分析。

(1)(5)(6) 式不会连续出现，且附着在 (3)(4) 式后，可以忽略。
(2)(3)(4) 式最多运行 \log 级别的次数。

所以，时间复杂度为稳定的 \operatorname O(\log(\max(a,b))) 。

Q.E.D.

优化·代码部分

现在，我们根据这个函数，写下它的代码。

int gcdt(int a,int b)
{
    if(b==0)//(1)式 
    {
        return a;
    }
    while(!(a&1))//(3)式 
    {
        a>>=1;
    }
    while(!(b&1))//(4)式 
    {
        b>>=1;
    }
    /*
    注意，根据证明，(2)式不会出现，(3)(4)式
    显然不会同时出现 
    */
    if(a<b)//(5)式 
    {
        return gcdt(a,b-a);
    }
    return gcdt(b,a-b);//(6)式 
}

int gcd2(int a,int b)//真正的函数入口
{
    int ans=1;
    while(!(a&1||b&1))//先行处理(2)式 
    {
        ans<<=1;
        a>>=1;
        b>>=1;
    }
    return ans*gcdt(a,b);//不要忘了乘回去 
}

顺带一提，上文所述优化也被称作Stein算法，可用于计算整数的最大公因数。

真的开始：扩展·更相减损法

推导

我们所要解的，是以下的方程：

我们所要利用的，式刚刚推出来的式子：

为了方便，设 a',b' 为进入式子后的 a,b 值。

设 x',y' 满足

下面，开始解！方！程！

(1) 式

在其中， x=1,y=0 显然是一组解。

(2) 式

则 a'=a \div 2,b'=b \div 2 。
由于

所以

所以，x=x',y=y' 是一组特解。

(3) 式

则 a'=a \div 2,b'=b 。
由于

所以

所以，x=\frac{x'}{2},y=y' 是一组特解。

(4) 式

则 a'=a,b'=b \div 2 。
由于

所以

所以，x=x',y=\frac{y'}{2} 是一组特解。

(5) 式

则 a'=b,b'=a 。
由于

所以

所以，x=y',y=x' 是一组特解。

(6) 式

则 a'=b,b'=a-b 。
由于

所以

所以，x=y',y=x'-y' 是一组特解。

漏洞与解决

刚刚的推导看似无懈可击，但对于 (3)(4) 式，有致命的漏洞。

x,y都必须要是整数！

而在 (3)(4) 式中，出现了形如 \frac{x'}{2} 这样的式子。

例如，在 (3) 式中，当 2 \mid x' 时， \frac{x'}{2} 是整数，无事发生。

但是 2 \nmid x' 时，\frac{x'}{2} 不为整数，算法崩溃了
……吗？

在 (3) 式内，根据取值

设整数 r_1,r_2 满足

于是，括号内两式为整数。

由上文

上下两式作差得

此时右式中的分式部分已经是最简分数，为使 r_1 , r_2 为整数，令

（其中， n 是一个整数）

我们想要 \frac{x'}{2}+\frac{r_1}{2} 为整数，所以

因为 2\nmid x',2\nmid b,2\nmid\gcd(a\div 2,b) ，
所以

所以

为了 卡常 方便，令

由

得

满足上式。

所以

那么 (4) 式又如何呢？

因为推导过程和 (3) 式相同，所以过程省略，结果如下：

那么，针对 (3)(4) 式的“女娲补天”完成。

代码（例）

我们注意到，(5) 式后面一定递归 (6) 式。
所以，我们可以通过此种方式减少递归，一定程度上加快算法。

我们以扩欧模板题之一 P1082为例

依照扩展更相减损法，在 C++14(GCC9) 情况下，代码如下：

#include <cstdio>

using namespace std;
#define int long long
int ans1,ans2;

void exgcd_asa(int a,int b)
{
    if(b==0)
    {
        ans1=1,ans2=0;return;
    }
    if(!(a&1))
    {
        exgcd_asa(a>>=1,b);
        if(ans1&1)
        {
            ans1+=b;
            ans2-=a;
        }
        ans1>>=1;return;
    }
    if(!(b&1))
    {
        exgcd_asa(a,b>>=1);
        if(ans2&1)
        {
            ans1-=b;
            ans2+=a;
        }
        ans2>>=1;return;
    }
    if(a<b)
    {
        exgcd_asa(a,b-a);
        int temp=ans2;
        ans1=ans1-temp;return;
    }
    exgcd_asa(b,a-b);
    int temp=ans1;
    ans1=ans2;
    ans2=temp-ans1;
    return;
}
void exasa(int a,int b)
{
    while(!(a&1||b&1))
    {
        a>>=1,b>>=1;
    }
    exgcd_asa(a,b);
}

signed main()
{
    int a,b;
    scanf("%lld%lld",&a,&b);
    exasa(a,b);
    printf("%lld\n",(ans1%b+b)%b);
    return 0;
}

后记