多項式 (polynomial)
(include/emthrm/math/polynomial.hpp)

• View this file on GitHub
• Last update: 2023-02-23 21:59:12+09:00
• Include: #include "emthrm/math/polynomial.hpp"

時間計算量

	時間計算量
加減算	$O(N)$
スカラー倍	$O(N)$
乗算	$O(N\log{N})$
除算	$O(N^2)$
剰余演算	$O(N^2)$
ホーナー法 (Horner’s rule)	$O(N)$
形式微分	$O(N)$
累乗
translate	$O(N\log{N})$

仕様

template <typename T>
struct Polynomial;

メンバ変数

名前	説明
std::vector<T> coef	係数列

メンバ関数

名前	効果・戻り値
explicit Polynomial(const int deg = 0);	$\mathrm{deg}$ 次まで係数列をもつオブジェクトを構築する。
explicit Polynomial(const std::vector<T>& coef);	係数列を $\mathrm{coef}$ とするオブジェクトを構築する。
Polynomial(const std::initializer_list<T> init);	初期化子リストを受け取るコンストラクタ
template <typename InputIter> explicit Polynomial(const InputIter first, const InputIter last);	イテレータ範囲コンストラクタ
inline const T& operator[](const int term) const; inline T& operator[](const int term);	${\lbrack x^{\mathrm{term}} \rbrack}f$
void resize(const int deg);	$\mathrm{deg}$ 次まで係数列をもつ。
void shrink();	正規化を行う。
int degree() const;	現在の次数
Polynomial& operator=(const std::vector<T>& coef_); Polynomial& operator=(const Polynomial& x);	代入演算子
Polynomial& operator+=(const Polynomial& x); Polynomial operator+(const Polynomial& x) const;	加算
Polynomial& operator-=(const Polynomial& x); Polynomial operator-(const Polynomial& x) const;	減算
Polynomial& operator*=(const T x); Polynomial& operator*=(const Polynomial& x); Polynomial operator*(const T x) const; Polynomial operator*(const Polynomial& x) const;	乗算
Polynomial& operator/=(const T x); Polynomial& operator/=(const Polynomial& x); Polynomial operator/(const T x) const; Polynomial operator/(const Polynomial& x) const;	除算
Polynomial& operator%=(const Polynomial& x); Polynomial operator%(const Polynomial& x) const;	剰余演算
std::pair<Polynomial, Polynomial> divide(Polynomial x) const;	$x$ で割った商とあまり
Polynomial& operator<<=(const int n); FormalPowerSeries operator<<(const int n) const;	$x^n f$
bool operator==(Polynomial x) const;	等値比較
Polynomial operator+() const;	$+{f}$
Polynomial operator-() const;	$-{f}$
T horner(const T x) const;	$f(x)$
Polynomial differential() const;	$f^{\prime}$
Polynomial pow(int exponent) const;	$f^{\mathrm{exponent}}$
Polynomial translate(const T c) const;	$f(x + c)$

静的メンバ関数

名前	効果
static void set_mult(const Mult mult);	乗算を定義する。

メンバ型

名前	説明
Mult	std::function<std::vector<T>(const std::vector<T>&, const std::vector<T>&)>

参考文献

• https://github.com/beet-aizu/library/blob/29e15f77befa18d06b1f61221c509b5a58f4cb4c/polynomial/polynomial.cpp
• https://github.com/primenumber/ProconLib/blob/dc175741e461e88f69350532b47fa57527b38fa3/Math/Polynomial.cpp

TODO

• http://www.math.kobe-u.ac.jp/Asir/ca.pdf
• 除算の高速化
  • https://github.com/spaghetti-source/algorithm/blob/master/math/polynomial_int.cc
  • http://people.csail.mit.edu/madhu/ST12/scribe/lect06.pdf
  • https://maspypy.com/%e5%a4%9a%e9%a0%85%e5%bc%8f%e3%83%bb%e5%bd%a2%e5%bc%8f%e7%9a%84%e3%81%b9%e3%81%8d%e7%b4%9a%e6%95%b0-%e9%ab%98%e9%80%9f%e3%81%ab%e8%a8%88%e7%ae%97%e3%81%a7%e3%81%8d%e3%82%8b%e3%82%82%e3%81%ae#toc10
  • https://judge.yosupo.jp/problem/division_of_polynomials
• DKA 法 (Durand–Kerner–Aberth method)
  • https://drive.google.com/file/d/1Fb5TPXzEQVKzYDwhNN2XZpWZWa49aCIY
