I couldn't understand the non-recursive implementation of the extended Euclidean algorithm used in ABC186-E, so I'll summarize what I researched and understood in this article. The extended Euclidean algorithm itself is omitted in detail because there were many easy-to-understand articles such as Extended Euclidean algorithm ~ How to solve the linear indefinite equation ax + by = c ~. A special feature on how to find "too much divided by"! ~ From the inverse element to the discrete logarithm ~](https://qiita.com/drken/items/3b4fdf0a78e7a138cd9a#3-%E5%89%B2%E3%82%8A%E7%AE%97-a--b) summarizes the principle of non-recursive implementation of the extended Euclidean algorithm.
Suddenly, let's express the extended Euclidean algorithm in the form of a matrix.
ax + by = \gcd(a,b)
When solving
a = \lfloor\frac{a}{b}\rfloor b + r
As, the problem of $ (a, b) $ is replaced with the problem of $ (b, r) $ and solved. Expressing this as a matrix,
\begin{bmatrix}
a \\
b
\end{bmatrix}
=
\begin{bmatrix}
k & 1 \\
1 & 0
\end{bmatrix}
\begin{bmatrix}
b \\
r
\end{bmatrix}
It will be. ($ \ Lfloor \ frac {a} {b} \ rfloor = k $). The extended Euclidean algorithm repeats this until it becomes $ (\ gcd (a, b), 0) $. Therefore, $ a = r_0, b = r_1, r = r_2, \ dots, \ gcd (a, b) = r_h $, and $ k_i = \ lfloor \ frac {r_ {i-1}} {r_i} \ rfloor If you write $, you can express it as follows.
\begin{bmatrix}
r_0 \\
r_1
\end{bmatrix}
=
\begin{bmatrix}
k_0 & 1 \\
1 & 0
\end{bmatrix}
\dots
\begin{bmatrix}
k_{h-1} & 1 \\
1 & 0
\end{bmatrix}
\begin{bmatrix}
r_h \\
0
\end{bmatrix} \tag{1}
here,
\begin{bmatrix}
k_i & 1 \\
1 & 0
\end{bmatrix}
The inverse matrix of
\begin{bmatrix}
0 & 1 \\
1 & -k_i
\end{bmatrix}
Therefore, the $ (1) $ formula can be transformed as follows.
\begin{bmatrix}
r_h \\
0
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 \\
1 & -k_{h-1}
\end{bmatrix}
\dots
\begin{bmatrix}
0 & 1 \\
1 & -k_0
\end{bmatrix}
\begin{bmatrix}
r_0 \\
r_1
\end{bmatrix} \tag{2}
Finally,
\begin{bmatrix}
x & y \\
u & v
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 \\
1 & -k_{h-1}
\end{bmatrix}
\dots
\begin{bmatrix}
0 & 1 \\
1 & -k_0
\end{bmatrix} \tag{3}
If you put, the $ (2) $ formula will be
\begin{bmatrix}
r_h \\
0
\end{bmatrix}
=
\begin{bmatrix}
x & y \\
u & v
\end{bmatrix}
\begin{bmatrix}
r_0 \\
r_1
\end{bmatrix}
This $ (x, y) $ becomes the solution $ (x, y) $ of the opening $ ax + by = \ gcd (a, b) $.
Consider an algorithm that actually finds $ (x, y) $ from the above matrix representation. Since it can be calculated based on the $ (3) $ formula, set the initial value to $ (x_0 = 1, y_0 = 0, u_0 = 0, v_0 = 1) $ (identity matrix).
\begin{bmatrix}
x_1 & y_1 \\
u_1 & v_1
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 \\
1 & -k_0
\end{bmatrix}
\begin{bmatrix}
x_0 & y_0 \\
u_0 & v_0
\end{bmatrix}
When you calculate
\begin{array}{llll}
x_1 = u_0 \\
y_1 = v_0 \\
u_1 = x_0 - k_0u_0 \\
v_1 = y_0 - k_0v_0
\end{array}
As a program,
x -= k[0]*u
y -= k[0]*v
x,u = u,x
y,v = v,y
It will be an image like. All you have to do is repeat this,
def extGCD(a,b):
x, y, u, v = 1, 0, 0, 1
while b:
k = a // b
x -= k * u
y -= k * v
x, u = u, x
y, v = v, y
a, b = b, a % b
return x, y
In this way, you can implement the extended Euclidean algorithm in a non-recursive way.
Like ABC186-E written at the beginning, I think that the extended Euclidean algorithm is often used to find the modular reciprocal in the case of competitive pros. When $ (a, b) $ are relatively prime, the solution $ x $ of $ ax + by = 1 $ becomes $ ax \ equiv 1 (\ mod b) $, so the extended Euclidean algorithm is used to calculate the modular reciprocal. You can use the reciprocal method.
In this case, the value of $ (v, y) $ in the last program of 2. is not particularly required, so it can be calculated by the following program. (Some post-loop processing has been added to make $ x $ a positive value.)
def modinv(a, b):
p = b
x, y, u, v = 1, 0, 0, 1
while b:
k = a // b
x -= k * u
y -= k * v
x, u = u, x
y, v = v, y
a, b = b, a % b
x %= p
if x < 0:
x += p
return x
Recommended Posts