結果
問題 | No.502 階乗を計算するだけ |
ユーザー |
|
提出日時 | 2024-01-30 21:07:10 |
言語 | C++17 (gcc 13.3.0 + boost 1.87.0) |
結果 |
TLE
|
実行時間 | - |
コード長 | 12,848 bytes |
コンパイル時間 | 2,850 ms |
コンパイル使用メモリ | 127,480 KB |
最終ジャッジ日時 | 2025-02-19 00:48:20 |
ジャッジサーバーID (参考情報) |
judge4 / judge1 |
(要ログイン)
ファイルパターン | 結果 |
---|---|
other | AC * 42 TLE * 10 |
ソースコード
#define PROBLEM "https://yukicoder.me/problems/no/502"#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];}template <typename ...Ds, std::enable_if_t<std::conjunction_v<std::is_integral<Ds>...>, std::nullptr_t> = nullptr>U polynom(const int n, const Ds& ...ds) {if (n < 0) return 0;ensure(n);int sumd = 0;U res = _fac[n];for (int d : { ds... }) {if (d < 0 or d > n) return 0;sumd += d;res *= _fac_inv[d];}if (sumd > n) return 0;res *= _fac_inv[n - sumd];return res;}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 suisennamespace 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 suisennamespace 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^24static constexpr long long MOD2 = 167772161; // 2^25static constexpr long long MOD3 = 469762049; // 2^26static 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^24static constexpr long long MOD2 = 167772161; // 2^25static constexpr long long MOD3 = 469762049; // 2^26static 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 suisennamespace suisen {// mod must be a prime numbertemplate <typename mint,std::enable_if_t<atcoder::internal::is_static_modint<mint>::value, std::nullptr_t> = nullptr>struct factorial_large {using value_type = mint;static constexpr int LOG_BLOCK_SIZE = 9;static constexpr int BLOCK_SIZE = 1 << LOG_BLOCK_SIZE;static constexpr int BLOCK_NUM = value_type::mod() >> LOG_BLOCK_SIZE;static inline int threshold = 2000000;static_assert(atcoder::internal::is_prime_constexpr(mint::mod()));static value_type fac(int n) {return n <= threshold ? factorial<mint>{}.fac(n) : _large_fac(n);}static value_type fac_inv(int n) {return n <= threshold ? factorial<mint>{}.fac_inv(n) : _large_fac(n).inv();}static value_type binom(int n, int r) {if (r < 0 or r > n) return 0;return fac(n) * fac_inv(r) * fac_inv(n - r);}template <typename ...Ds, std::enable_if_t<std::conjunction_v<std::is_integral<Ds>...>, std::nullptr_t> = nullptr>static value_type polynom(const int n, const Ds& ...ds) {if (n < 0) return 0;long long sumd = 0;value_type res = fac(n);for (int d : { ds... }) {if (d < 0 or d > n) return 0;sumd += d;res *= fac_inv(d);}if (sumd > n) return 0;res *= fac_inv(n - sumd);return res;}static value_type perm(int n, int r) {if (r < 0 or r > n) return 0;return fac(n) * fac_inv(n - r);}private:static inline std::vector<value_type> _block_fac{};static void _build() {if (_block_fac.size()) {return;}std::vector<value_type> f{ 1 };f.reserve(BLOCK_SIZE);for (int i = 0; i < LOG_BLOCK_SIZE; ++i) {std::vector<value_type> g = shift_of_sampling_points<value_type>(f, 1 << i, 3 << i, arbitrary_mod_convolution<value_type>);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) {f[j] = get(2 * j) * get(2 * j + 1) * ((2 * j + 1) << i);}}// f_B(x) = (x+1) * ... * (x+B-1)if (BLOCK_NUM > BLOCK_SIZE) {std::vector<value_type> g = shift_of_sampling_points<value_type>(f, BLOCK_SIZE, BLOCK_NUM - BLOCK_SIZE, arbitrary_mod_convolution<value_type>);std::move(g.begin(), g.end(), std::back_inserter(f));} else {f.resize(BLOCK_NUM);}for (int i = 0; i < BLOCK_NUM; ++i) {f[i] *= value_type(i + 1) * BLOCK_SIZE;}// f[i] = (i*B + 1) * ... * (i*B + B)f.insert(f.begin(), 1);for (int i = 1; i <= BLOCK_NUM; ++i) {f[i] *= f[i - 1];}_block_fac = std::move(f);}static value_type _large_fac(int n) {_build();value_type res;int q = n / BLOCK_SIZE, r = n % BLOCK_SIZE;if (2 * r <= BLOCK_SIZE) {res = _block_fac[q];for (int i = 0; i < r; ++i) {res *= value_type::raw(n - i);}} else if (q != factorial_large::BLOCK_NUM) {res = _block_fac[q + 1];value_type den = 1;for (int i = 1; i <= BLOCK_SIZE - r; ++i) {den *= value_type::raw(n + i);}res /= den;} else {// Wilson's theoremres = value_type::mod() - 1;value_type den = 1;for (int i = value_type::mod() - 1; i > n; --i) {den *= value_type::raw(i);}res /= den;}return res;}};} // namespace suisenint main() {int n;std::cin >> n;std::cout << suisen::factorial_large<mint>{}.fac(n).val() << '\n';}