• Bairstow’s method
  • https://en.wikipedia.org/wiki/Bairstow%27s_method
• 微積分
  • https://topcoder-g-hatena-ne-jp.jag-icpc.org/pepsin-amylase/20141208/

Verified with

• 数学/多項式 (test/math/polynomial.test.cpp)

Code

#ifndef EMTHRM_MATH_POLYNOMIAL_HPP_
#define EMTHRM_MATH_POLYNOMIAL_HPP_

#include <algorithm>
#include <cassert>
#include <functional>
#include <initializer_list>
#include <iterator>
#include <numeric>
#include <utility>
#include <vector>

namespace emthrm {

template <typename T>
struct Polynomial {
  std::vector<T> coef;

  explicit Polynomial(const int deg = 0) : coef(deg + 1, 0) {}
  explicit Polynomial(const std::vector<T>& coef) : coef(coef) {}
  Polynomial(const std::initializer_list<T> init)
      : coef(init.begin(), init.end()) {}
  template <typename InputIter>
  explicit Polynomial(const InputIter first, const InputIter last)
      : coef(first, last) {}

  inline const T& operator[](const int term) const { return coef[term]; }
  inline T& operator[](const int term) { return coef[term]; }

  using Mult = std::function<std::vector<T>(const std::vector<T>&,
                                            const std::vector<T>&)>;
  static void set_mult(const Mult mult) { get_mult() = mult; }

  void resize(const int deg) { coef.resize(deg + 1, 0); }
  void shrink() {
    while (coef.size() > 1 && coef.back() == 0) coef.pop_back();
  }
  int degree() const { return std::ssize(coef) - 1; }

  Polynomial& operator=(const std::vector<T>& coef_) {
    coef = coef_;
    return *this;
  }
  Polynomial& operator=(const Polynomial& x) = default;

  Polynomial& operator+=(const Polynomial& x) {
    const int deg_x = x.degree();
    if (deg_x > degree()) resize(deg_x);
    for (int i = 0; i <= deg_x; ++i) {
      coef[i] += x[i];
    }
    return *this;
  }
  Polynomial& operator-=(const Polynomial& x) {
    const int deg_x = x.degree();
    if (deg_x > degree()) resize(deg_x);
    for (int i = 0; i <= deg_x; ++i) {
      coef[i] -= x[i];
    }
    return *this;
  }
  Polynomial& operator*=(const T x) {
    for (T& e : coef) e *= x;
    return *this;
  }
  Polynomial& operator*=(const Polynomial& x) {
    return *this = get_mult()(coef, x.coef);
  }
  Polynomial& operator/=(const T x) {
    assert(x != 0);
    return *this *= static_cast<T>(1) / x;
  }
  std::pair<Polynomial, Polynomial> divide(Polynomial x) const {
    x.shrink();
    Polynomial rem = *this;
    const int n = rem.degree(), m = x.degree(), deg = n - m;
    if (deg < 0) return {Polynomial{0}, rem};
    Polynomial quo(deg);
    for (int i = 0; i <= deg; ++i) {
      quo[deg - i] = rem[n - i] / x[m];
      for (int j = 0; j <= m; ++j) {
        rem[n - i - j] -= x[m - j] * quo[deg - i];
      }
    }
    rem.resize(deg);
    return {quo, rem};
  }
  Polynomial& operator/=(const Polynomial& x) {
    return *this = divide(x).first;
  }
  Polynomial& operator%=(const Polynomial& x) {
    return *this = divide(x).second;
  }
  Polynomial& operator<<=(const int n) {
    coef.insert(coef.begin(), n, 0);
    return *this;
  }

  bool operator==(Polynomial x) const {
    x.shrink();
    Polynomial y = *this;
    y.shrink();
    return x.coef == y.coef;
  }

  Polynomial operator+() const { return *this; }
  Polynomial operator-() const {
    Polynomial res = *this;
    for (T& e : res.coef) e = -e;
    return res;
  }

  Polynomial operator+(const Polynomial& x) const {
    return Polynomial(*this) += x;
  }
  Polynomial operator-(const Polynomial& x) const {
    return Polynomial(*this) -= x;
  }
  Polynomial operator*(const T x) const {
    return Polynomial(*this) *= x;
  }
  Polynomial operator*(const Polynomial& x) const {
    return Polynomial(*this) *= x;
  }
  Polynomial operator/(const T x) const { return Polynomial(*this) /= x; }
  Polynomial operator/(const Polynomial& x) const {
    return Polynomial(*this) /= x;
  }
  Polynomial operator%(const Polynomial& x) const {
    return Polynomial(*this) %= x;
  }
  Polynomial operator<<(const int n) const { return Polynomial(*this) <<= n; }

