"""
    extended_euclidean(a, b)

The function extended_euclidean takes two integers,
a and b, as input and returns a tuple (g, x, y), where g is
the greatest common divisor of a and b, and x and y are the
Bezout's coefficients such that ax + by = g (i.e a*x = g mod b)

The algorithm uses recursion to compute the gcd and the
coefficients. The base case is when b is zero, in which case the gcd
is a and the coefficients are 1 and 0. Otherwise, it recursively calls
itself with b and a % b as the new inputs, and computes the gcd and
coefficients based on the results of the recursive call.
"""
function extended_euclidean(a::Integer, b::Integer)
    if b == 0
        return (a, 1, 0)
    else
        g, x, y = extended_euclidean(b, a % b)
        # return (g, y, x - (a ÷ b) * y)
        return (g, y, x - div(a, b) * y)
    end
end

"""
    modular_inverse(e, phi)

Given integer e and phi, gcd(e, phi) == 1, find such d that d*e = 1 mod phi
"""
function modular_inverse(e::Integer, phi::Integer)
    g, x, _ = extended_euclidean(e, phi)
    g == 1 || error("$e and $phi have the greatest common divisor > 1")
    if x < 0
        x += phi
    end
    return x
end