#define PROBLEM "https://yukicoder.me/problems/no/502" #include #include using mint = atcoder::modint1000000007; #include #include #include #include namespace suisen { template 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 ...>, 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 _fac; static std::vector _fac_inv; }; template std::vector factorial::_fac{ 1 }; template std::vector factorial::_fac_inv{ 1 }; } // namespace suisen namespace suisen { template , Convolve, std::vector, std::vector>, std::nullptr_t> = nullptr> std::vector shift_of_sampling_points(const std::vector& ys, mint t, int m, const Convolve &convolve) { const int n = ys.size(); factorial fac(std::max(n, m)); std::vector b = [&] { std::vector 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 b = convolve(f, g); b.resize(n); return b; }(); std::vector e = [&] { std::vector 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 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 f(m); for (int i = 0; i < m; ++i) f[i] = fac.fac_inv(i); std::vector res = convolve(e, f); res.resize(m); for (int i = 0; i < m; ++i) res[i] *= fac.fac(i); return res; } template std::vector shift_of_sampling_points(const std::vector& ys, mint t, int m) { auto convolve = [&](const std::vector &f, const std::vector &g) { return atcoder::convolution(f, g); }; return shift_of_sampling_points(ys, t, m, convolve); } } // namespace suisen #include #include namespace suisen::internal { template std::vector convolution_naive(const std::vector& a, const std::vector& b) { const int n = a.size(), m = b.size(); std::vector 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 * = nullptr> std::vector arbitrary_mod_convolution(const std::vector& a, const std::vector& b) { int n = int(a.size()), m = int(b.size()); if constexpr (atcoder::internal::is_static_modint::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(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 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(a2, b2); auto c2 = atcoder::convolution(a2, b2); auto c3 = atcoder::convolution(a2, b2); const long long m1m2 = mint(M1M2).val(); std::vector 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(c2[i] - x1) * INV_M1_MOD2).val(); long long x3 = (atcoder::static_modint(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 &a, const std::vector &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(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(a, b); auto c2 = atcoder::convolution(a, b); auto c3 = atcoder::convolution(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(c2[i] - x1) * INV_M1_MOD2).val(); int x3 = (atcoder::static_modint(c3[i] - x1 - x2 * MOD1) * INV_M1M2_MOD3).val(); c[i] = x1 + x2 * MOD1 + __uint128_t(x3) * M1M2; } return c; } } // namespace suisen namespace suisen { // mod must be a prime number template ::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(long long n) { if (n >= mint::mod()) return 0; return n <= threshold ? factorial{}.fac(n) : _large_fac(n); } static value_type fac_inv(long long n) { assert(n < (long long) mint::mod()); return n <= threshold ? factorial{}.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 ...>, 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 _block_fac{}; static void _build() { if (_block_fac.size()) { return; } std::vector f{ 1 }; f.reserve(BLOCK_SIZE); for (int i = 0; i < LOG_BLOCK_SIZE; ++i) { std::vector g = shift_of_sampling_points(f, 1 << i, 3 << i, arbitrary_mod_convolution); 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 g = shift_of_sampling_points(f, BLOCK_SIZE, BLOCK_NUM - BLOCK_SIZE, arbitrary_mod_convolution); 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 theorem res = 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 suisen int main() { long long n; std::cin >> n; std::cout << suisen::factorial_large{}.fac(n).val() << '\n'; }