多項式補間 (polynomial interpolation)
(include/emthrm/math/formal_power_series/polynomial_interpolation.hpp)

• View this file on GitHub
• Last update: 2023-02-25 16:35:06+09:00
• Include: #include "emthrm/math/formal_power_series/polynomial_interpolation.hpp"

ラグランジュ補間 (Lagrange interpolation)

$1 \leq i < j \leq N,\ x_i \neq x_j$ を満たす $(x_i, y_i)$ に対して $f(x_i) = y_i$ を満たす $N - 1$ 次以下の多項式 $f$ を求める。

ラグランジュの補間多項式 (interpolation polynomial in the Lagrange form)

\[f(x) = \sum_{i = 1}^N f(x_i) \prod_{j \neq i} \dfrac{x - x_j}{x_i - x_j} = \sum_{i = 1}^N \dfrac{f(x_i)}{g^{\prime}(x_i)} \prod_{j \neq i} (x - x_j) \text{ where } g(x) = \prod_{i = 1}^N (x - x_i).\]

時間計算量

	時間計算量
評価版	$O(N^2)$
評価版2	$O(N)$
多項式補間	$O(N(\log{N})^2)$

仕様

評価版

名前	戻り値	要件
template <typename T> T lagrange_interpolation(const std::vector<T>& x, const std::vector<T>& y, const T t);	$f(x_i) = y_i$ を満たす $f(t)$	$x_i$ は互いに異なる。
template <typename T> T lagrange_interpolation(const std::vector<T>& y, const T t);	$f(i) = y_i$ を満たす $f(t)$	$t < 0,\ N \leq t$

多項式補間

名前	戻り値	備考
template <template <typename> class C, typename T> C<T> polynomial_interpolation(const std::vector<T>& x, const std::vector<T>& y);	$f(x_i) = y_i$ を満たす $f$	C は冪級数を表す構造体である。

参考文献

ラグランジュ補間

• https://mathtrain.jp/hokan
• https://tubo28.me/compprog/algorithm/lagrange_interpolation/

多項式補間

• https://github.com/ei1333/library/blob/c5a89ad7cb3855ebb448da00a629504e26adce24/math/fps/polynomial-interpolation.cpp

TODO

• スプライン補間 (spline interpolation)
  • https://en.wikipedia.org/wiki/Spline_interpolation
  • https://docs.google.com/presentation/d/1n-RR-AjQN8oq0DkDKBFi9Ccg53HsDUHKtOsC8z7S-To
• ニュートン補間 (Newtonian interpolation)
  • https://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%88%E3%83%B3%E8%A3%9C%E9%96%93
  • https://twitter.com/noshi91/status/1160191749416898560
• shift of sampling points of polynomial
  • https://ei1333.hateblo.jp/entry/2021/07/08/221742
  • https://hly1204.github.io/library/math/formal_power_series/sample_points_shift.hpp
  • https://judge.yosupo.jp/problem/shift_of_sampling_points_of_polynomial

Submissons

• 評価版
• 評価版2
• 多項式補間

Depends on

• multipoint evaluation (include/emthrm/math/formal_power_series/multipoint_evaluation.hpp)

Verified with

• 数学/形式的冪級数/多項式補間 (test/math/formal_power_series/polynomial_interpolation.test.cpp)

Code

#ifndef EMTHRM_MATH_FORMAL_POWER_SERIES_POLYNOMIAL_INTERPOLATION_HPP_
#define EMTHRM_MATH_FORMAL_POWER_SERIES_POLYNOMIAL_INTERPOLATION_HPP_

#include <vector>

#include "emthrm/math/formal_power_series/multipoint_evaluation.hpp"

namespace emthrm {

template <template <typename> class C, typename T>
C<T> polynomial_interpolation(const std::vector<T>& x,
                              const std::vector<T>& y) {
  MultipointEvaluation<C, T> m(x);
  m.build(m.subproduct_tree[1].differential());
  const int n = x.size();
  const auto f = [&y, &m, n](auto f, const int node) -> C<T> {
    return node >= n ? C<T>{y[node - n] / m.f_x[node - n]} :
                       f(f, node << 1) * m.subproduct_tree[(node << 1) + 1]
                       + f(f, (node << 1) + 1) * m.subproduct_tree[node << 1];
  };
  return f(f, 1);
}

}  // namespace emthrm

#endif  // EMTHRM_MATH_FORMAL_POWER_SERIES_POLYNOMIAL_INTERPOLATION_HPP_

#line 1 "include/emthrm/math/formal_power_series/polynomial_interpolation.hpp"



#include <vector>

#line 1 "include/emthrm/math/formal_power_series/multipoint_evaluation.hpp"



#include <algorithm>
#include <iterator>
#line 7 "include/emthrm/math/formal_power_series/multipoint_evaluation.hpp"

namespace emthrm {

template <template <typename> class C, typename T>
struct MultipointEvaluation {
  std::vector<T> f_x;
  std::vector<C<T>> subproduct_tree;

  explicit MultipointEvaluation(const std::vector<T> &xs)
      : f_x(xs.size()), subproduct_tree(xs.size() << 1), n(xs.size()) {
    std::transform(xs.begin(), xs.end(), std::next(subproduct_tree.begin(), n),
                   [](const T& x) -> C<T> { return C<T>{-x, 1}; });
    for (int i = n - 1; i > 0; --i) {
      subproduct_tree[i] =
          subproduct_tree[i << 1] * subproduct_tree[(i << 1) + 1];
    }
  }

  void build(const C<T>& f) { dfs(f, 1); }

 private:
  const int n;

  void dfs(C<T> f, int node) {
    f %= subproduct_tree[node];
    if (node < n) {
      dfs(f, node << 1);
      dfs(f, (node << 1) + 1);
    } else {
      f_x[node - n] = f[0];
    }
  }
};

}  // namespace emthrm


#line 7 "include/emthrm/math/formal_power_series/polynomial_interpolation.hpp"

namespace emthrm {

template <template <typename> class C, typename T>
C<T> polynomial_interpolation(const std::vector<T>& x,
                              const std::vector<T>& y) {
  MultipointEvaluation<C, T> m(x);
  m.build(m.subproduct_tree[1].differential());
  const int n = x.size();
  const auto f = [&y, &m, n](auto f, const int node) -> C<T> {
    return node >= n ? C<T>{y[node - n] / m.f_x[node - n]} :
                       f(f, node << 1) * m.subproduct_tree[(node << 1) + 1]
                       + f(f, (node << 1) + 1) * m.subproduct_tree[node << 1];
  };
  return f(f, 1);
}

}  // namespace emthrm

Back to top page