多項式 Half-GCD（軽量版）

components/polynomial_half_gcd_light.hpp · GitHub で見る

64 KiBサイト · 軽量実装

高速実装（256 KiB向け）はこちら

fps_light_detail::PolynomialHalfGCD<T> は、軽量多項式 GCD が使う内部部品です。通常は公開する多項式 GCD のヘッダーを使います。F は LightFormalPowerSeries<T>、係数の法は 998244353 です。

型・関数	内容・返り値
Matrix	成分 a を行優先で持つ 2×2 多項式行列。初期値は単位行列
Matrix::identity	単位行列として処理を省略できるか
high(a,k)	定数項側から k 係数を除いた多項式。k >= 0
multiply(a,b)	行列積 a*b
step(matrix,q)	互除法の行列 [0,1;1,-q] を左から掛ける
apply(matrix,a,b)	行列を多項式ペアに作用させて上書きし、末尾の 0 を除く
half(a,b)	後ろの多項式の係数数を元の a.size()/2 以下にする変換行列
reduce(a,b,matrix = nullptr)	b が零多項式になるまで縮小。matrix を渡すと同じ変換を左から累積する

half は末尾の 0 を除いた多項式について a.size() > b.size() を前提とします。b は空でも構いません。作用後は先頭の多項式の係数数が元の半分より大きく、後ろは半分以下になります。行列式は ±1 です。reduce の入力も末尾の 0 を除きます。両方が空の入力も扱えます。

Matrix::a を直接書き換える場合は identity を正しく設定してください。NTT キャッシュや専用 SIMD 処理を持つ高速実装はこちら。

最大係数数を \(n\)、畳み込みの計算量を \(M(n)\) として、half・reduce は \(O(M(n)\log n)\)、行列積・行列の作用・step は \(O(M(n))\)、high は \(O(n)\) 時間です。追加空間は \(O(n)\) です。

依存ライブラリ

• FPS の互換ヘッダー（軽量版）
• 軽量 FPS の基本演算
• NTT（998244353）（軽量版）
• 軽量 FPS の疎演算
• 疎 FPS の項収集・演算種別
• 多項式の積・打ち切り積の共通処理
• 多項式の共通積木
• 形式的冪級数（軽量版）
• modint998244353

検証コード

• test/fps_light_polynomial_extended_gcd.test.cpp
• test/fps_light_polynomial_gcd.test.cpp
• test/fps_light_polynomial_inv_mod.test.cpp

コード

#ifndef CP_COMPONENTS_POLYNOMIAL_HALF_GCD_LIGHT_HPP
#define CP_COMPONENTS_POLYNOMIAL_HALF_GCD_LIGHT_HPP

#include "../fps/fps_light.hpp"
#include <array>

namespace fps_light_detail {
template<class T>
struct PolynomialHalfGCD {
	using F = LightFormalPowerSeries<T>;
	struct Matrix {
		std::array<F, 4> a{F{1}, F{}, F{}, F{1}};
		bool identity = true;
	};
	static F high(const F &a, int k) {
		assert(k >= 0);
		return k >= int(a.size()) ? F{} : F(a.begin() + k, a.end());
	}
	static Matrix multiply(const Matrix &a, const Matrix &b) {
		if (a.identity) return b;
		if (b.identity) return a;
		Matrix result;
		result.identity = false;
		for (int i = 0; i < 2; ++i) for (int j = 0; j < 2; ++j) {
			result.a[2 * i + j] = a.a[2 * i] * b.a[j] + a.a[2 * i + 1] * b.a[2 + j];
			result.a[2 * i + j].shrink();
		}
		return result;
	}
	static void step(Matrix &matrix, const F &q) {
		F c = matrix.a[0] - q * matrix.a[2], d = matrix.a[1] - q * matrix.a[3];
		c.shrink(), d.shrink();
		matrix.a[0] = std::move(matrix.a[2]), matrix.a[1] = std::move(matrix.a[3]);
		matrix.a[2] = std::move(c), matrix.a[3] = std::move(d);
		matrix.identity = false;
	}
	static void apply(const Matrix &matrix, F &a, F &b) {
		if (matrix.identity) return;
		F c = matrix.a[0] * a + matrix.a[1] * b;
		b = matrix.a[2] * a + matrix.a[3] * b, a = std::move(c);
		a.shrink(), b.shrink();
	}
	static Matrix half(F a, F b) {
		assert(a.size() > b.size());
		const int m = int(a.size()) / 2;
		Matrix result;
		if (int(b.size()) <= m) return result;
		if (a.size() <= 64) {
			while (int(b.size()) > m) {
				auto [q, r] = a.divmod(b);
				step(result, q), a = std::move(b), b = std::move(r);
			}
			return result;
		}
		result = half(high(a, m), high(b, m));
		apply(result, a, b);
		if (int(b.size()) <= m) return result;
		auto [q, r] = a.divmod(b);
		step(result, q), a = std::move(b), b = std::move(r);
		if (int(b.size()) <= m) return result;
		const int k = 2 * m + 1 - int(a.size());
		return multiply(half(high(a, k), high(b, k)), result);
	}
	static void reduce(F &a, F &b, Matrix *matrix = nullptr) {
		while (!b.empty()) {
			if (a.size() > b.size()) {
				auto part = half(a, b);
				apply(part, a, b);
				if (matrix) *matrix = multiply(part, *matrix);
				if (b.empty()) break;
			}
			auto [q, r] = a.divmod(b);
			if (matrix) step(*matrix, q);
			a = std::move(b), b = std::move(r);
		}
	}
};
}

#endif

Back to top page