Convert To Newton Basis
(library/polynomial/convert_to_newton_basis.hpp)

View this file on GitHub
Last update: 2023-05-21 01:49:53+09:00
Include: #include "library/polynomial/convert_to_newton_basis.hpp"

Convert To Newton Basis

多項式 $\displaystyle f(x) = \sum _ {i = 0} ^ {n - 1} a _ i x ^ i$ および長さ $n$ の列 $\lbrace p _ i \rbrace _ {i = 0} ^ {n - 1}$ に対して、次を満たす列 $\lbrace b _ i \rbrace _ {i = 0} ^ {n - 1}$ を $O(n (\log n) ^ 2)$ 時間で計算するライブラリ。

\[f(x) = \sum _ {i = 0} ^ {n - 1} b _ i \prod _ {j = 0} ^ {i - 1} (x - p _ j).\]

アルゴリズム

分割統治法による。$0\leq l \lt r\leq n$ なる整数組 $(l, r)$ に対して $\displaystyle f _ {l, r}(x) \coloneqq \sum _ {i = l} ^ {r - 1} b _ i \prod _ {j = l} ^ {i - 1} (x - p _ j)$ と定義すると、整数 $m\in\lbrack l, r)$ に対して次が成り立つ。

\[\begin{aligned} f _ {l, r}(x) &= \sum _ {i = l} ^ {m - 1} b _ i \prod _ {j = l} ^ {i - 1} (x - p _ j) + \left(\prod _ {i = l} ^ {m - 1} (x - p _ i)\right) \sum _ {i = m} ^ {r - 1} b _ i \prod _ {j = m} ^ {i - 1} (x - p _ j) \newline &= f _ {l, m}(x) + g _ {l, m}(x) f _ {m, r}(x). \quad \left(g _ {l, m}(x) \coloneqq \prod _ {i = l} ^ {m - 1} (x - p _ i)\right) \end{aligned}\]

ここで $\deg g _ {l, m} = m - l,\ \deg f _ {l, m} = m - l - 1$ に注目すると、$f _ {m, r}(x)$ と $f _ {l, m}(x)$ はそれぞれ $f _ {l, r}(x)$ を $g _ {l, m}(x)$ で割った商と余りに対応することが分かる。

従って、二分木状に $g _ {l, r}$ をボトムアップに計算し (subproduct tree)、トップダウンに $f _ {l, r}$ を計算することで $O(n (\log n) ^ 2)$ 時間で全ての $i = 0, 1, \ldots, n - 1$ に対する $f _ {i, i + 1}(x) \equiv b _ i$ を計算することができる。

Required by

library/polynomial/rook_polynomial.hpp

Verified with

Code

#ifndef SUISEN_CONVERT_TO_NEWTON_BASIS
#define SUISEN_CONVERT_TO_NEWTON_BASIS

#include <tuple>
#include <vector>

namespace suisen {
    // Returns b=(b_0,...,b_{N-1}) s.t. f(x) = Sum[i=0,N-1] b_i Prod[j=0,i-1](x - p_j)
    template <typename FPSType>
    std::vector<typename FPSType::value_type> convert_to_newton_basis(const FPSType& f, const std::vector<typename FPSType::value_type>& p) {
        const int n = p.size();
        assert(f.size() == n);

        int m = 1;
        while (m < n) m <<= 1;

        std::vector<FPSType> seg(2 * m);
        for (int i = 0; i < m; ++i) {
            seg[m + i] = { i < n ? -p[i] : 0, 1 };
        }
        for (int i = m - 1; i > 0; --i) {
            if (((i + 1) & -(i + 1)) == (i + 1)) continue; // i = 2^k - 1
            seg[i] = seg[2 * i] * seg[2 * i + 1];
        }

        seg[1] = f;
        for (int i = 1; i < m; ++i) {
            std::tie(seg[2 * i + 1], seg[2 * i]) = seg[i].div_mod(seg[2 * i]);
        }

        std::vector<typename FPSType::value_type> b(n);
        for (int i = 0; i < n; ++i) {
            b[i] = seg[m + i].safe_get(0);
        }
        return b;
    }
} // namespace suisen


#endif // SUISEN_CONVERT_TO_NEWTON_BASIS

#line 1 "library/polynomial/convert_to_newton_basis.hpp"



#include <tuple>
#include <vector>

namespace suisen {
    // Returns b=(b_0,...,b_{N-1}) s.t. f(x) = Sum[i=0,N-1] b_i Prod[j=0,i-1](x - p_j)
    template <typename FPSType>
    std::vector<typename FPSType::value_type> convert_to_newton_basis(const FPSType& f, const std::vector<typename FPSType::value_type>& p) {
        const int n = p.size();
        assert(f.size() == n);

        int m = 1;
        while (m < n) m <<= 1;

        std::vector<FPSType> seg(2 * m);
        for (int i = 0; i < m; ++i) {
            seg[m + i] = { i < n ? -p[i] : 0, 1 };
        }
        for (int i = m - 1; i > 0; --i) {
            if (((i + 1) & -(i + 1)) == (i + 1)) continue; // i = 2^k - 1
            seg[i] = seg[2 * i] * seg[2 * i + 1];
        }

        seg[1] = f;
        for (int i = 1; i < m; ++i) {
            std::tie(seg[2 * i + 1], seg[2 * i]) = seg[i].div_mod(seg[2 * i]);
        }

        std::vector<typename FPSType::value_type> b(n);
        for (int i = 0; i < n; ++i) {
            b[i] = seg[m + i].safe_get(0);
        }
        return b;
    }
} // namespace suisen

Convert To Newton Basis (library/polynomial/convert_to_newton_basis.hpp)