FPS 逆数のブロック変換

components/fps_inverse_blocked.hpp · GitHub で見る

共通ページ

fps::inv() が使う内部部品です。通常は公開 FPS API から利用します。Linux x86_64 / GCC 用に定義され、呼び出し側で AVX2 対応を実行時に確認します。他の環境では既存の FPS 実装を使います。

操作	内容・返り値	計算量
native_inverse(in,input,out,n)	逆数の先頭 n 係数を空の out へ追加	\(O(n\log n)\)
append_scaled(a,out,scale)	a の係数を scale 倍して out へ追加	予約済み: 追加係数数に比例、再確保時: 既存要素数も加算

関数は fps_detail::blocked_poly 内にあります。通常は FPS の f.inv(n) を使います。

fps_detail::blocked_poly::native_inverse(in,input,out,n) は、in の input 個の Mint 係数から、逆数 \(1/f\bmod x^n\) を求め、空の vector<Mint> である out に n 項を追加します。1 <= n <= 2^23、input >= 1、定数項が非零であることを前提とします。内部ワードの範囲は \([0,2p)\)、\(p=998244353\) です。入力の範囲外を 0 とし、次数 n 以降は読みません。入力は変更せず、out と重ねません。公開配列に padding やアラインメントは要求しません。

計算は共有 Transform と Newton 更新を使います。最大変換長の半分まで逆数を求めた後、最後の更新の作業領域に入力の内部コピーを再利用します。ここで \(R=2^{32}\bmod p\) とします。入力の内部ワードを通常の整数として扱うと \(Rf\) となるため、初期の高々 32 項の逆数を \(C=R^2\) 倍し、全更新で \(C/(Rf)=R/f\) の表現を保ちます。Newton 更新の 2 回目の積の倍率へ \(-C^{-1}=53925088\bmod p\) を組み込みます。結果はすでに元の FPS の Montgomery 表現なので、出力前の全係数の倍率変換を省きます。

append_scaled(a,out,scale) は、通常の内部ワード列 a に scale を掛け、結果を Mint の内部ワードとして out の末尾へ追加する補助関数です。a と out は重ねず、a の値と scale は \([0,p)\) とします。入力は変更しません。

公開 FPS 側では 512 項以上でこの部品を使い、疎な入力に対する既存の漸化式を優先します。2 冪の直後の少数項は、倍の長さの変換を増やさず内積で補います。

作業領域は RAII の Arena が所有し、mmap が使えない場合は通常の整列配列へ戻ります。確保できなければ例外を送出します。

計算量: native_inverse は \(O(n\log n)\) 時間、\(O(n)\) 追加空間です。append_scaled は追加する項数に対して \(O(n)\) 時間です。

依存ライブラリ

• FPS のブロック変換
• FPS の結果配列への一括変換
• FPS の作業領域
• modint998244353 の共通剰余演算
• modint998244353

検証コード

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

コード

#ifndef CP_COMPONENTS_FPS_INVERSE_BLOCKED_HPP
#define CP_COMPONENTS_FPS_INVERSE_BLOCKED_HPP

#include "fps_blocked.hpp"

#if defined(__linux__) && defined(__x86_64__) && defined(__GNUC__) && !defined(__clang__)
#pragma GCC push_options
_Pragma("GCC target(\"avx2\")")
#pragma GCC optimize("O3")

namespace fps_detail::blocked_poly {
inline void append_scaled(Words a, std::vector<mint> &out, Word scale) {
	const detail::Factor factor(scale);
	if (a.size() > out.max_size() - out.size()) throw std::length_error("fps output");
	prepare_fps_output(out, out.size() + a.size());
	alignas(32) Word values[64];
	for (size_t i = 0; i < a.size(); i += 64) {
		const auto count = std::min<size_t>(64, a.size() - i);
		size_t j = 0;
		for (; j + 8 <= count; j += 8)
			detail::store(values + j, detail::reduce(detail::times(detail::load_unaligned(a.data() + i + j), factor), mint::mod()));
		for (; j < count; ++j) values[j] = blocked_ntt::detail::multiply_mod(a[i + j], scale);
		append_fps_words<true>(out, values, count);
	}
}

inline void native_inverse(const mint *in, size_t input, std::vector<mint> &out, size_t n) {
	const int lg = inverse_log(n);
	const size_t len = size_t(1) << lg, k = len / 2;
	const size_t c = std::min(input, n), scratch_words = std::max(len, inverse_scratch(k));
	Arena arena(Transform::words(lg) + Arena::footprint(len) + Arena::footprint(scratch_words) + Arena::footprint(k));
	const Transform t(arena, lg);
	auto a = arena.take(len), scratch = arena.take(scratch_words), g = arena.take(k);
	copy_coefficients(in, a.first(c));
	const ConstWords f = a.first(c);
	const size_t target = std::min(n, k), base = std::min(target, detail::kInverseBase);
	detail::inverse_direct(f, g.first(base));
	for (size_t i = 0; i < base; ++i) g[i] = blocked_ntt::detail::multiply_mod(g[i], 932051910);
	for (size_t count = base; count < target; count *= 2) {
		const size_t size = 2 * count;
		inverse_step(t, f, g.first(count), g.subspan(count, std::min(count, target - count)), scratch.first(size),
				scratch.subspan(Arena::footprint(size), size), 53925088);
	}
	if (n <= k) return append_fps_words<true>(out, g.data(), n);
	inverse_step(t, f, g, a.subspan(k, n - k), scratch.first(len), a, 53925088);
	prepare_fps_output(out, n);
	append_fps_words<true>(out, g.data(), k);
	append_fps_words<true>(out, a.data() + k, n - k);
}
}

#pragma GCC pop_options
#endif
#endif

Back to top page