Extended Euclidean Algorithm in modern and readable C++

Question

In Python the Extended Euclidean Algorithm (egcd) could be written as follows:

def egcd(a, b):
  if b == 0:
    return (a, 1, 0)
  else:
    (d, tmp, s) = egcd(b, a%b)
    return (d, s, tmp - (a//b) * s)

I want to write a native and modern C++ version of the egcd. This is what I have so far:

template<typename A, typename B, typename T = std::common_type_t<A, B>>
constexpr std::tuple<T, T, T> egcd(A a, B b) {
  static_assert(std::is_integral<A>::value, "arguments to egcd are integers");
  static_assert(std::is_integral<B>::value, "arguments to egcd are integers");

  if (b == 0) {
    return std::make_tuple(a, 1, 0);
  } else {
    auto triple = egcd(b, a%b);
    return std::make_tuple(std::get<0>(triple), std::get<2>(triple), std::get<1>(triple) - (a/b) * std::get<2>(triple))
  }
}

As you can see using <tuple> makes the code very long and hard (compare with the Python version..) to read.

Any suggestions to make this code more readable?

Answer 1

The code is rather tiny so there isn't much to say, but I still have a couple of remarks:

• Since you're using C++17, you can take advantage of variable templates to simplify your static assertions a bit:

  static_assert(std::is_integral_v<A>, "arguments to egcd are integers");
  static_assert(std::is_integral_v<B>, "arguments to egcd are integers");
  

• The LWG issues 2733 and 2759 highlight the fact that any possibly cv-qualified bool satisfies the std::is_integral trait, but the GCD and LCM operations don't make much sense for boolean types. Consequently, the extended operations don't make much sense for boolean types eithers. You could add the following static assertions to make it clear:

  static_assert(not std::is_same_v<std::decay_t<A>, bool>, "arguments to egcd are not booleans");
  static_assert(not std::is_same_v<std::decay_t<B>, bool>, "arguments to egcd are not booleans");
  

• Since you're always computing both a / b and a % b, you can compute std::div instead which computes both at once and returns a structure with quot and rem members. Actually, the algorithm used to compute the remainder of a and b also computes the quotient of a and b, so functions like std::div exist in several programming languages so that users can take advantage of it instead of discarding a computed value and computing it again.

• They're not widely available at the moment, but if you intend to keep returning an std::tuple, a C++17 decomposition declaration would probably make your code more readable by providing explicit names instead of the opaque std::get<position>:

  if (b == 0) {
      return std::make_tuple(a, 1, 0);
  } else {
      auto divmod = std::div(a, b);
      auto [gcd, bezout_x, bezout_y] = egcd(b, divmod.rem);
      return std::make_tuple(gcd, bezout_y, bezout_x - divmod.quot * bezout_y)
  }
  

Answer 2

So, if my brief look at Wikipedia is correct, the algorithm produces a "Bézout's identity", which happens to be two numbers. EDIT: and the gcd.

Don't represent this as a tuple. Create a simple object with two integers. It's easy to write, offers more type-safety and is more expressive.

I'll write some pseudo code to show you what I mean.

func (a int, b int) BagOStuff {
    ...[Complex stuff goes here]
    return BagOStuff( item1: foo, item2: bar, item3: zim);
}

Not great. Even though my fantasy tuple's are tenser than C++ they still suck. I've missed the opportunity to use the object I'm returning to express my intent, and, without digging into the [complex stuff] I can't even remember what number was meant to do what.

func (a int, b int) EgcdResult {
    ...[Complex stuff goes here]
    return EgcdResult( gcd: foo, benzoutX: bar, benzoutY: zim);
}

type EgcdResult struct {
    gcd int
    benzoutX int
    benzoutY int
}

This is much better. As a caller I know exactly what to do with this.

Answer 3

Firstly, you didn't provide a main() so I had to improvise:

#include <iostream>
int main()
{
    auto result = egcd(29664648L, 8974651687326L);
    std::cout << std::get<0>(result) << std::endl;
}

I also had to add

#include <tuple>
#include <type_traits>

and add a missing semicolon to make your code syntactically valid. In future, please compile exactly the code you post.

---

If you're willing to live on the bleeding edge, then Concepts can replace your static_assert() lines:

constexpr std::tuple<T, T, T> egcd(A a, B b)
    requires std::is_integral<A>::value
          && std::is_integral<B>::value
{

If you have (and enable) Library Fundamentals V1, you can replace std::is_integral<X>::value with std::experimental::is_integral_v<X>.

---

To reduce the clumsiness of all those std::get calls, it would be nice to be able to std::tie the result into variables like this:

    T gcd, rx, ry;
    std::tie(gcd, rx, ry) = egcd(b, a%b);
    return std::make_tuple(gcd, rx, rx - (a/b) * ry);

but that's not allowed in a constexpr function.

However, since the elements are all of the same type, you may be better off with a std::array, giving:

#include <array>
#include <type_traits>

template<typename A, typename B, typename T = std::common_type_t<A, B>>
constexpr std::array<T, 3> egcd(A a, B b)
    requires std::is_integral<A>::value && std::is_integral<B>::value
{
    if (b == 0) {
        return { a, 1, 0 };
    } else {
        const auto triple = egcd(b, a%b);
        return { triple[0], triple[1], triple[1] - (a/b) * triple[2] };
    }
}

We didn't need to change main() because the signature of std::get(std::array) exactly parallels that of std::get(std::tuple).