線形漸化式の第 k 項（高速版）

fps/linear_recurrence.hpp · GitHub で見る

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

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

#include "lib/fps/linear_recurrence.hpp"

std::vector<mint> initial{0, 1}, coefficients{1, 1};
auto value = linear_recurrence(initial, coefficients, 10); // 55

初期値と係数は vector<T> で渡します。法は 998244353 です。

係数が未知で先頭の数項が分かっている場合は、線形漸化式の復元で係数を求められます。

連続する複数の項は 有理 FPS の区間係数取得 でまとめて計算できます。

initial は {a[0],a[1],...,a[d-1]} の順、coefficients[0] は一つ前の項に掛ける係数です。例は a[0]=0, a[1]=1, a[n]=a[n-1]+a[n-2] の Fibonacci 数列です。k は 0 始まりなので、10 を渡すと a[10]=55 を返します。戻り値も mint です。

操作	内容・返り値	計算量
linear_recurrence(initial,c,k)	\(a_n=c_0a_{n-1}+\cdots+c_{d-1}a_{n-d}\) の \(a_k\)。初期値・係数は各 \(d\) 個、\(k\ge0\)。\(d=0\) は 0	\(O(d\log d\log k)\)

係数末尾の零は内部で除き、その次数に対応する末尾の初期値を使います。指定済みの初期値と、漸化式を適用し始める位置は維持します。

漸化式の次数を \(d\) として、時間 \(O(d\log d\log k)\)、追加空間 \(O(d)\) です。

依存ライブラリ

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

検証コード

• test/fps_linear_recurrence.test.cpp

コード

#ifndef CP_FPS_LINEAR_RECURRENCE_HPP
#define CP_FPS_LINEAR_RECURRENCE_HPP

#include "bostan_mori.hpp"

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

template<class T>
T linear_recurrence(const std::vector<T> &initial, const std::vector<T> &coefficients, long long k) {
	using F = FormalPowerSeries<T>;
	using O = fps_detail::Operations<T>;
	assert(k >= 0 && initial.size() == coefficients.size());
	const int d = static_cast<int>(initial.size());
	if (k < d) return initial[k];
	int order = d;
	auto zero = [](T x) {
		if constexpr (std::is_same_v<T, mint>) return x == T(0);
		else return x.val() == 0;
	};
	while (order && zero(coefficients[order - 1])) --order;
	if (order == 0) return 0;
	if constexpr (!std::is_same_v<T, mint>) {
		typename O::W a(order), c(order);
		for (int i = 0; i < order; ++i) {
			a[i] = mint::raw(initial[d - order + i].val());
			c[i] = mint::raw(coefficients[i].val());
		}
		return T::raw(linear_recurrence<mint>(a, c, k - d + order).val());
	}
	if (order == 1) return initial[d - 1] * coefficients[0].pow(k - d + 1);
	F q(order + 1);
	q[0] = 1;
	for (int i = 0; i < order; ++i) q[i + 1] = -coefficients[i];
	F p = F::multiply(F(initial.end() - order, initial.end()), q).pre(order);
	return bostan_mori(std::move(p), std::move(q), k - d + order);
}

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

#endif

Back to top page