拡張ユークリッドの互除法についてはここ[drken2019]を参照してください．
特に断りがなければ，数式中のすべての記号は整数として扱います．
$ax+by=\gcd(a,b)$ を満たす $(x,y)$ を求めるアルゴリズムについて，drken2019の再帰関数実装から非再帰版実装を考えるための方法を記述する．
drken2019での実装を以下のようになっている．
// 返り値: a と b の最大公約数
// ax + by = gcd(a, b) を満たす (x, y) が格納される
long long extGCD(long long a, long long b, long long &x, long long &y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
long long d = extGCD(b, a%b, y, x);
y -= a/b * x;
return d;
}
この実装では $x,y$ を参照としているため，状態の更新順序が把握しづらくなっている．
これをもう少し数学っぽく書き換える．
具体的には拡張ユークリッド互除法を行う関数を $f:\mathbb{Z}^2\mapsto\mathbb{Z}^3$ とし，$(g,x,y)=f(a,b)$ となる関数に見えるようにする．
#include <stdint.h> // int64_t
#include <utility> // tuple
auto f(int64_t a, int64_t b)->std::tuple<int64_t,int64_t,int64_t>{
if(b==0) return {a, 1, 0};
const auto [g,x,y] = f(b, a%b);
return {g, y, x-a/b*y};
}
参照が無くなったことで，処理の依存が少しほどけてわかりやすくなった気がする．
再帰関数をループに展開するときに考えなければならないことは，スタックに情報を持たなければならないかどうかである．
つまり，callerがcalleeから戻ってきたときに保持している情報を用いて計算を再開するかどうかである．
この処理は $f(a,b)\to f(b,a\mod b)\to \cdots\to f(\gcd(a,b),0)\to\cdots\to f(b,a\mod b)\to f(a,b)$ のように動作するが，このとき戻ってきた先で保存しておかなければならない情報はreturn {g,y,x-a/b*y};を見てわかる通り引数である．
つまり最終的に何らかの方法を用いて x-a/b*yが前進的に計算できれば非再帰化できるということになる．
最終的に必要な$(x,y)$のために，各ステップ $i$ における剰余を初期の$a,b$についての線形結合で表す変数$r_i$を導入する．
ユークリッドの互除法の手続きは以下のようになっている．
\begin{align*}
a&=q_0b+r_2\\
b&=q_1r_2+r_3\\
&\ \ \vdots\\
\gcd(a,b)&=q_kr_k+0
\end{align*}
これを一般化すると $r_{i-1}=q_ir_i+r_{i+1}$ となり，移項して $r_{i+1}=r_{i-1}-q_ir_i$ となる．
$a,b$ は剰余ではないが，最初に式に出てくるため初期状態として採用し$r_0=a, r_1=b$とする．
$r_i$は任意の$i$について $ax+by=r_i$ で表すことができるので(ユークリッドの互除法の定義より成り立つ．詳しくはdrken2019の2.ユークリッドの互除法の原理で)，以下のように$x,y$が求まる．
\begin{align*}
ax_0+by_0&=r_0\\
ax_0+by_0&=a\\
x_0&=1,\ y_0=0\\
ax_1+by_1&=r_1\\
ax_1+by_1&=b\\
x_1&=0,\ y_1=1
\end{align*}
$r_{i+1}=r_{i-1}-q_ir_i$ を $r_i=ax_i+by_i$を用いて表すと以下のようになる．
\begin{align*}
r_{i+1}&=(ax_{i-1}+by_{i-1})-q_i(ax_i+by_i)\\
&=a(x_{i-1}-q_ix_i)+b(y_{i-1}-q_iy_i)
\end{align*}
したがって
\begin{align*}
x_{i+1}&=x_{i-1}-q_ix_i\\
y_{i+1}&=y_{i-1}-q_iy_i
\end{align*}
となる．
漸化式なので行列で表すと以下のようになる．
\begin{align*}
Q_i&=\begin{pmatrix}
0&1\\
1&-q_i
\end{pmatrix}\\
R_i&=\begin{pmatrix}
r_{i-1}\\
r_{i}
\end{pmatrix},\
Z_i=\begin{pmatrix}
x_{i-1}&y_{i-1}\\
x_i&y_i
\end{pmatrix}\\
R_{i+1}&=Q_iR_i=Q_i\cdots Q_1R_1\\
Z_{i+1}&=Q_iZ_{i}=Q_i\cdots Q_1Z_1
\end{align*}
またこのときの初期状態$R_1, Z_1$は以下のようになる．
\begin{align*}
R_1&=\begin{pmatrix} a\\b\end{pmatrix}\\
Z_1&=\begin{pmatrix}1&0\\0&1\end{pmatrix}
\end{align*}
以下のプログラムでは，今までの数式に対して以下のように変数が対応している．
a, b は $r_{i-1},\ r_i$.
x, u は $x_{i-1},\ x_i$.
y, v は $y_{i-1},\ y_i$.
#include <type_traits>
template<class T>
constexpr T exgcd(T a, T b, T& x, T& y) {
auto assign = [&](T& s, T& t, T u, T v) -> void {
s = u; t = v;
};
// 初期状態 i=0, i=1 の係数設定
x = 1; y = 0; // x_{i-1}, y_{i-1}
T u = 0, v = 1; // x_i, y_i
// b (r_i) が 0 になるまで前進的に計算
while (b) {
T k = a / b; // q_i
T r = a % b; // r_{i+1}
// 状態遷移行列の乗算に相当
assign(a, b, b, r); // 剰余の更新(普通のユークリッドの互除法)
assign(x, u, u, x - u * k); // 係数 x の更新
assign(y, v, v, y - v * k); // 係数 y の更新
}
// 符号の正規化（GCDは正整数として返す）
if constexpr (std::is_signed_v<T>) {
if (a < 0) {
a = -a;
x = -x;
y = -y;
}
}
return a;
}
ここで，最後にx, y，つまり$x_{i-1},\ y_{i-1}$ を参照しているので，正しくないような気がするかもしれない．
ただ，$ax_i+by_i=r_i$ なので， $r_i=0$ のとき正しい値ではない．
$r_i=0$ になるとき， $r_{i-1}=\gcd(a,b)$ になるので，$x_{i-1}, y_{i-1}$ を用いるのが正しい．
また，これを用いて逆元を求める場合は，以下のようなプログラムになる．
template<class T>
constexpr T mod_inv(T a,T m){
T x,y,g=exgcd(a,m,x,y);
if(g!=1)return -1;
return(x%m+m)%m;
}
g!=1のとき，最大公約数が $1$ ではないので，法と値が互いに素ではないので逆元は存在しない．
C++では $0$ の方向へ丸められるので，xが負の数になったときは x%m+m で正しい値になる．
そして，正の数のときこの値は欲しい値よりm大きいので%mでそれを消す．
条件分岐や?:を使えばコンパイラはcmovなどの効率的な命令を吐き出すはずなので，除算よりはそちらのほうがいいだろう．
また，mが定数なら定数除算最適化で多少高速化が期待できるのでNTTPを積極的に使った方がいい．