拡張ユークリッドの互除法 (extgcd)
(library/math/extgcd.hpp)

View this file on GitHub
Last update: 2025-03-03 00:33:48+09:00
Include: #include "library/math/extgcd.hpp"

ユークリッドの互除法なので最大公約数も求められるとよい気はするが、まあ $\gcd(a, b)$ なんて stl に存在するので、最大公約数無視して 1 次不定方程式を解くやつを extgcd として実装する。

解きたい問題

非負整数 $a, b, c$ が与えられたとき、 1 次不定方程式 $ax + by = c$ となる整数解 $(x, y)$ の 1 つを求める。

解法

ユークリッドの互除法同様に処理を行うことで、

\[\begin{eqnarray*} a &=& b q_0 + r_0 \\ b &=& r_0 q_1 + r_1 \\ r_0 &=& r_1 q_2 + r_2 \\ &\vdots& \\ r_{n-3} &=& r_{n-2} q_{n-1} + r_{n-1} \\ r_{n-2} &=& r_{n-1} q_{n} + r_{n} \\ r_{n-1} &=& r_{n} q_{n+1} + 0 \\ \end{eqnarray*}\]

となる $r_n$ が求まる。このとき、$r_n$ は $a$ と $b$ の最大公約数である。

この操作を逆方向に行うことによって、

\[\begin{eqnarray*} r_{n} &=& r_{n-2} - r_{n-1} q_{n} \\ &=& r_{n-2} - (r_{n-3} - r_{n-2} q_{n-1}) q_{n} \\ &\vdots& \\ &=& x a + y b \end{eqnarray*}\]

となるような整数 $x, y$ の 1 つが求まる。

$c$ が $\gcd(a, b)$ の倍数ではない場合、整数解は存在しない。一方、$c$ が $\gcd(a, b)$ の倍数である場合、上で求まった $x, y$ を $\dfrac{c}{\gcd(a, b)}$ 倍すればよい。

実装

再帰関数で実装できる。

「$a, b, c$ を受け取って $ax + by = c$ となる $x, y$ を返す関数 $f$」としておく。 $a = bq + r$ から $b = r q’ + r’$ に遷移することを考えると、「$bx’ + ry’ = c$」つまり $f(b, r, c) = f(b, a \bmod b, c)$ を求めることを繰り返していく。最終的に $b = 0$ になったら $ax + by = c = ax$ となるから、 $(x, y) = (c / a, 0)$ を返してやればよい。

そこから逆算していく。ある状態 $ax + by = \gcd(a, b)$ において、 $bx’ + ry’ = \gcd(a, b)$ が再帰で求まっている。このとき、 $r = a - bq$ を使うと、 $ax + by = bx’ + ry’ = bx’ + (a - bq)y’ = ay’ + b(x’ - qy’)$ となるから、 $(x, y) = (y’, x’ - qy’)$ を返してやればよい。

使い方

$\text{extgcd}(a, b, c)$: $ax + by = c$ となる整数解 $(x, y)$ の 1 つを optional<pair<T, T>> で返す。存在しない場合は nullopt を返す。計算量は $O(\log \min(a, b))$ 。

Required by

Modint (library/data-structure/modint.hpp)

Verified with

Code

#pragma once

#include <concepts>
#include <optional>
#include <utility>
using namespace std;

namespace asalib {
  namespace math {
    // Returns a pair (x, y) such that ax + by = c
    template<integral T>
    constexpr optional<pair<T, T>> extgcd(T a, T b, T c) {
      if (b == 0) {
        if (c % a != 0) return nullopt;
        return make_pair(c / a, 0);
      }
      auto res = extgcd(b, a % b, c);
      if (!res) return nullopt;
      auto [x, y] = *res;
      return make_pair(y, x - (a / b) * y);
    }
  }  // namespace math
}  // namespace asalib

#line 2 "library/math/extgcd.hpp"

#include <concepts>
#include <optional>
#include <utility>
using namespace std;

namespace asalib {
  namespace math {
    // Returns a pair (x, y) such that ax + by = c
    template<integral T>
    constexpr optional<pair<T, T>> extgcd(T a, T b, T c) {
      if (b == 0) {
        if (c % a != 0) return nullopt;
        return make_pair(c / a, 0);
      }
      auto res = extgcd(b, a % b, c);
      if (!res) return nullopt;
      auto [x, y] = *res;
      return make_pair(y, x - (a / b) * y);
    }
  }  // namespace math
}  // namespace asalib

拡張ユークリッドの互除法 (extgcd) (library/math/extgcd.hpp)

解きたい問題

解法

実装

使い方

Required by

Verified with

Code

拡張ユークリッドの互除法 (extgcd)
(library/math/extgcd.hpp)