Factorial mod 998244353

math/factorial_mod.hpp · GitHub で見る

共通ページ · 軽量実装なし

大きな非負整数 \(N\) に対して、\(N! \bmod 998244353\) を求めます。法を固定して用意した小さな表から最寄りの階乗を取り出し、短い区間の積で補正します。実行時の前計算は不要です。

#include "lib/math/factorial_mod.hpp"
#include <cassert>
#include <cstdint>

int main() {
	assert(factorial_mod(0).val() == 1);
	assert(factorial_mod(10).val() == 3628800);
	assert(factorial_mod(998244352).val() == 998244352);
	assert(factorial_mod(998244353).val() == 0);
	assert(factorial_mod(std::uint64_t(1) << 60).val() == 0);
}

factorial_mod(std::uint64_t n) は、剰余を既存の mint 型で返します。整数値は .val() で取得できます。引数に既定値はなく、負の値は渡さないでください。std::uint64_t の全域を受け付け、n >= 998244353 なら即座に 0 を返します。0! は 1 です。別の法や動的法には対応しません。

保持する表は components/factorial_mod_table.hpp の 15233 個の uint32_t で、table[k] は (65536*k)! mod 998244353 です。最後の項は (998244352)! で、値は 998244352 になります。表は定数として約 60 KiB を使用し、起動時の生成やヒープ領域の確保はありません。

表の間隔を \(B=65536\)、法を \(p=998244353\) とすると、n < p の計算時間は \(O(B+\log p)\)、最大 32768 個の整数の積と、必要なら 1 回の逆元計算です。追加空間は \(O(1)\)、固定表の保持は \(O(p/B)\) です。n >= p と表に一致する値は定数時間です。多数の連続した小さい添字に対する階乗・逆階乗の全表が必要な場合は、Comb の階乗表も利用できます。

依存ライブラリ

• Factorial Mod Table
• modint998244353

検証コード

• test/factorial_mod.test.cpp

コード

#ifndef FACTORIAL_MOD_HPP
#define FACTORIAL_MOD_HPP

#include "mint.hpp"
#include "../components/factorial_mod_table.hpp"
#include <algorithm>
#include <cstdint>

inline mint factorial_mod(std::uint64_t n) {
	if (n >= mint::mod()) return 0;
	constexpr std::uint32_t step = 65536;
	std::uint32_t block = std::uint32_t((n + step / 2) / step), anchor = block * step;
	std::uint32_t first = std::uint32_t(std::min<std::uint64_t>(n, anchor)) + 1;
	std::uint32_t last = std::uint32_t(std::max<std::uint64_t>(n, anchor));
	std::uint64_t a = 1, b = 1, c = 1, d = 1;
	for (; last >= first && last - first >= 3; first += 4) {
		a = a * first % 998244353;
		b = b * (first + 1) % 998244353;
		c = c * (first + 2) % 998244353;
		d = d * (first + 3) % 998244353;
	}
	for (; first <= last; ++first) a = a * first % 998244353;
	mint product = mint::raw(std::uint32_t((a * b % 998244353) * (c * d % 998244353) % 998244353));
	mint result = mint::raw(factorial_mod_detail::table[block]);
	return anchor > n ? result / product : result * product;
}

#endif

Back to top page