多項式 GCD（高速版）

fps/polynomial_gcd.hpp · GitHub で見る

全体サイト（256 KiB向け） · 高速実装

軽量実装（64 KiB）はこちら

#include "lib/fps/polynomial_gcd.hpp"

fps a{mint(-1), 0, 1};
fps b{mint(-1), 1};
auto g = polynomial_gcd(a, b);              // x-1
auto [d, u, v] = polynomial_extended_gcd(a, b);
auto inverse = polynomial_inv_mod(fps{1, 1}, fps{1, 0, 1});
if (inverse) {
	fps h = *inverse;
	auto remainder = (fps{1, 1} * h) % fps{1, 0, 1}; // {1}
}

関数・型	内容・返り値	計算量
polynomial_gcd(a,b)	最高次係数を 1 にした最大公約多項式	\(O(n\log^2 n)\)
polynomial_extended_gcd(a,b)	PolynomialExtendedGCD<T>{gcd,x,y}。a*x+b*y=gcd を満たす	\(O(n\log^2 n)\)
PolynomialExtendedGCD<T>	gcd・x・y はいずれも FormalPowerSeries<T>。構造化束縛が可能	—（型）
polynomial_inv_mod(a,modulus)	a*h ≡ 1 mod modulus を満たす次数最小の代表元 h を optional<FormalPowerSeries<T>> で返す。存在しなければ nullopt	\(O(n\log^2 n)\)

引数は FPS の FormalPowerSeries<T> で、係数は定数項から並べます。法は 998244353 です。既定の fps・mint に加え、同じ法で val()・raw() など FPS の要件を満たす別の modint 型も使えます。整数係数・浮動小数点係数の GCD には対応しません。

入力の末尾の 0 は除き、返り値にも不要な末尾の 0 は含みません。零多項式は空配列です。入力は変更しません。

gcd(0,0) は零多項式とし、拡張 GCD も三つの空配列を返します。片方だけが 0 の場合は、もう一方を最高次係数 1 に正規化します。非零の定数同士の GCD は fps{1} です。

polynomial_inv_mod の法多項式は非零にします。法が非零の定数の場合、代表元は零多項式で、空配列を持つ optional を返します。それ以外では GCD が 1 のときに限り逆元が存在し、返り値の次数は法多項式の次数未満です。optional が値を持つかどうかで存在を判定してください。FPS の f.inv() が求める \(x^n\) を法とする逆数とは別の操作です。

大きい入力では Half-GCD により互除法の操作をまとめ、小さい入力では途中の除算を省いた互除法を使います。変換行列は必要な上段だけを求め、NTT を再利用します。行列の作用には次数の低下を利用した巡回畳み込みを使います。必要な NTT 長は \(2^{23}\) 以下に収めます。最大係数数が \(2^{22}\) 以下ならこの条件を満たします。

入力の最大係数数を \(n\) として、各関数は \(O(n\log^2 n)\) 時間、追加空間は \(O(n)\) です。

依存ライブラリ

• FPS のブロック変換
• FPS の融合畳み込み
• FPS の定数係数畳み込み
• FPS の非対称畳み込み
• FPS と多項式除算の内部変換
• FPS の指数・冪の内部変換
• FPS の一点評価
• FPS 逆数のブロック変換
• FPS の内部演算
• FPS の結果配列への一括変換
• 疎 FPS の共通演算
• 疎 FPS の次数圧縮・前処理
• 疎 FPS の項収集・演算種別
• FPS 平方根のブロック変換
• FPS の作業領域
• modint998244353 の共通剰余演算
• modint998244353 の共有畳み込み
• NTT（998244353）（高速版）
• 多項式 Half-GCD の専用変換
• 多項式 Half-GCD（高速版）
• 形式的冪級数（高速版）
• modint998244353

検証コード

• test/polynomial_inv_mod.test.cpp

コード

#ifndef CP_FPS_POLYNOMIAL_GCD_HPP
#define CP_FPS_POLYNOMIAL_GCD_HPP

#include "fps.hpp"
#include "../components/polynomial_half_gcd.hpp"
#include "../components/polynomial_gcd_blocked.hpp"
#include <optional>

#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC push_options
#pragma GCC optimize("O3")
#endif

namespace fps_detail {
inline PolynomialHalfGCD::Matrix polynomial_gcd_jump(const PolynomialHalfGCD::F& a, const PolynomialHalfGCD::F& b) {
#if defined(__linux__) && defined(__x86_64__) && defined(__GNUC__) && !defined(__clang__)
	if (a.size() >= 129 && a.size() <= (1u << 23) && NTT::vectorized()) {
		PolynomialHalfGCD::Matrix out;
		blocked_poly::native_gcd_matrix(a.data(), a.size(), b.data(), b.size(), 3,
			[&out](const blocked_poly::gcd_detail::Matrix& matrix) {
				out.shift = int(matrix.progress), out.identity = matrix.progress == 0;
				for (int j = 0; j < 2; ++j) {
					out.a[j].clear();
					out.a[j].reserve(matrix.size[j]);
					blocked_poly::division_detail::append_words(out.a[j], matrix.entry[j], matrix.size[j], blocked_ntt::detail::kR);
				}
			});
		return out;
	}
#endif
	return PolynomialHalfGCD::jump(a, b, int(a.size()) - 1, true);
}
}

