多項式の共通積木

components/polynomial_product_tree.hpp · GitHub で見る

共通ページ

polynomial_detail::ProductTree<F> は、FPS の型 F を受け取る多点評価・補間の内部部品です。高速 FPS の積木と軽量 FPS の積木が共有します。通常は各版の多点評価・補間の公開関数を使ってください。

各葉に \(x-x_i\) を置き、親には子の多項式の積を保存します。評価では子へ多項式の剰余を渡し、補間では積多項式の導関数を評価して各点の重みを求め、子の結果を結合します。積・除算・微分などは F へ委譲します。

内部 API	内容・返り値
ProductTree<F>(xs)	評価点から構築。T = F::value_type とし、xs は vector<T>
tree.evaluate(f)	点列と同じ順序の評価値を vector<T> で返す。点の重複を許す
tree.interpolate(ys)	点列に対応する値から、点数と同じ係数数の F を返す
tree.count / tree.size	点数／点数以上の最小の 2 冪。空列の size は 1
tree.product	1 始まりの節点配列。添字 0 は未使用

product[1] は \(\displaystyle\prod_{i=0}^{m-1}(x-x_i)\) です。空の点列では定数 1 を保持し、評価・補間の結果は空列です。補間では点数と値の個数が一致し、各点における積多項式の導関数の値が可逆である必要があります。体の係数なら点が互いに異なることと同値です。構築後の product・size・count は変更しないでください。

点数は \(2^{23}\) 以下です。点数が上限と一致するときは、最上位の積を分割して計算し、畳み込みを上限内に収めます。ProductTree<F, Limit> で畳み込みの上限を別に指定することもできます。評価する多項式の長さと、剰余計算などに必要な長さも F の上限に収めてください。係数型と畳み込みの条件は使う FPS に従います。

点数を \(m\)、評価する多項式の係数数を \(n\) として、構築・補間は \(O(m\log^2 m)\) 時間、木の保持空間は \(O(m\log m)\) です。評価は \(O((n+m)\log^2(n+m))\) 時間・\(O(n+m)\) 追加空間、補間の追加空間は \(O(m)\) です。追加空間に構築済みの木は含めません。

依存ライブラリ

• 多項式の積・打ち切り積の共通処理

検証コード

• test/fps_interpolate.test.cpp
• test/fps_light_compose.test.cpp
• test/fps_light_compositional_inverse.test.cpp
• test/fps_light_convolution.test.cpp
• test/fps_light_divmod.test.cpp
• test/fps_light_exp.test.cpp
• test/fps_light_interpolate.test.cpp
• test/fps_light_inv.test.cpp
• test/fps_light_linear_recurrence.test.cpp
• test/fps_light_log.test.cpp
• test/fps_light_multipoint_eval.test.cpp
• test/fps_light_polynomial_extended_gcd.test.cpp
• test/fps_light_polynomial_gcd.test.cpp
• test/fps_light_polynomial_inv_mod.test.cpp
• test/fps_light_pow.test.cpp
• test/fps_light_power_projection.test.cpp
• test/fps_light_rational_coefficients.test.cpp
• test/fps_light_sqrt.test.cpp
• test/fps_light_taylor_shift.test.cpp
• test/multipoint_eval.test.cpp
• test/online_convolution_light.test.cpp

コード

#ifndef CP_COMPONENTS_POLYNOMIAL_PRODUCT_TREE_HPP
#define CP_COMPONENTS_POLYNOMIAL_PRODUCT_TREE_HPP

#include "polynomial_multiply.hpp"
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <utility>
#include <vector>

namespace polynomial_detail {
template<class F, std::size_t Limit = 1U << 23>
struct ProductTree : Products<F, Limit> {
	using Products<F, Limit>::multiply;
	using Products<F, Limit>::multiply_low;
	static_assert(Limit >= 2);
	using T = typename F::value_type;
	std::size_t size = 1, count;
	std::vector<F> product;
	static void trim(F& f) { while (!f.empty() && f.back() == T(0)) f.pop_back(); }


	explicit ProductTree(const std::vector<T>& xs) : count(xs.size()) {
		assert(count <= Limit);
		while (size < count) size <<= 1;
		product.resize(2 * size);
		for (std::size_t i = 0; i < size; ++i)
			product[size + i] = i < count ? F{-xs[i], T(1)} : F{T(1)};
		for (std::size_t i = size; --i;) {
			const auto& a = product[2 * i];
			const auto& b = product[2 * i + 1];
			product[i] = b.size() == 1 ? a : multiply(a, b);
		}
	}

	std::vector<T> evaluate(const F& f) const {
		std::vector<T> result(count);
		if (!count) return result;
		auto remainder = [](const F& a, const F& b) {
			return a.size() < b.size() ? a : a.divmod(b).second;
		};
		auto visit = [&](auto&& self, std::size_t node, const F& value) -> void {
			if (value.empty()) return;
			if (node >= size) {
				result[node - size] = value[0];
				return;
			}
			self(self, node * 2, remainder(value, product[node * 2]));
			if (product[node * 2 + 1].size() > 1)
				self(self, node * 2 + 1, remainder(value, product[node * 2 + 1]));
		};
		F value = f;
		trim(value);
		visit(visit, 1, remainder(value, product[1]));
		return result;
	}

	F interpolate(const std::vector<T>& ys) const {
		assert(ys.size() == count);
		if (!count) return {};
		F derivative = product[1].diff();
		trim(derivative);
		std::vector<T> denominator = evaluate(derivative), prefix(count + 1, T(1));
		for (std::size_t i = 0; i < count; ++i) {
			assert(denominator[i] != T(0));
			prefix[i + 1] = prefix[i] * denominator[i];
		}
		T inverse = prefix[count].inv();
		std::vector<F> value(2 * size);
		for (std::size_t i = count; i--;) {
			value[size + i] = F{ys[i] * inverse * prefix[i]};
			inverse *= denominator[i];
		}
		for (std::size_t i = size; --i;) {
			value[i] = F::multiply(value[2 * i], product[2 * i + 1])
				+ F::multiply(value[2 * i + 1], product[2 * i]);
			value[2 * i] = F();
			value[2 * i + 1] = F();
		}
		value[1].resize(count);
		return std::move(value[1]);
	}
};
}

#endif

Back to top page