Bézout's identity (tools/bezout.hpp)

template <typename T>
std::optional<std::tuple<T, T, T, T>> bezout(T a, T b, T c);

It solves $ax + by = c$ for $x$ and $y$ over the integers. It returns four integers $p, q, r, s$ if the solutions of the equation exist. They mean that $x$ can be denoted as $pm + q$ and $y$ can be denoted as $rm + s$ for any integer $m$.

Constraints

$a \neq 0$
$b \neq 0$
$ac$ is in <T>.
$bc$ is in <T>.

Time Complexity

$O(\log(\min(|a|, |b|)))$

Note

The solutions of the equation exist if and only if $c \equiv 0 \pmod{\gcd(a, b)}$. If the solutions exist, the following equations hold when we denote a particular solution of $a x’ + b y’ = \gcd(a, b)$ as $(x’_0, y’_0)$.

\[\begin{align*} p &= -\frac{b}{\gcd(a, b)}\\ q &= \frac{c}{\gcd(a, b)} x'_0\\ r &= \frac{a}{\gcd(a, b)}\\ s &= \frac{c}{\gcd(a, b)} y'_0 \end{align*}\]

If $a, b, c, x, y \geq 0$, $\left\lceil -\frac{c y’_0}{a} \right\rceil \leq m \leq \left\lfloor \frac{c x’_0}{b} \right\rfloor$ holds.

License

Author

anqooqie

Depends on

Verified with

tests/bezout.test.cpp

Code

#ifndef TOOLS_BEZOUT_HPP
#define TOOLS_BEZOUT_HPP

#include <optional>
#include <tuple>
#include "tools/extgcd.hpp"

namespace tools {
  template <typename T>
  ::std::optional<::std::tuple<T, T, T, T>> bezout(const T a, const T b, const T c) {
    assert(a != T(0));
    assert(b != T(0));
    const auto [x0, y0, gcd] = ::tools::extgcd(a, b);
    return c % gcd == T(0) ? ::std::make_optional(::std::make_tuple(-b / gcd, c / gcd * x0, a / gcd, c / gcd * y0)) : ::std::nullopt;
  }
}

#endif

#line 1 "tools/bezout.hpp"



#include <optional>
#include <tuple>
#line 1 "tools/extgcd.hpp"



#line 5 "tools/extgcd.hpp"
#include <utility>
#include <cassert>
#include <algorithm>
#line 1 "tools/abs.hpp"



namespace tools {
  constexpr float abs(const float x) {
    return x < 0 ? -x : x;
  }
  constexpr double abs(const double x) {
    return x < 0 ? -x : x;
  }
  constexpr long double abs(const long double x) {
    return x < 0 ? -x : x;
  }
  constexpr int abs(const int x) {
    return x < 0 ? -x : x;
  }
  constexpr long abs(const long x) {
    return x < 0 ? -x : x;
  }
  constexpr long long abs(const long long x) {
    return x < 0 ? -x : x;
  }
  constexpr unsigned int abs(const unsigned int x) {
    return x;
  }
  constexpr unsigned long abs(const unsigned long x) {
    return x;
  }
  constexpr unsigned long long abs(const unsigned long long x) {
    return x;
  }
}


#line 9 "tools/extgcd.hpp"

namespace tools {

  template <typename T>
  ::std::tuple<T, T, T> extgcd(T prev_r, T r) {
    const bool prev_r_is_neg = prev_r < T(0);
    const bool r_is_neg = r < T(0);

    if (prev_r_is_neg) prev_r = -prev_r;
    if (r_is_neg) r = -r;

    #ifndef NDEBUG
    const T a = prev_r;
    const T b = r;
    #endif

    T prev_s(1);
    T prev_t(0);
    T s(0);
    T t(1);
    while (r != T(0)) {
      const T q = prev_r / r;
      ::std::tie(prev_r, r) = ::std::make_pair(r, prev_r - q * r);
      ::std::tie(prev_s, s) = ::std::make_pair(s, prev_s - q * s);
      ::std::tie(prev_t, t) = ::std::make_pair(t, prev_t - q * t);
    }

    if (prev_r_is_neg) prev_s = -prev_s;
    if (r_is_neg) prev_t = -prev_t;

    assert(::tools::abs(prev_s) <= ::std::max(b / prev_r / T(2), T(1)));
    assert(::tools::abs(prev_t) <= ::std::max(a / prev_r / T(2), T(1)));
    return ::std::make_tuple(prev_s, prev_t, prev_r);
  }
}


#line 7 "tools/bezout.hpp"

namespace tools {
  template <typename T>
  ::std::optional<::std::tuple<T, T, T, T>> bezout(const T a, const T b, const T c) {
    assert(a != T(0));
    assert(b != T(0));
    const auto [x0, y0, gcd] = ::tools::extgcd(a, b);
    return c % gcd == T(0) ? ::std::make_optional(::std::make_tuple(-b / gcd, c / gcd * x0, a / gcd, c / gcd * y0)) : ::std::nullopt;
  }
}