template<class T>
struct PolynomialExtendedGCD {
	FormalPowerSeries<T> gcd, x, y;
};

template<class T>
FormalPowerSeries<T> polynomial_gcd(FormalPowerSeries<T> a, FormalPowerSeries<T> b) {
	using H = fps_detail::PolynomialHalfGCD;
	using O = fps_detail::Operations<T>;
	if constexpr (!std::is_same_v<T, mint>)
		return O::export_work(polynomial_gcd(O::import_work(a), O::import_work(b)));
	else {
		O::trim(a), O::trim(b);
		if (a.size() < b.size()) a.swap(b);
		if (!b.empty() && a.size() == b.size()) {
			auto r = a.divmod(b).second;
			a = std::move(b), b = std::move(r);
		}
		if (!b.empty()) {
			auto matrix = fps_detail::polynomial_gcd_jump(a, b);
			if (matrix.shift == int(a.size()) - 1) return {T(1)};
			a = H::apply_first(matrix, a, b);
		}
		if (!a.empty()) a /= a.back();
		return a;
	}
}

template<class T>
PolynomialExtendedGCD<T> polynomial_extended_gcd(FormalPowerSeries<T> a, FormalPowerSeries<T> b) {
	using H = fps_detail::PolynomialHalfGCD;
	using O = fps_detail::Operations<T>;
	if constexpr (!std::is_same_v<T, mint>) {
		auto result = polynomial_extended_gcd(O::import_work(a), O::import_work(b));
		return {O::export_work(result.gcd), O::export_work(result.x), O::export_work(result.y)};
	} else {
		O::trim(a), O::trim(b);
		if (a.empty() && b.empty()) return {};
		typename H::Matrix matrix;
		while (!b.empty() && a.size() <= b.size()) {
			auto [q, r] = a.divmod(b);
			H::step(matrix, q), a = std::move(b), b = std::move(r);
		}
		if (!b.empty()) {
			auto part = fps_detail::polynomial_gcd_jump(a, b);
			a = H::apply_first(part, a, b);
			matrix = H::multiply(part, matrix, true);
		}
		T inverse = a.back().inv();
		a *= inverse, matrix.a[0] *= inverse, matrix.a[1] *= inverse;
		return {std::move(a), std::move(matrix.a[0]), std::move(matrix.a[1])};
	}
}

template<class T>
std::optional<FormalPowerSeries<T>> polynomial_inv_mod(FormalPowerSeries<T> a, FormalPowerSeries<T> modulus) {
	fps_detail::Operations<T>::trim(modulus);
	assert(!modulus.empty());
	if (modulus.size() == 1) return FormalPowerSeries<T>{};
	using H = fps_detail::PolynomialHalfGCD;
	using O = fps_detail::Operations<T>;
	if constexpr (!std::is_same_v<T, mint>) {
		auto result = polynomial_inv_mod(O::import_work(a), O::import_work(modulus));
		if (!result) return std::nullopt;
		return O::export_work(*result);
	} else {
#if defined(__linux__) && defined(__x86_64__) && defined(__GNUC__) && !defined(__clang__)
		if (modulus.size() >= 129 && modulus.size() <= (1u << 23) && a.size() <= (1u << 23) && fps_detail::NTT::vectorized()) {
			FormalPowerSeries<T> out;
			if (!fps_detail::blocked_poly::native_inv_mod(a.data(), a.size(), modulus.data(), modulus.size(), out)) return std::nullopt;
			return out;
		}
#endif
		if (a.size() >= modulus.size()) a = a.divmod(modulus).second;
		O::trim(a);
		if (a.empty()) return std::nullopt;
		auto matrix = H::jump(modulus, a, int(modulus.size()) - 1, true);
		if (matrix.shift != int(modulus.size()) - 1) return std::nullopt;
		T constant = (matrix.a[0].empty() ? T(0) : matrix.a[0][0] * modulus[0]) +
			(matrix.a[1].empty() ? T(0) : matrix.a[1][0] * a[0]);
		matrix.a[1] /= constant;
		return std::move(matrix.a[1]);
	}
}

#if defined(__GNUC__) && !defined(__clang__) && (defined(__x86_64__) || defined(__i386__))
#pragma GCC pop_options
#endif

#endif

Back to top page