Polynomial Gcd
(fps/polynomial-gcd.hpp)

• View this file on GitHub
• Last update: 2025-11-06 12:30:44+09:00
• Include: #include "fps/polynomial-gcd.hpp"

多項式 $f(x),g(x)$ に対する $\gcd(f(x),g(x))$ を $n=\max(\deg f,\deg g)$ として $O(n(\log n)^2)$ 時間で計算する．

また $f(x)h(x)\equiv 1\pmod{g(x)}$ となる多項式 $h(x)$ も同じ計算量で求められる．

Half-GCD 法

Euclid の互除法をそのまま適用すると，一回の反復で次数が $1$ しか減少しない場合に $\Omega(n^2\log n)$ 時間かかる． 解決するのが Half-GCD 法．

• [https://codeforces.com/blog/entry/101850]
• [https://dl.acm.org/doi/10.1145/800125.804045]

$f,g$ の GCD は互除法で求められる．

• $f_0=f,f_1=g$ とする（$\deg f\geq\deg g$）．
• $i=1,2,\dots$ について順に $f_{i-1}$ を $f_i$ で割った商を $q_i$ とし，$f_{i+1}=f_i-f_{i+1}q_i$ と定めていく．
• $f_{i+1}=0$ となったとき $\operatorname{GCD}(f,g)=f_i$ である．

過程は行列で表示できる：$M_i=\begin{pmatrix}0&1\1&-q_i\end{pmatrix}$ とすると $\begin{pmatrix}f_i\f_{i+1}\end{pmatrix}=M_i\cdots M_1\begin{pmatrix}f_0\f_1\end{pmatrix}$ である．

$\deg f_t \gt \frac{1}{2}\deg f\geq \deg f_{t+1}$ となる $t$ について $\operatorname{HGCD}(f,g)=M_t\cdots M_1$ とおく．

HGCD は再帰的に計算できる（証明は後述）．ただし $\deg f\geq\deg g$ とし，f >> k は $f\operatorname{div}x^k$ を意味する．

def HGCD(f, g):
  if deg(g) <= deg(f) / 2: return I
  k = floor(deg(f) / 2) + 1
  M = HGCD(f >> k, g >> k)
  (f, g)^T = M * (f, g)^T
  if deg(g) <= deg(f) / 2: return M
  (通常の互除法の 1 step を行う)
  if deg(g) <= deg(f) / 2: return M
  l = 2 * k - deg(f)
  return HGCD(f >> l, g >> l) * M

2 回の HGCD ではどちらも次数が半分以下になっている． また多項式除算を 1 回行っている． 従って $T(n)=2T(n/2)+n\log n$ を解いて $T(n)=n(\log n)^2$ であるから計算量は $O(n(\log n)^2)$ 時間である．

HGCD を用いれば以下のアルゴリズムで GCD が求められる． ループごとに $f$ の次数が半減するので $T(n)+T(n/2)+T(n/4)+\cdots=O(n(\log n)^2)$ 時間．

def GCD(f, g):
  while g != 0:
    (f, g) = HGCD(f, g) * (f, g)
    if g == 0: break
    (f, g) = (g, f % g)
  return f

正当性

形式的 Laurent 級数（有限個の負冪の項を許す）として考えると，$f$ を $g$ で割った商 $q$ は

\[q(x^{-1})=\frac{f(x^{-1})}{g(x^{-1})}\bmod x\]

と表せる．

補題. 多項式 $f,g\ (\deg f-\deg g=k)$ に対し $f$ を $g$ で割った商 $q$ は $f,g$ のそれぞれ上から $k+1$ 次の係数にのみ依存する（特に $\deg f$ の値によらない）．

証明

両辺に $x^k=x^{\deg f-\deg g}$ をかけると $$\begin{align*} x^kq(x^{-1}) &=\frac{x^{\deg f}f(x^{-1})}{x^{\deg g}g(x^{-1})}\bmod x^{k+1}\\ &=\frac{x^{\deg f}f(x^{-1})\bmod x^{k+1}}{x^{\deg g}g(x^{-1})\bmod x^{k+1}}\bmod x^{k+1} \end{align*}$$ であるのでよい． (証明終)

---

定理. $(f,g)$ に対し互除法を適用して得られる行列を $M_1,M_2,\dots$ とする．非負整数 $k$ に対し，ある $t$ が存在して $\operatorname{HGCD}(f\operatorname{div}x^k,g\operatorname{div}x^k)=M_t\cdots M_1$ となる．

証明

$n=\deg f$ とする． HGCD 側は $\deg f'_i\gt\frac{1}{2}(n-k)$ の間操作を行う． これをみたす $i$ に対して以下を示すことができればよい． - $\deg f_i=\deg f'_i+k$ - $f_i$ と $f'_i$ のそれぞれ上 $\deg f_{i-1}-\lfloor\frac{n+k}{2}\rfloor+1=\deg f'_{i-1}-\lfloor\frac{n-k}{2}\rfloor+1$ 次が一致する (TODO: 証明する，実験では成り立ちそう & これより次数評価は緩められなさそう) (証明終)

Depends on

• fps/formal-power-series.hpp

Verified with

• verify/fps/LC_inv_of_polynomials.test.cpp

Code

#pragma once
#include "fps/formal-power-series.hpp"

namespace polynomial_gcd {
template <class mint>
using FPS = FormalPowerSeries<mint>;
template <class mint>
using P = pair<FPS<mint>, FPS<mint>>;
template <class mint>
struct Mat {
  using fps = FPS<mint>;
  fps a00, a01, a10, a11;
  Mat() = default;
  Mat(const fps& a00_, const fps& a01_, const fps& a10_, const fps& a11_) : a00(a00_), a01(a01_), a10(a10_), a11(a11_) {}
  Mat& operator*=(const Mat& r) {
    fps b00 = a00 * r.a00 + a01 * r.a10;
    fps b01 = a00 * r.a01 + a01 * r.a11;
    fps b10 = a10 * r.a00 + a11 * r.a10;
    fps b11 = a10 * r.a01 + a11 * r.a11;
    swap(a00, b00);
    swap(a01, b01);
    swap(a10, b10);
    swap(a11, b11);
    a00.shrink();
    a01.shrink();
    a10.shrink();
    a11.shrink();
    return *this;
  }
  static Mat I() { return Mat(fps{mint(1)}, fps(), fps(), fps{mint(1)}); }
  Mat operator*(const Mat& r) const { return Mat(*this) *= r; }
};
template <typename mint>
P<mint> operator*(const Mat<mint>& m, const P<mint>& a) {
  using fps = FPS<mint>;
  fps b0 = m.a00 * a.first + m.a01 * a.second;
  fps b1 = m.a10 * a.first + m.a11 * a.second;
  b0.shrink();
  b1.shrink();
  return {b0, b1};
};

template <class mint>
void StepGCD(Mat<mint>& mat, P<mint>& p) {
  auto q = p.first / p.second;
  mat.a00 -= mat.a10 * q, mat.a00.shrink();
  mat.a01 -= mat.a11 * q, mat.a01.shrink();
  swap(mat.a00, mat.a10);
  swap(mat.a01, mat.a11);
  p.first -= p.second * q, p.first.shrink();
  swap(p.first, p.second);
}
template <class mint>
Mat<mint> HalfGCD(P<mint> p) {
  int n = (int)p.first.size() - 1, m = (int)p.second.size() - 1;
  int k = n / 2 + 1;
  if (m < k) return Mat<mint>::I();
  auto mat = HalfGCD(P<mint>{p.first >> k, p.second >> k});
  p = mat * p;
  if (p.second.size() <= k) return mat;
  StepGCD(mat, p);
  if (p.second.size() <= k) return mat;
  int l = 2 * k - ((int)p.first.size() - 1);
  return HalfGCD(P<mint>{p.first >> l, p.second >> l}) * mat;
}
template <class mint>
Mat<mint> ExtGCD(const FPS<mint>& f, const FPS<mint>& g) {
  P<mint> p{f, g};
  p.first.shrink();
  p.second.shrink();
  if (p.first.size() < p.second.size()) {
    auto mat = ExtGCD(p.second, p.first);
    swap(mat.a00, mat.a01);
    swap(mat.a10, mat.a11);
    return mat;
  }
  auto res = Mat<mint>::I();
  if (p.second.empty()) return res;
  while (true) {
    Mat<mint> mat = HalfGCD(p);
    p = mat * p;
    if (p.second.empty()) return mat * res;
    StepGCD(mat, p);
    if (p.second.empty()) return mat * res;
    res = mat * res;
  }
}

template <class mint>
FPS<mint> PolynomialGcd(FPS<mint> f, FPS<mint> g) {
  auto m = ExtGCD(f, g);
  auto h = m.a00 * f + m.a01 * g;
  h.shrink();
  if (!h.empty()) {
    mint c = h.back().inv();
    for (auto& v : h) v *= c;
  }
  return h;
}
// f*h=1 (mod g)
template <class mint>
optional<FPS<mint>> PolynomialInv(FPS<mint> f, FPS<mint> g) {
  auto m = ExtGCD(f, g);
  auto gcd = m.a00 * f + m.a01 * g;
  gcd.shrink();
  if (gcd.size() != 1) return nullopt;
  mint c = gcd[0].inv();
  auto h = m.a00;
  for (auto& v : h) v *= c;
  return h;
}
};  // namespace polynomial_gcd
using polynomial_gcd::PolynomialGcd;
using polynomial_gcd::PolynomialInv;
/**
 * @brief Polynomial Gcd
 * @docs docs/fps/polynomial-gcd.md
 */

#line 2 "fps/formal-power-series.hpp"

template <class mint>
struct FormalPowerSeries : vector<mint> {
  using vector<mint>::vector;
  using FPS = FormalPowerSeries;
  FormalPowerSeries(const vector<mint>& r) : vector<mint>(r) {}
  FormalPowerSeries(vector<mint>&& r) : vector<mint>(std::move(r)) {}
  FPS& operator=(const vector<mint>& r) {
    vector<mint>::operator=(r);
    return *this;
  }
  FPS& operator+=(const FPS& r) {
    if (r.size() > this->size()) this->resize(r.size());
    for (int i = 0; i < (int)r.size(); i++) (*this)[i] += r[i];
    return *this;
  }
  FPS& operator+=(const mint& r) {
    if (this->empty()) this->resize(1);
    (*this)[0] += r;
    return *this;
  }
  FPS& operator-=(const FPS& r) {
    if (r.size() > this->size()) this->resize(r.size());
    for (int i = 0; i < (int)r.size(); i++) (*this)[i] -= r[i];
    return *this;
  }
  FPS& operator-=(const mint& r) {
    if (this->empty()) this->resize(1);
    (*this)[0] -= r;
    return *this;
  }
  FPS& operator*=(const mint& v) {
    for (int k = 0; k < (int)this->size(); k++) (*this)[k] *= v;
    return *this;
  }
  FPS& operator/=(const FPS& r) {
    if (this->size() < r.size()) {
      this->clear();
      return *this;
    }
    int n = this->size() - r.size() + 1;
    if ((int)r.size() <= 64) {
      FPS f(*this), g(r);
      g.shrink();
      mint coeff = g.at(g.size() - 1).inv();
      for (auto& x : g) x *= coeff;
      int deg = (int)f.size() - (int)g.size() + 1;
      int gs = g.size();
      FPS quo(deg);
      for (int i = deg - 1; i >= 0; i--) {
        quo[i] = f[i + gs - 1];
        for (int j = 0; j < gs; j++) f[i + j] -= quo[i] * g[j];
      }
      *this = quo * coeff;
      this->resize(n, mint(0));
      return *this;
    }
    return *this = ((*this).rev().pre(n) * r.rev().inv(n)).pre(n).rev();
  }
  FPS& operator%=(const FPS& r) {
    *this -= *this / r * r;
    shrink();
    return *this;
  }
  FPS operator+(const FPS& r) const { return FPS(*this) += r; }
  FPS operator+(const mint& v) const { return FPS(*this) += v; }
  FPS operator-(const FPS& r) const { return FPS(*this) -= r; }
  FPS operator-(const mint& v) const { return FPS(*this) -= v; }
  FPS operator*(const FPS& r) const { return FPS(*this) *= r; }
  FPS operator*(const mint& v) const { return FPS(*this) *= v; }
  FPS operator/(const FPS& r) const { return FPS(*this) /= r; }
  FPS operator%(const FPS& r) const { return FPS(*this) %= r; }
  FPS operator-() const {
    FPS ret(this->size());
    for (int i = 0; i < (int)this->size(); i++) ret[i] = -(*this)[i];
    return ret;
  }
  void shrink() {
    while (this->size() && this->back() == mint(0)) this->pop_back();
  }
  FPS rev() const {
    FPS ret(*this);
    reverse(begin(ret), end(ret));
    return ret;
  }
  FPS dot(FPS r) const {
    FPS ret(min(this->size(), r.size()));
    for (int i = 0; i < (int)ret.size(); i++) ret[i] = (*this)[i] * r[i];
    return ret;
  }
  FPS pre(int sz) const {
    return FPS(begin(*this), begin(*this) + min((int)this->size(), sz));
  }
  FPS operator>>=(int sz) {
    assert(sz >= 0);
    if ((int)this->size() <= sz)
      this->clear();
    else
      this->erase(this->begin(), this->begin() + sz);
    return *this;
  }
  FPS operator>>(int sz) const {
    if ((int)this->size() <= sz) return {};
    FPS ret(*this);
    ret.erase(ret.begin(), ret.begin() + sz);
    return ret;
  }
  FPS operator<<=(int sz) {
    assert(sz >= 0);
    this->insert(this->begin(), sz, mint(0));
    return *this;
  }
  FPS operator<<(int sz) const {
    FPS ret(*this);
    ret.insert(ret.begin(), sz, mint(0));
    return ret;
  }
  FPS diff() const {
    const int n = (int)this->size();
    FPS ret(max(0, n - 1));
    mint one(1), coeff(1);
    for (int i = 1; i < n; i++) {
      ret[i - 1] = (*this)[i] * coeff;
      coeff += one;
    }
    return ret;
  }
  FPS integral() const {
    const int n = (int)this->size();
    FPS ret(n + 1);
    ret[0] = mint(0);
    if (n > 0) ret[1] = mint(1);
    auto mod = mint::get_mod();
    for (int i = 2; i <= n; i++) ret[i] = (-ret[mod % i]) * (mod / i);
    for (int i = 0; i < n; i++) ret[i + 1] *= (*this)[i];
    return ret;
  }
  mint eval(mint x) const {
    mint r = 0, w = 1;
    for (auto& v : *this) r += w * v, w *= x;
    return r;
  }
  FPS log(int deg = -1) const {
    assert((*this)[0] == mint(1));
    if (deg == -1) deg = (int)this->size();
    return (this->diff() * this->inv(deg)).pre(deg - 1).integral();
  }
  FPS pow(int64_t k, int deg = -1) const {
    const int n = (int)this->size();
    if (deg == -1) deg = n;
    if (k == 0) {
      FPS ret(deg);
      if (deg) ret[0] = 1;
      return ret;
    }
    for (int i = 0; i < n; i++) {
      if ((*this)[i] != mint(0)) {
        mint rev = mint(1) / (*this)[i];
        FPS ret = (((*this * rev) >> i).log(deg) * k).exp(deg);
        ret *= (*this)[i].pow(k);
        ret = (ret << (i * k)).pre(deg);
        if ((int)ret.size() < deg) ret.resize(deg, mint(0));
        return ret;
      }
      if (__int128_t(i + 1) * k >= deg) return FPS(deg, mint(0));
    }
    return FPS(deg, mint(0));
  }

  static void* ntt_ptr;
  static void set_ntt();
  FPS& operator*=(const FPS& r);
  FPS middle_product(const FPS& r) const;
  void ntt();
  void intt();
  void ntt_doubling();
  static int ntt_root();
  FPS inv(int deg = -1) const;
  FPS exp(int deg = -1) const;
};
template <typename mint>
void* FormalPowerSeries<mint>::ntt_ptr = nullptr;
#line 3 "fps/polynomial-gcd.hpp"

namespace polynomial_gcd {
template <class mint>
using FPS = FormalPowerSeries<mint>;
template <class mint>
using P = pair<FPS<mint>, FPS<mint>>;
template <class mint>
struct Mat {
  using fps = FPS<mint>;
  fps a00, a01, a10, a11;
  Mat() = default;
  Mat(const fps& a00_, const fps& a01_, const fps& a10_, const fps& a11_) : a00(a00_), a01(a01_), a10(a10_), a11(a11_) {}
  Mat& operator*=(const Mat& r) {
    fps b00 = a00 * r.a00 + a01 * r.a10;
    fps b01 = a00 * r.a01 + a01 * r.a11;
    fps b10 = a10 * r.a00 + a11 * r.a10;
    fps b11 = a10 * r.a01 + a11 * r.a11;
    swap(a00, b00);
    swap(a01, b01);
    swap(a10, b10);
    swap(a11, b11);
    a00.shrink();
    a01.shrink();
    a10.shrink();
    a11.shrink();
    return *this;
  }
  static Mat I() { return Mat(fps{mint(1)}, fps(), fps(), fps{mint(1)}); }
  Mat operator*(const Mat& r) const { return Mat(*this) *= r; }
};
template <typename mint>
P<mint> operator*(const Mat<mint>& m, const P<mint>& a) {
  using fps = FPS<mint>;
  fps b0 = m.a00 * a.first + m.a01 * a.second;
  fps b1 = m.a10 * a.first + m.a11 * a.second;
  b0.shrink();
  b1.shrink();
  return {b0, b1};
};

template <class mint>
void StepGCD(Mat<mint>& mat, P<mint>& p) {
  auto q = p.first / p.second;
  mat.a00 -= mat.a10 * q, mat.a00.shrink();
  mat.a01 -= mat.a11 * q, mat.a01.shrink();
  swap(mat.a00, mat.a10);
  swap(mat.a01, mat.a11);
  p.first -= p.second * q, p.first.shrink();
  swap(p.first, p.second);
}
template <class mint>
Mat<mint> HalfGCD(P<mint> p) {
  int n = (int)p.first.size() - 1, m = (int)p.second.size() - 1;
  int k = n / 2 + 1;
  if (m < k) return Mat<mint>::I();
  auto mat = HalfGCD(P<mint>{p.first >> k, p.second >> k});
  p = mat * p;
  if (p.second.size() <= k) return mat;
  StepGCD(mat, p);
  if (p.second.size() <= k) return mat;
  int l = 2 * k - ((int)p.first.size() - 1);
  return HalfGCD(P<mint>{p.first >> l, p.second >> l}) * mat;
}
template <class mint>
Mat<mint> ExtGCD(const FPS<mint>& f, const FPS<mint>& g) {
  P<mint> p{f, g};
  p.first.shrink();
  p.second.shrink();
  if (p.first.size() < p.second.size()) {
    auto mat = ExtGCD(p.second, p.first);
    swap(mat.a00, mat.a01);
    swap(mat.a10, mat.a11);
    return mat;
  }
  auto res = Mat<mint>::I();
  if (p.second.empty()) return res;
  while (true) {
    Mat<mint> mat = HalfGCD(p);
    p = mat * p;
    if (p.second.empty()) return mat * res;
    StepGCD(mat, p);
    if (p.second.empty()) return mat * res;
    res = mat * res;
  }
}

template <class mint>
FPS<mint> PolynomialGcd(FPS<mint> f, FPS<mint> g) {
  auto m = ExtGCD(f, g);
  auto h = m.a00 * f + m.a01 * g;
  h.shrink();
  if (!h.empty()) {
    mint c = h.back().inv();
    for (auto& v : h) v *= c;
  }
  return h;
}
// f*h=1 (mod g)
template <class mint>
optional<FPS<mint>> PolynomialInv(FPS<mint> f, FPS<mint> g) {
  auto m = ExtGCD(f, g);
  auto gcd = m.a00 * f + m.a01 * g;
  gcd.shrink();
  if (gcd.size() != 1) return nullopt;
  mint c = gcd[0].inv();
  auto h = m.a00;
  for (auto& v : h) v *= c;
  return h;
}
};  // namespace polynomial_gcd
using polynomial_gcd::PolynomialGcd;
using polynomial_gcd::PolynomialInv;
/**
 * @brief Polynomial Gcd
 * @docs docs/fps/polynomial-gcd.md
 */

Back to top page