  T horner(const T x) const {
    return std::accumulate(
        coef.rbegin(), coef.rend(), static_cast<T>(0),
        [x](const T l, const T r) -> T { return l * x + r; });
  }

  Polynomial differential() const {
    const int deg = degree();
    assert(deg >= 0);
    Polynomial res(std::max(deg - 1, 0));
    for (int i = 1; i <= deg; ++i) {
      res[i - 1] = coef[i] * i;
    }
    return res;
  }

  Polynomial pow(int exponent) const {
    Polynomial res{1}, base = *this;
    for (; exponent > 0; exponent >>= 1) {
      if (exponent & 1) res *= base;
      base *= base;
    }
    return res;
  }

  Polynomial translate(const T c) const {
    const int n = coef.size();
    std::vector<T> fact(n, 1), inv_fact(n, 1);
    for (int i = 1; i < n; ++i) {
      fact[i] = fact[i - 1] * i;
    }
    inv_fact[n - 1] = static_cast<T>(1) / fact[n - 1];
    for (int i = n - 1; i > 0; --i) {
      inv_fact[i - 1] = inv_fact[i] * i;
    }
    std::vector<T> g(n), ex(n);
    for (int i = 0; i < n; ++i) {
      g[i] = coef[i] * fact[i];
    }
    std::reverse(g.begin(), g.end());
    T pow_c = 1;
    for (int i = 0; i < n; ++i) {
      ex[i] = pow_c * inv_fact[i];
      pow_c *= c;
    }
    const std::vector<T> conv = get_mult()(g, ex);
    Polynomial res(n - 1);
    for (int i = 0; i < n; ++i) {
      res[i] = conv[n - 1 - i] * inv_fact[i];
    }
    return res;
  }

 private:
  static Mult& get_mult() {
    static Mult mult = [](const std::vector<T>& a, const std::vector<T>& b)
        -> std::vector<T> {
      const int n = a.size(), m = b.size();
      std::vector<T> res(n + m - 1, 0);
      for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
          res[i + j] += a[i] * b[j];
        }
      }
      return res;
    };
    return mult;
  }
};

}  // namespace emthrm

#endif  // EMTHRM_MATH_POLYNOMIAL_HPP_

#line 1 "include/emthrm/math/polynomial.hpp"



#include <algorithm>
#include <cassert>
#include <functional>
#include <initializer_list>
#include <iterator>
#include <numeric>
#include <utility>
#include <vector>

namespace emthrm {

template <typename T>
struct Polynomial {
  std::vector<T> coef;

  explicit Polynomial(const int deg = 0) : coef(deg + 1, 0) {}
  explicit Polynomial(const std::vector<T>& coef) : coef(coef) {}
  Polynomial(const std::initializer_list<T> init)
      : coef(init.begin(), init.end()) {}
  template <typename InputIter>
  explicit Polynomial(const InputIter first, const InputIter last)
      : coef(first, last) {}

  inline const T& operator[](const int term) const { return coef[term]; }
  inline T& operator[](const int term) { return coef[term]; }

  using Mult = std::function<std::vector<T>(const std::vector<T>&,
                                            const std::vector<T>&)>;
  static void set_mult(const Mult mult) { get_mult() = mult; }

  void resize(const int deg) { coef.resize(deg + 1, 0); }
  void shrink() {
    while (coef.size() > 1 && coef.back() == 0) coef.pop_back();
  }
  int degree() const { return std::ssize(coef) - 1; }

  Polynomial& operator=(const std::vector<T>& coef_) {
    coef = coef_;
    return *this;
  }
  Polynomial& operator=(const Polynomial& x) = default;

