結果

問題 No.502 階乗を計算するだけ
ユーザー suisen
提出日時 2023-11-12 23:23:49
言語 C++17
(gcc 13.3.0 + boost 1.87.0)
結果
AC  
実行時間 71 ms / 1,000 ms
コード長 11,113 bytes
コンパイル時間 2,838 ms
コンパイル使用メモリ 128,128 KB
最終ジャッジ日時 2025-02-17 21:46:45
ジャッジサーバーID
(参考情報)
judge4 / judge1
このコードへのチャレンジ
(要ログイン)
ファイルパターン 結果
other AC * 52
権限があれば一括ダウンロードができます

ソースコード

diff #
プレゼンテーションモードにする

#define PROBLEM "https://judge.yosupo.jp/problem/factorial"
#include <iostream>
#include <atcoder/modint>
using mint = atcoder::modint1000000007;
#include <utility>
#include <vector>
#include <atcoder/convolution>
#include <cassert>
namespace suisen {
template <typename T, typename U = T>
struct factorial {
factorial() = default;
factorial(int n) { ensure(n); }
static void ensure(const int n) {
int sz = _fac.size();
if (n + 1 <= sz) return;
int new_size = std::max(n + 1, sz * 2);
_fac.resize(new_size), _fac_inv.resize(new_size);
for (int i = sz; i < new_size; ++i) _fac[i] = _fac[i - 1] * i;
_fac_inv[new_size - 1] = U(1) / _fac[new_size - 1];
for (int i = new_size - 1; i > sz; --i) _fac_inv[i - 1] = _fac_inv[i] * i;
}
T fac(const int i) {
ensure(i);
return _fac[i];
}
T operator()(int i) {
return fac(i);
}
U fac_inv(const int i) {
ensure(i);
return _fac_inv[i];
}
U binom(const int n, const int r) {
if (n < 0 or r < 0 or n < r) return 0;
ensure(n);
return _fac[n] * _fac_inv[r] * _fac_inv[n - r];
}
U perm(const int n, const int r) {
if (n < 0 or r < 0 or n < r) return 0;
ensure(n);
return _fac[n] * _fac_inv[n - r];
}
private:
static std::vector<T> _fac;
static std::vector<U> _fac_inv;
};
template <typename T, typename U>
std::vector<T> factorial<T, U>::_fac{ 1 };
template <typename T, typename U>
std::vector<U> factorial<T, U>::_fac_inv{ 1 };
} // namespace suisen
namespace suisen {
template <typename mint, typename Convolve,
std::enable_if_t<std::is_invocable_r_v<std::vector<mint>, Convolve, std::vector<mint>, std::vector<mint>>, std::nullptr_t> = nullptr>
std::vector<mint> shift_of_sampling_points(const std::vector<mint>& ys, mint t, int m, const Convolve &convolve) {
const int n = ys.size();
factorial<mint> fac(std::max(n, m));
std::vector<mint> b = [&] {
std::vector<mint> f(n), g(n);
for (int i = 0; i < n; ++i) {
f[i] = ys[i] * fac.fac_inv(i);
g[i] = (i & 1 ? -1 : 1) * fac.fac_inv(i);
}
std::vector<mint> b = convolve(f, g);
b.resize(n);
return b;
}();
std::vector<mint> e = [&] {
std::vector<mint> c(n);
mint prd = 1;
std::reverse(b.begin(), b.end());
for (int i = 0; i < n; ++i) {
b[i] *= fac.fac(n - i - 1);
c[i] = prd * fac.fac_inv(i);
prd *= t - i;
}
std::vector<mint> e = convolve(b, c);
e.resize(n);
return e;
}();
std::reverse(e.begin(), e.end());
for (int i = 0; i < n; ++i) {
e[i] *= fac.fac_inv(i);
}
std::vector<mint> f(m);
for (int i = 0; i < m; ++i) f[i] = fac.fac_inv(i);
std::vector<mint> res = convolve(e, f);
res.resize(m);
for (int i = 0; i < m; ++i) res[i] *= fac.fac(i);
return res;
}
template <typename mint>
std::vector<mint> shift_of_sampling_points(const std::vector<mint>& ys, mint t, int m) {
auto convolve = [&](const std::vector<mint> &f, const std::vector<mint> &g) { return atcoder::convolution(f, g); };
return shift_of_sampling_points(ys, t, m, convolve);
}
} // namespace suisen
#include <atcoder/convolution>
#include <iostream>
namespace suisen::internal {
template <typename T, typename R = T>
std::vector<R> convolution_naive(const std::vector<T>& a, const std::vector<T>& b) {
const int n = a.size(), m = b.size();
std::vector<R> c(n + m - 1);
if (n < m) {
for (int j = 0; j < m; j++) for (int i = 0; i < n; i++) c[i + j] += R(a[i]) * b[j];
} else {
for (int i = 0; i < n; i++) for (int j = 0; j < m; j++) c[i + j] += R(a[i]) * b[j];
}
return c;
}
} // namespace suisen
namespace suisen {
template <typename mint, atcoder::internal::is_modint_t<mint>* = nullptr>
std::vector<mint> arbitrary_mod_convolution(const std::vector<mint>& a, const std::vector<mint>& b) {
int n = int(a.size()), m = int(b.size());
if constexpr (atcoder::internal::is_static_modint<mint>::value) {
int maxz = 1;
while (not ((mint::mod() - 1) & maxz)) maxz <<= 1;
int z = 1;
while (z < n + m - 1) z <<= 1;
if (z <= maxz) return atcoder::convolution<mint>(a, b);
}
if (n == 0 or m == 0) return {};
if (std::min(n, m) <= 120) return internal::convolution_naive(a, b);
static constexpr long long MOD1 = 754974721; // 2^24
static constexpr long long MOD2 = 167772161; // 2^25
static constexpr long long MOD3 = 469762049; // 2^26
static constexpr long long M1M2 = MOD1 * MOD2;
static constexpr long long INV_M1_MOD2 = atcoder::internal::inv_gcd(MOD1, MOD2).second;
static constexpr long long INV_M1M2_MOD3 = atcoder::internal::inv_gcd(M1M2, MOD3).second;
std::vector<int> a2(n), b2(m);
for (int i = 0; i < n; ++i) a2[i] = a[i].val();
for (int i = 0; i < m; ++i) b2[i] = b[i].val();
auto c1 = atcoder::convolution<MOD1>(a2, b2);
auto c2 = atcoder::convolution<MOD2>(a2, b2);
auto c3 = atcoder::convolution<MOD3>(a2, b2);
const long long m1m2 = mint(M1M2).val();
std::vector<mint> c(n + m - 1);
for (int i = 0; i < n + m - 1; ++i) {
// Garner's Algorithm
// X = x1 + x2 * m1 + x3 * m1 * m2
// x1 = c1[i], x2 = (c2[i] - x1) / m1 (mod m2), x3 = (c3[i] - x1 - x2 * m1) / m2 (mod m3)
long long x1 = c1[i];
long long x2 = (atcoder::static_modint<MOD2>(c2[i] - x1) * INV_M1_MOD2).val();
long long x3 = (atcoder::static_modint<MOD3>(c3[i] - x1 - x2 * MOD1) * INV_M1M2_MOD3).val();
c[i] = x1 + x2 * MOD1 + x3 * m1m2;
}
return c;
}
std::vector<__uint128_t> convolution_int(const std::vector<int> &a, const std::vector<int> &b) {
int n = int(a.size()), m = int(b.size());
auto check_nonnegative = [](int e) { return e >= 0; };
assert(std::all_of(a.begin(), a.end(), check_nonnegative));
assert(std::all_of(b.begin(), b.end(), check_nonnegative));
if (n == 0 or m == 0) return {};
if (std::min(n, m) <= 120) return internal::convolution_naive<int, __uint128_t>(a, b);
static constexpr long long MOD1 = 754974721; // 2^24
static constexpr long long MOD2 = 167772161; // 2^25
static constexpr long long MOD3 = 469762049; // 2^26
static constexpr long long M1M2 = MOD1 * MOD2;
static constexpr long long INV_M1_MOD2 = atcoder::internal::inv_gcd(MOD1, MOD2).second;
static constexpr long long INV_M1M2_MOD3 = atcoder::internal::inv_gcd(M1M2, MOD3).second;
auto c1 = atcoder::convolution<MOD1>(a, b);
auto c2 = atcoder::convolution<MOD2>(a, b);
auto c3 = atcoder::convolution<MOD3>(a, b);
std::vector<__uint128_t> c(n + m - 1);
for (int i = 0; i < n + m - 1; ++i) {
// Garner's Algorithm
// X = x1 + x2 * m1 + x3 * m1 * m2
// x1 = c1[i], x2 = (c2[i] - x1) / m1 (mod m2), x3 = (c3[i] - x1 - x2 * m1) / m2 (mod m3)
int x1 = c1[i];
int x2 = (atcoder::static_modint<MOD2>(c2[i] - x1) * INV_M1_MOD2).val();
int x3 = (atcoder::static_modint<MOD3>(c3[i] - x1 - x2 * MOD1) * INV_M1M2_MOD3).val();
c[i] = x1 + x2 * MOD1 + __uint128_t(x3) * M1M2;
}
return c;
}
} // namespace suisen
namespace suisen {
template <typename mint, atcoder::internal::is_static_modint_t<mint>* = nullptr>
struct FactorialLargePrimeMod {
private:
static constexpr int _p = mint::mod();
static constexpr int _log_b = 15;
static constexpr int _b = 1 << _log_b;
static constexpr int _q = _p >> _log_b;
public:
static constexpr int block_size = _b;
static constexpr int log_block_size = _log_b;
static constexpr int block_num = _q;
FactorialLargePrimeMod() = delete;
static mint fac(long long n) {
if (_p <= n) return 0;
build();
const int q = n >> _log_b, r = n & (_b - 1);
// n! = (qb)! * (n-r+1)(n-r+2)...(n)
mint ans = _block_fact[q];
for (int j = 0; j < r; ++j) {
ans *= mint::raw(n - j);
}
return ans;
}
private:
static inline std::vector<mint> _block_fact{};
static inline bool _built = false;
static void build() {
if (std::exchange(_built, true)) return;
const auto convolve = [&](const std::vector<mint> &f, const std::vector<mint> &g) {
return arbitrary_mod_convolution(f, g);
};
// f_d(x) := (dx+1)*...*(dx+d-1)
// Suppose that we have f_d(0),...,f_d(d-1). (Note that (deg f_d)+1=d)
// f_{2d}(x) = ((2dx+1)*...*(2dx+d-1)) * (2dx+d) * (((2dx+d)+1)* ...*((2dx+d)+d-1))
// = f_d(2x) * f_d(2x+1) * (2dx+d)
// We can calculate f_{2d}(0), ..., f_{2d}(2d-1) from f_d(0), f_d(1), ..., f_d(4d-2), f_d(4d-1)
std::vector<mint> f{ 1 };
f.reserve(_b);
for (int i = 0; i < _log_b; ++i) {
std::vector<mint> g = shift_of_sampling_points<mint>(f, 1 << i, 3 << i, convolve);
const auto get = [&](int j) { return j < (1 << i) ? f[j] : g[j - (1 << i)]; };
f.resize(2 << i);
for (int j = 0; j < 2 << i; ++j) {
// (2*j+1)*2^i <= 2^(2*_log_b) + 2^(_log_b-1) < 2^31 holds if _log_b <= 15
f[j] = get(2 * j) * get(2 * j + 1) * ((2 * j + 1) << i);
}
}
// f_B(x) = (x+1) * ... * (x+B-1)
if (_q > _b) {
std::vector<mint> g = shift_of_sampling_points<mint>(f, _b, _q - _b, convolve);
std::move(g.begin(), g.end(), std::back_inserter(f));
} else {
f.resize(_q);
}
for (int i = 0; i < _q; ++i) {
f[i] *= mint(i + 1) * _b;
}
// f[i] = (i*B + 1) * ... * (i*B + B)
_block_fact = std::move(f);
_block_fact.insert(_block_fact.begin(), 1);
for (int i = 1; i <= _q; ++i) {
_block_fact[i] *= _block_fact[i - 1];
}
}
};
} // namespace suisen
int main() {
using Factorial = suisen::FactorialLargePrimeMod<mint>;
int n;
std::cin >> n;
std::cout << Factorial::fac(n).val() << '\n';
}
הההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההה
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
0