軽量 FPS の疎演算

components/fps_sparse_scalar.hpp · GitHub で見る

共通ページ

fps_detail::ScalarSparseSeries は、非零係数が少ない FPS の逆数・対数・指数・冪を計算する内部部品です。軽量版 FPS から利用します。

Terms・Kind・collect・get の引数・返り値・前提条件は 通常版の疎演算 と同じです。係数型は mint です。

次数の最大公約数で係数列を圧縮し、指数の計算では次数と係数の積をあらかじめ求めます。非零係数より前の区間では漸化式の計算を省きます。

非零項が 1 つの冪は二項係数の漸化式で求めます。小さい非負の冪は項の組合せを直接列挙します。

非零係数数を \(k\)、求める係数数を \(n\) として、時間 \(O(nk)\)、追加空間 \(O(n+k)\) です。逆元表の準備時間・空間は別途必要です。

依存ライブラリ

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

検証コード

• 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/online_convolution_light.test.cpp

コード

#ifndef FPS_SPARSE_SCALAR_HPP
#define FPS_SPARSE_SCALAR_HPP

#include "fps_sparse_terms.hpp"

namespace fps_detail {
struct ScalarSparseSeries : SparseTerms {
	template <Kind kind>
	static Series get(Terms terms, int n, mint initial, mint exponent, const mint *(*inverses)(int)) {
		assert(n >= 0 && terms.size() <= 64 && (kind == Kind::inverse || unsigned(n) < mint::mod()));
		if (!n) return {};
		Series out(n);
		out[0] = initial;
		if (terms.empty()) 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, initial}};
				for (unsigned k = 0; k < e; ++k) {
					const auto previous = power;
					for (auto [i, x] : previous) for (auto [j, y] : terms)
						if (i < n - j) power.emplace_back(i + j, x * y);
				}
				out[0] = 0;
				for (auto [i, x] : power) out[i] += x;
				return out;
			}
		}
		int stride = 0;
		for (auto [j, x] : terms) stride = std::gcd(stride, j);
		const int count = (n - 1) / stride + 1;
		for (auto &[j, x] : terms) j /= stride;
		if constexpr (kind == Kind::exponential) for (auto &[j, x] : terms) x *= mint(j);
		const mint *inv = kind == Kind::inverse ? nullptr : inverses(count);
		if (terms.size() == 1) {
			mint value = kind == Kind::logarithm ? mint(-1) : initial;
			mint ratio = kind == Kind::exponential ? terms[0].second : -terms[0].second;
			for (int i = 1; i < count; ++i) {
				value *= ratio;
				if constexpr (kind == Kind::power) value *= (mint(i - 1) - exponent) * inv[i];
				if constexpr (kind == Kind::exponential) value *= inv[i];
				out[i * stride] = value;
			}
		} else for (int i = terms[0].first; i < count; ++i) {
			mint sum = 0;
			for (auto [j, x] : terms) {
				if (j > i) break;
				if constexpr (kind == Kind::power) x *= (exponent + mint(1)) * mint(j) - mint(i);
				if constexpr (kind == Kind::inverse || kind == Kind::logarithm) x = -x;
				sum += x * out[(i - j) * stride];
				if constexpr (kind == Kind::logarithm) if (j == i) sum -= mint(j) * x;
			}
			if constexpr (kind == Kind::power || kind == Kind::exponential) sum *= inv[i];
			out[i * stride] = sum;
		}
		if constexpr (kind == Kind::logarithm) for (int i = 1; i < count; ++i) out[i * stride] *= inv[i];
		return out;
	}
};
}
#endif

Back to top page