疎 FPS の次数圧縮・前処理

components/fps_sparse_core.hpp · GitHub で見る

共通ページ

fps_detail::SparseBase<Engine> は、疎 FPS の次数の最大公約数による圧縮・小さな整数冪の計算を行う内部部品です。通常版の疎演算 から使います。項の型・収集処理は fps_sparse_terms.hpp で軽量版と共有します。

Terms・Kind・collect・get の引数・返り値・前提条件は通常版と同じです。Engine は solve<kind>・geometric<exponential>・scale を提供します。get は次数を圧縮して演算を委譲し、結果を元の次数へ戻します。SIMD 演算には依存しません。

非零係数数を \(k\)、求める係数数を \(n\) として、演算を含む時間は \(O(nk)\)、追加空間は \(O(n+k)\) です。

依存ライブラリ

• 疎 FPS の項収集・演算種別
• modint998244353

検証コード

• test/convolution.test.cpp
• test/fps_bostan_mori.test.cpp
• test/fps_compose.test.cpp
• test/fps_compose_large.test.cpp
• test/fps_compositional_inverse.test.cpp
• test/fps_compositional_inverse_large.test.cpp
• test/fps_divmod.test.cpp
• test/fps_exp.test.cpp
• test/fps_exp_sparse.test.cpp
• test/fps_interpolate.test.cpp
• test/fps_inv.test.cpp
• test/fps_inv_sparse.test.cpp
• test/fps_linear_recurrence.test.cpp
• test/fps_log.test.cpp
• test/fps_log_sparse.test.cpp
• test/fps_pow.test.cpp
• test/fps_pow_sparse.test.cpp
• test/fps_rational_coefficients.test.cpp
• test/fps_sqrt.test.cpp
• test/fps_sqrt_sparse.test.cpp
• test/fps_taylor_shift.test.cpp
• test/multipoint_eval.test.cpp
• test/polynomial_inv_mod.test.cpp
• test/power_projection.test.cpp

コード

#ifndef FPS_SPARSE_CORE_HPP
#define FPS_SPARSE_CORE_HPP

#include "fps_sparse_terms.hpp"
#include <algorithm>
#include <cassert>
#include <numeric>
#include <utility>
#include <vector>

namespace fps_detail {
template <class E> struct SparseBase : SparseTerms {
	static Terms multiply(const Terms &a, const Terms &b, int n) {
		Terms result;
		for (auto [i, x] : a) for (auto [j, y] : b) if (i < n - j) result.emplace_back(i + j, x * y);
		std::sort(result.begin(), result.end(), [](auto a, auto b) { return a.first < b.first; });
		std::size_t used = 0;
		for (auto [i, x] : result) {
			if (used && result[used - 1].first == i) result[used - 1].second += x;
			else result[used++] = {i, x};
		}
		result.resize(used);
		result.erase(std::remove_if(result.begin(), result.end(), [](auto term) { return term.second == mint(0); }), result.end());
		return result;
	}
	template <Kind kind>
	static Series get(Terms terms, int n, mint initial, mint exponent, const mint *(*inverses)(int)) {
		assert(n >= 0 && terms.size() <= 64);
		if constexpr (kind != Kind::inverse) assert(unsigned(n) < mint::mod());
		if (!n) return {};
		if (terms.empty()) { Series out(n); out[0] = initial; return out; }
		if constexpr (kind == Kind::power) {
			unsigned e = exponent.val();
			if (e <= 4 && terms.size() <= 12 && (n >= 4096 || terms.size() <= 3)) {
				Terms power{{0, 1}}, f{{0, 1}};
				f.insert(f.end(), terms.begin(), terms.end());
				while (e) {
					if (e & 1) power = multiply(power, f, n);
					e >>= 1;
					if (e) f = multiply(f, f, n);
				}
				Series out(n);
				for (auto [i, x] : power) out[i] = x * initial;
				return out;
			}
		}
		int stride = 0;
		for (auto term : terms) stride = std::gcd(stride, term.first);
		const int count = (n - 1) / stride + 1;
		for (auto &term : terms) term.first /= stride;
		const mint *inv = nullptr;
		if constexpr (kind != Kind::inverse) inv = inverses(count);
		Series compressed;
		if (terms.size() == 1 && kind != Kind::power) {
			compressed.resize(count);
			if constexpr (kind == Kind::exponential) {
				E::template geometric<true>(compressed, initial, terms[0].second, inv);
			} else {
				E::geometric(compressed, kind == Kind::logarithm ? mint(-1) : initial, -terms[0].second);
				if constexpr (kind == Kind::logarithm) compressed[0] = 0, E::scale(compressed, inv);
			}
		} else compressed = E::template solve<kind>(terms, count, initial, exponent, inv);
		if (stride == 1) return compressed;
		Series out(n);
		for (int i = 0; i < count; ++i) out[i * stride] = compressed[i];
		return out;
	}
};
}
#endif

Back to top page