非再帰で拡張ユークリッドの互除法を書く

n 番煎じだと思いますが，既存の説明で私はあまりよく理解できていなかったので刺さる人がいればいいなと思って書きました．

問題

整数 $a$ , $b$ が与えられる． $a$ と $b$ の最大公約数 $g=\gcd(a, b)$ と $ax+by=g$ を満たす整数 $x$ , $y$ を求めよ．

再帰解法

$a$ を $b$ で割ったときの商と余りをそれぞれ $p$ , $q$ とおく．このとき $a=bp+q$ であり，方程式に代入して $(1)$ を得る．

$\displaystyle \begin{align} ax+by&=g\\ (bp+q)x+by&=g\\ b(px+y)+qx&=g\tag{1} \end{align}$

$\gcd(b,q)=\gcd(a,b)$ なので， $a\leftarrow b$ , $b\leftarrow q$ として再帰的にこの問題を解けば元の解 $(x,\ y)$ が復元できます．具体的には， $bx'+qy'=g$ を満たす整数 $x'$ , $y'$ に対し， $x=y'$ , $y=x'-py'$ とすればよいです．

再帰の終端条件は $b=0$ であり，このとき $(x,\ y,\ g)=(1,\ 0,\ g)$ とすればよいです．また，再帰の深さは $O(\log\min(a,b))$ であることが示されます．この辺りの説明は qiita.com が詳しいので分らない場合は参照してください．

非再帰解法

上の再帰解法において，再帰の深さが $i$ のときに解く方程式を $a_i x+b_i y=g$ ，解を $(x,\ y)=(x_i,\ y_i)$ とします．また， $a_i=b_i p_i + q_i$ とします． $b _ i\neq 0$ のとき，上述の解の復元方法から，次の漸化式が立ちます．

$\displaystyle \begin{cases} x_i=y_{i+1}\\ y_i=x_{i+1}-p_i y_{i+1} \end{cases}$

また，再帰の終端の深さを $n$ とすると， $x_n=1$ , $y_n=0$ です．よくある漸化式と異なり，添え字の大きい方から $(x_i,\ y_i)$ 求まっていくことに注意です．一方で， $a_i$ , $b_i$ は次のように添え字の小さい方から求まっていきます．このことが非再帰で拡張ユークリッドの互除法を書くことを難しくしています．

$\displaystyle \begin{cases} a_i=b_{i-1}\\ b_i=q_{i-1} \end{cases}$

ここで， $(x_i,\ y_i)$ に対する遷移を行列で表現すると以下のようになります．

$\displaystyle \begin{pmatrix} 0&1\\ 1&-p_i \end{pmatrix} \begin{pmatrix} x_{i+1}\\ y_{i+1} \end{pmatrix}= \begin{pmatrix} x_i\\ y_i \end{pmatrix}$

これを $i=0,\ldots,n-1$ に対して適用することで，次の $(2)$ を得ます．

$\displaystyle \begin{pmatrix} 0&1\\ 1&-p_0 \end{pmatrix} \begin{pmatrix} 0&1\\ 1&-p_1 \end{pmatrix}\cdots \begin{pmatrix} 0&1\\ 1&-p_{n-1} \end{pmatrix} \begin{pmatrix} x_n\\ y_n \end{pmatrix} = \begin{pmatrix} x_0\\ y_0 \end{pmatrix} \tag{2}$

$a_0=a$ , $b_0=b$ であり， $x_0$ , $y_0$ は $ax+by=g$ を満たす整数解であること，加えて $x_n=1$ , $y_n=0$ であることに注意します．

$a_i$ , $b_i$ , $p_i$ は添え字の小さい順に求まっていくので，行列を左から順にかけていく *1 ことで，めでたく非再帰で $(x,\ y)$ が求められます．

実装

C++17 で動くと思います．バグってたらすみません．

利便性のため求まる $\gcd$ が非負になるように補正をかけていますが，必須ではありません．

// return {{x,y},g} s.t. ax+by=g=gcd(a,b)>=0. 
std::pair<std::pair<long long, long long>, long long> ext_gcd(long long a, long long b) {
    /** 計算中の行列を
     *   x y
     *   z w
     *  としている
     */
    long long x = 1, y = 0;
    long long z = 0, w = 1;
    long long tmp;
    while (b) {
        long long p = a / b, q = a % b;
        tmp = x - y * p; x = y; y = tmp;
        tmp = z - w * p; z = w; w = tmp;
        a = b; b = q;
    }
    // gcd(a,b)>=0 となるように補正を掛ける (任意)
    if (a >= 0) return {{x, z}, a};
    else return {{-x, -z}, -a};
}

*1: 行列の乗算は結合則を満たすのでこのようなことが出来ます