  Polynomial& operator+=(const Polynomial& x) {
    const int deg_x = x.degree();
    if (deg_x > degree()) resize(deg_x);
    for (int i = 0; i <= deg_x; ++i) {
      coef[i] += x[i];
    }
    return *this;
  }
  Polynomial& operator-=(const Polynomial& x) {
    const int deg_x = x.degree();
    if (deg_x > degree()) resize(deg_x);
    for (int i = 0; i <= deg_x; ++i) {
      coef[i] -= x[i];
    }
    return *this;
  }
  Polynomial& operator*=(const T x) {
    for (T& e : coef) e *= x;
    return *this;
  }
  Polynomial& operator*=(const Polynomial& x) {
    return *this = get_mult()(coef, x.coef);
  }
  Polynomial& operator/=(const T x) {
    assert(x != 0);
    return *this *= static_cast<T>(1) / x;
  }
  std::pair<Polynomial, Polynomial> divide(Polynomial x) const {
    x.shrink();
    Polynomial rem = *this;
    const int n = rem.degree(), m = x.degree(), deg = n - m;
    if (deg < 0) return {Polynomial{0}, rem};
    Polynomial quo(deg);
    for (int i = 0; i <= deg; ++i) {
      quo[deg - i] = rem[n - i] / x[m];
      for (int j = 0; j <= m; ++j) {
        rem[n - i - j] -= x[m - j] * quo[deg - i];
      }
    }
    rem.resize(deg);
    return {quo, rem};
  }
  Polynomial& operator/=(const Polynomial& x) {
    return *this = divide(x).first;
  }
  Polynomial& operator%=(const Polynomial& x) {
    return *this = divide(x).second;
  }
  Polynomial& operator<<=(const int n) {
    coef.insert(coef.begin(), n, 0);
    return *this;
  }

  bool operator==(Polynomial x) const {
    x.shrink();
    Polynomial y = *this;
    y.shrink();
    return x.coef == y.coef;
  }

  Polynomial operator+() const { return *this; }
  Polynomial operator-() const {
    Polynomial res = *this;
    for (T& e : res.coef) e = -e;
    return res;
  }

  Polynomial operator+(const Polynomial& x) const {
    return Polynomial(*this) += x;
  }
  Polynomial operator-(const Polynomial& x) const {
    return Polynomial(*this) -= x;
  }
  Polynomial operator*(const T x) const {
    return Polynomial(*this) *= x;
  }
  Polynomial operator*(const Polynomial& x) const {
    return Polynomial(*this) *= x;
  }
  Polynomial operator/(const T x) const { return Polynomial(*this) /= x; }
  Polynomial operator/(const Polynomial& x) const {
    return Polynomial(*this) /= x;
  }
  Polynomial operator%(const Polynomial& x) const {
    return Polynomial(*this) %= x;
  }
  Polynomial operator<<(const int n) const { return Polynomial(*this) <<= n; }

  T horner(const T x) const {
    return std::accumulate(
        coef.rbegin(), coef.rend(), static_cast<T>(0),
        [x](const T l, const T r) -> T { return l * x + r; });
  }

  Polynomial differential() const {
    const int deg = degree();
    assert(deg >= 0);
    Polynomial res(std::max(deg - 1, 0));
    for (int i = 1; i <= deg; ++i) {
      res[i - 1] = coef[i] * i;
    }
    return res;
  }

  Polynomial pow(int exponent) const {
    Polynomial res{1}, base = *this;
    for (; exponent > 0; exponent >>= 1) {
      if (exponent & 1) res *= base;
      base *= base;
    }
    return res;
  }

  Polynomial translate(const T c) const {
    const int n = coef.size();
    std::vector<T> fact(n, 1), inv_fact(n, 1);
    for (int i = 1; i < n; ++i) {
      fact[i] = fact[i - 1] * i;
    }
    inv_fact[n - 1] = static_cast<T>(1) / fact[n - 1];
    for (int i = n - 1; i > 0; --i) {
      inv_fact[i - 1] = inv_fact[i] * i;
    }
    std::vector<T> g(n), ex(n);
    for (int i = 0; i < n; ++i) {
      g[i] = coef[i] * fact[i];
    }
    std::reverse(g.begin(), g.end());
    T pow_c = 1;
    for (int i = 0; i < n; ++i) {
      ex[i] = pow_c * inv_fact[i];
      pow_c *= c;
    }
    const std::vector<T> conv = get_mult()(g, ex);
    Polynomial res(n - 1);
    for (int i = 0; i < n; ++i) {
      res[i] = conv[n - 1 - i] * inv_fact[i];
    }
    return res;
  }

 private:
  static Mult& get_mult() {
    static Mult mult = [](const std::vector<T>& a, const std::vector<T>& b)
        -> std::vector<T> {
      const int n = a.size(), m = b.size();
      std::vector<T> res(n + m - 1, 0);
      for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
          res[i + j] += a[i] * b[j];
        }
      }
      return res;
    };
    return mult;
  }
};

}  // namespace emthrm

Back to top page