結果
| 問題 |
No.502 階乗を計算するだけ
|
| コンテスト | |
| ユーザー |
 hly1204
|
| 提出日時 | 2021-07-26 16:22:40 |
| 言語 | C++17(clang) (17.0.6 + boost 1.87.0) |
| 結果 |
AC
|
| 実行時間 | 48 ms / 1,000 ms |
| コード長 | 25,271 bytes |
| コンパイル時間 | 5,077 ms |
| コンパイル使用メモリ | 144,436 KB |
| 実行使用メモリ | 30,556 KB |
| 最終ジャッジ日時 | 2024-07-22 06:54:20 |
| 合計ジャッジ時間 | 7,474 ms |
|
ジャッジサーバーID （参考情報） |
judge5 / judge4 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| other | AC * 52 |
ソースコード
#line 1 "remote_test\\yuki\\math\\502.0.test.cpp"
#define PROBLEM "https://yukicoder.me/problems/no/502"
#include <iostream>
#line 1 "math\\formal_power_series\\arbitrary_modulo_convolution.hpp"
/**
* @brief arbitrary modulo convolution / 任意模数卷积
*
*/
#include <cstdint>
#include <type_traits>
#include <vector>
#line 1 "modint\\Montgomery_modint.hpp"
/**
* @brief Montgomery modint / Montgomery 取模类
* @docs docs/modint/Montgomery_modint.md
*/
#line 12 "modint\\Montgomery_modint.hpp"
namespace lib {
/**
* @brief Montgomery 取模类
* @see https://nyaannyaan.github.io/library/modint/montgomery-modint.hpp
* @author Nyaan
* @tparam mod 为奇数且大于 1
*/
template <std::uint32_t mod>
class MontgomeryModInt {
public:
using i32 = std::int32_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;
using m32 = MontgomeryModInt;
using value_type = u32;
static constexpr u32 get_mod() { return mod; }
static constexpr u32 get_primitive_root_prime() {
u32 tmp[32] = {};
int cnt = 0;
const u32 phi = mod - 1;
u32 m = phi;
for (u32 i = 2; i * i <= m; ++i) {
if (m % i == 0) {
tmp[cnt++] = i;
do { m /= i; } while (m % i == 0);
}
}
if (m != 1) tmp[cnt++] = m;
for (m32 res = 2;; res += 1) {
bool f = true;
for (int i = 0; i < cnt && f; ++i) f &= res.pow(phi / tmp[i]) != 1;
if (f) return u32(res);
}
}
constexpr MontgomeryModInt() = default;
~MontgomeryModInt() = default;
template <typename T, std::enable_if_t<std::is_integral_v<T>, int> = 0>
constexpr MontgomeryModInt(T v) : v_(reduce(u64(v % i32(mod) + i32(mod)) * r2)) {}
constexpr MontgomeryModInt(const m32 &) = default;
constexpr u32 get() const { return norm(reduce(v_)); }
template <typename T, std::enable_if_t<std::is_integral_v<T>, int> = 0>
explicit constexpr operator T() const {
return T(get());
}
constexpr m32 operator-() const {
m32 res;
res.v_ = (mod2 & -(v_ != 0)) - v_;
return res;
}
constexpr m32 inv() const {
i32 x1 = 1, x3 = 0, a = get(), b = mod;
while (b != 0) {
i32 q = a / b, x1_old = x1, a_old = a;
x1 = x3, x3 = x1_old - x3 * q, a = b, b = a_old - b * q;
}
return m32(x1);
}
constexpr m32 &operator=(const m32 &) = default;
constexpr m32 &operator+=(const m32 &rhs) {
v_ += rhs.v_ - mod2;
v_ += mod2 & -(v_ >> 31);
return *this;
}
constexpr m32 &operator-=(const m32 &rhs) {
v_ -= rhs.v_;
v_ += mod2 & -(v_ >> 31);
return *this;
}
constexpr m32 &operator*=(const m32 &rhs) {
v_ = reduce(u64(v_) * rhs.v_);
return *this;
}
constexpr m32 &operator/=(const m32 &rhs) { return operator*=(rhs.inv()); }
friend constexpr m32 operator+(const m32 &lhs, const m32 &rhs) { return m32(lhs) += rhs; }
friend constexpr m32 operator-(const m32 &lhs, const m32 &rhs) { return m32(lhs) -= rhs; }
friend constexpr m32 operator*(const m32 &lhs, const m32 &rhs) { return m32(lhs) *= rhs; }
friend constexpr m32 operator/(const m32 &lhs, const m32 &rhs) { return m32(lhs) /= rhs; }
friend constexpr bool operator==(const m32 &lhs, const m32 &rhs) {
return norm(lhs.v_) == norm(rhs.v_);
}
friend constexpr bool operator!=(const m32 &lhs, const m32 &rhs) {
return norm(lhs.v_) != norm(rhs.v_);
}
friend std::istream &operator>>(std::istream &is, m32 &rhs) {
i32 x;
is >> x;
rhs = m32(x);
return is;
}
friend std::ostream &operator<<(std::ostream &os, const m32 &rhs) { return os << rhs.get(); }
constexpr m32 pow(u64 y) const {
m32 res(1), x(*this);
for (; y != 0; y >>= 1, x *= x)
if (y & 1) res *= x;
return res;
}
private:
static constexpr u32 get_r() {
u32 two = 2, iv = mod * (two - mod * mod);
iv *= two - mod * iv;
iv *= two - mod * iv;
return iv * (mod * iv - two);
}
static constexpr u32 reduce(u64 x) { return (x + u64(u32(x) * r) * mod) >> 32; }
static constexpr u32 norm(u32 x) { return x - (mod & -((mod - 1 - x) >> 31)); }
u32 v_;
static constexpr u32 r = get_r();
static constexpr u32 r2 = -u64(mod) % mod;
static constexpr u32 mod2 = mod << 1;
static_assert((mod & 1) == 1, "mod % 2 == 0\n");
static_assert(-r * mod == 1, "???\n");
static_assert((mod & (3U << 30)) == 0, "mod >= (1 << 30)\n");
static_assert(mod != 1, "mod == 1\n");
};
// 别名
template <std::uint32_t mod>
using MontModInt = MontgomeryModInt<mod>;
} // namespace lib
#line 1 "math\\formal_power_series\\NTT_crt.hpp"
/**
* @brief NTT prime crt / NTT 素数用中国剩余定理
*
*/
#line 11 "math\\formal_power_series\\NTT_crt.hpp"
namespace lib {
template <std::uint32_t M0, std::uint32_t M1, std::uint32_t M2,
std::enable_if_t<(M0 < M1) && ((M0 | M1 | M2) < (1U << 31)) &&
((M0 & M1 & M2 & 1) == 1) && (M0 != M1) && (M0 != M2) && (M1 != M2),
int> = 0>
class NTTCRT3 {
public:
using u32 = std::uint32_t;
using u64 = std::uint64_t;
constexpr NTTCRT3(u32 mod) : m_(mod), M0M1_mod_m_(static_cast<u64>(M0) * M1 % mod) {}
~NTTCRT3() = default;
constexpr u32 operator()(u32 a, u32 b, u32 c) const {
// x mod M0 = a, x mod M1 = b, x mod M2 = c
// a + k0M0 = b + k1M1 => k0 = (b - a) / M0 (mod M1)
// x = a + k0M0 (mod M0M1) => a + k0M0 + k01M0M1 = c + k2M2
// => k01 = (c - (a + k0M0)) / (M0M1) (mod M2)
// => x mod M0M1M2 = a + k0M0 + k01M0M1
u32 k0 = b - a;
if (static_cast<int>(k0) < 0) k0 += M1;
k0 = static_cast<u64>(k0) * M0_inv_M1_ % M1;
u64 d = a + static_cast<u64>(k0) * M0;
u32 k01 = c - d % M2;
if (static_cast<int>(k01) < 0) k01 += M2;
// NTT 模数都小于 (1U << 31) 所以在这里可以使用加法后再取模
return (d + static_cast<u64>(k01) * M0M1_inv_M2_ % M2 * M0M1_mod_m_) % m_;
}
static constexpr u32 get_inv(u32 x, u32 mod) {
u32 res = 1;
for (u32 e = mod - 2; e != 0; e >>= 1) {
if (e & 1) res = static_cast<u64>(res) * x % mod;
x = static_cast<u64>(x) * x % mod;
}
return res;
}
private:
u32 m_, M0M1_mod_m_;
static constexpr u32 M0_inv_M1_ = get_inv(M0, M1);
static constexpr u32 M0M1_inv_M2_ = get_inv(static_cast<u64>(M0) * M1 % M2, M2);
};
} // namespace lib
#line 1 "math\\formal_power_series\\convolution.hpp"
/**
* @brief convolution / 卷积
*
*/
#include <cassert>
#line 11 "math\\formal_power_series\\convolution.hpp"
#line 1 "math\\formal_power_series\\radix_2_NTT.hpp"
/**
* @brief radix-2 NTT / 基-2 数论变换
* @docs docs/math/formal_power_series/radix_2_NTT.md
*/
#include <algorithm>
#line 13 "math\\formal_power_series\\radix_2_NTT.hpp"
#line 1 "traits\\modint.hpp"
/**
* @brief modint traits / 取模类萃取
*
*/
namespace lib {
template <typename mod_t>
struct modint_traits {
using type = typename mod_t::value_type;
static constexpr type get_mod() { return mod_t::get_mod(); }
static constexpr type get_primitive_root_prime() { return mod_t::get_primitive_root_prime(); }
};
} // namespace lib
#line 15 "math\\formal_power_series\\radix_2_NTT.hpp"
namespace lib {
/**
* @note 必须用 NTT 友好的模数！！！
*/
template <typename mod_t>
class NTT {
public:
NTT() = delete;
static void set_root(int len) {
static int lim = 0;
static constexpr mod_t g(modint_traits<mod_t>::get_primitive_root_prime());
if (lim == 0) {
constexpr int offset = 20;
rt.resize(1 << offset);
irt.resize(1 << offset);
rt[0] = irt[0] = mod_t(1);
mod_t g_t = g.pow(modint_traits<mod_t>::get_mod() >> (offset + 1)), ig_t = g_t.inv();
rt[1 << (offset - 1)] = g_t, irt[1 << (offset - 1)] = ig_t;
for (int i = offset - 2; i >= 0; --i) rt[1 << i] = (g_t *= g_t), irt[1 << i] = (ig_t *= ig_t);
lim = 1;
}
for (; (lim << 1) < len; lim <<= 1) {
mod_t g = rt[lim], ig = irt[lim];
for (int i = lim + 1, e = lim << 1; i < e; ++i)
rt[i] = rt[i - lim] * g, irt[i] = irt[i - lim] * ig;
}
}
static void dft(int n, mod_t *x) {
for (int j = 0, l = n >> 1; j != l; ++j) {
mod_t u = x[j], v = x[j + l];
x[j] = u + v, x[j + l] = u - v;
}
for (int i = n >> 1; i >= 2; i >>= 1) {
for (int j = 0, l = i >> 1; j != l; ++j) {
mod_t u = x[j], v = x[j + l];
x[j] = u + v, x[j + l] = u - v;
}
for (int j = i, l = i >> 1, m = 1; j != n; j += i, ++m) {
mod_t root = rt[m];
for (int k = 0; k != l; ++k) {
mod_t u = x[j + k], v = x[j + k + l] * root;
x[j + k] = u + v, x[j + k + l] = u - v;
}
}
}
}
static void idft(int n, mod_t *x) {
for (int i = 2; i < n; i <<= 1) {
for (int j = 0, l = i >> 1; j != l; ++j) {
mod_t u = x[j], v = x[j + l];
x[j] = u + v, x[j + l] = u - v;
}
for (int j = i, l = i >> 1, m = 1; j != n; j += i, ++m) {
mod_t root = irt[m];
for (int k = 0; k != l; ++k) {
mod_t u = x[j + k], v = x[j + k + l];
x[j + k] = u + v, x[j + k + l] = (u - v) * root;
}
}
}
mod_t iv(mod_t(n).inv());
for (int j = 0, l = n >> 1; j != l; ++j) {
mod_t u = x[j] * iv, v = x[j + l] * iv;
x[j] = u + v, x[j + l] = u - v;
}
}
static void even_dft(int n, mod_t *x) {
static constexpr mod_t IT(mod_t(2).inv());
for (int i = 0, j = 0; i != n; i += 2, ++j) x[j] = IT * (x[i] + x[i + 1]);
}
static void odd_dft(int n, mod_t *x) {
static constexpr mod_t IT(mod_t(2).inv());
for (int i = 0, j = 0; i != n; i += 2, ++j) x[j] = IT * irt[j] * (x[i] - x[i + 1]);
}
static void dft_doubling(int n, mod_t *x) {
static constexpr mod_t g(modint_traits<mod_t>::get_primitive_root_prime());
std::copy_n(x, n, x + n);
idft(n, x + n);
mod_t k(1), t(g.pow((modint_traits<mod_t>::get_mod() - 1) / (n << 1)));
for (int i = 0; i != n; ++i) x[n + i] *= k, k *= t;
dft(n, x + n);
}
private:
static inline std::vector<mod_t> rt, irt;
};
std::uint32_t get_ntt_len(std::uint32_t n) {
--n;
n |= n >> 1;
n |= n >> 2;
n |= n >> 4;
n |= n >> 8;
return (n | n >> 16) + 1;
}
/**
* @brief 接收一个多项式，返回二进制翻转后的 DFT 序列，即 x(1), x(-1) 等，
* 对于下标 i 和 i^1 必然是两个互为相反数的点值
*
* @tparam mod_t
* @param n
* @param x
*/
template <typename mod_t>
void dft(int n, mod_t *x) {
NTT<mod_t>::set_root(n);
NTT<mod_t>::dft(n, x);
}
/**
* @brief 接收二进制翻转后的 DFT 序列，返回多项式序列 mod (x^n - 1)
*
* @tparam mod_t
* @param n
* @param x
*/
template <typename mod_t>
void idft(int n, mod_t *x) {
NTT<mod_t>::set_root(n);
NTT<mod_t>::idft(n, x);
}
template <typename mod_t>
void dft(std::vector<mod_t> &x) {
NTT<mod_t>::set_root(x.size());
NTT<mod_t>::dft(x.size(), x.data());
}
template <typename mod_t>
void idft(std::vector<mod_t> &x) {
NTT<mod_t>::set_root(x.size());
NTT<mod_t>::idft(x.size(), x.data());
}
} // namespace lib
#line 13 "math\\formal_power_series\\convolution.hpp"
namespace lib {
/**
* @brief NTT 模数卷积
* @tparam mod_t NTT 友好的模数类
*/
template <typename mod_t>
std::vector<mod_t> convolve(const std::vector<mod_t> &x, const std::vector<mod_t> &y) {
int n = x.size(), m = y.size();
if (std::min(n, m) <= 32) {
std::vector<mod_t> res(n + m - 1, mod_t(0));
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) res[i + j] += x[i] * y[j];
return res;
}
int len = get_ntt_len(n + m - 1);
std::vector<mod_t> res(len);
std::copy_n(x.begin(), n, res.begin());
std::fill(res.begin() + n, res.end(), mod_t(0));
dft(res);
if (&x == &y) {
for (int i = 0; i < len; ++i) res[i] *= res[i];
} else {
std::vector<mod_t> y_tmp(len);
std::copy_n(y.begin(), m, y_tmp.begin());
std::fill(y_tmp.begin() + m, y_tmp.end(), mod_t(0));
dft(y_tmp);
for (int i = 0; i < len; ++i) res[i] *= y_tmp[i];
}
idft(res);
res.resize(n + m - 1);
return res;
}
/**
* @brief NTT 模数循环卷积
* @param x
* @param y
* @param cyclen 必须为 2 的幂次！
* @return std::vector<mod_t> convolve(x, y) mod (x^cyclen - 1)
*/
template <typename mod_t>
std::vector<mod_t> convolve_cyclic(const std::vector<mod_t> &x, const std::vector<mod_t> &y,
int cyclen) {
assert((cyclen & (cyclen - 1)) == 0);
int n = x.size(), m = y.size(), mask = cyclen - 1;
if (cyclen >= n + m - 1) return convolve(x, y);
if (std::min(n, m) <= 32) {
std::vector<mod_t> res(cyclen, mod_t(0));
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j) res[(i + j) & mask] += x[i] * y[j];
return res;
}
std::vector<mod_t> res(cyclen, mod_t(0));
for (int i = 0; i < n; ++i) res[i & mask] += x[i];
dft(res);
if (&x == &y) {
for (int i = 0; i < cyclen; ++i) res[i] *= res[i];
} else {
std::vector<mod_t> y_tmp(cyclen, mod_t(0));
for (int i = 0; i < m; ++i) y_tmp[i & mask] += y[i];
dft(y_tmp);
for (int i = 0; i < cyclen; ++i) res[i] *= y_tmp[i];
}
idft(res);
return res;
}
} // namespace lib
#line 16 "math\\formal_power_series\\arbitrary_modulo_convolution.hpp"
namespace lib {
/**
* @brief 任意模数卷积
* @note 只适用于模数为 32 位
*/
template <typename Int>
std::vector<Int> convolve_mod(const std::vector<Int> &x, const std::vector<Int> &y,
std::uint32_t mod) {
using u32 = std::uint32_t;
using Type = std::conditional_t<(sizeof(Int) <= 4), std::uint32_t, std::uint64_t>;
static constexpr u32 M0 = 880803841, M1 = 897581057, M2 = 998244353;
NTTCRT3<M0, M1, M2> crt(mod);
using mod_t0 = MontModInt<M0>;
using mod_t1 = MontModInt<M1>;
using mod_t2 = MontModInt<M2>;
int n = x.size(), m = y.size();
std::vector<mod_t0> x0(n), y0(m);
std::vector<mod_t1> x1(n), y1(m);
std::vector<mod_t2> x2(n), y2(m);
for (int i = 0; i < n; ++i) {
u32 v = static_cast<Type>(x[i]);
x0[i] = v, x1[i] = v, x2[i] = v;
}
for (int i = 0; i < m; ++i) {
u32 v = static_cast<Type>(y[i]);
y0[i] = v, y1[i] = v, y2[i] = v;
}
auto res0 = convolve(x0, y0);
auto res1 = convolve(x1, y1);
auto res2 = convolve(x2, y2);
int nm = res0.size();
std::vector<Int> res;
res.reserve(nm);
for (int i = 0; i < nm; ++i) res.emplace_back(crt(u32(res0[i]), u32(res1[i]), u32(res2[i])));
return res;
}
template <typename Int>
std::vector<Int> convolve_cyclic_mod(const std::vector<Int> &x, const std::vector<Int> &y,
std::uint32_t mod, int cyclen) {
using u32 = std::uint32_t;
using Type = std::conditional_t<(sizeof(Int) <= 4), std::uint32_t, std::uint64_t>;
static constexpr u32 M0 = 880803841, M1 = 897581057, M2 = 998244353;
NTTCRT3<M0, M1, M2> crt(mod);
using mod_t0 = MontModInt<M0>;
using mod_t1 = MontModInt<M1>;
using mod_t2 = MontModInt<M2>;
int n = x.size(), m = y.size();
std::vector<mod_t0> x0(n), y0(m);
std::vector<mod_t1> x1(n), y1(m);
std::vector<mod_t2> x2(n), y2(m);
for (int i = 0; i < n; ++i) {
u32 v = static_cast<Type>(x[i]);
x0[i] = v, x1[i] = v, x2[i] = v;
}
for (int i = 0; i < m; ++i) {
u32 v = static_cast<Type>(y[i]);
y0[i] = v, y1[i] = v, y2[i] = v;
}
auto res0 = convolve_cyclic(x0, y0, cyclen);
auto res1 = convolve_cyclic(x1, y1, cyclen);
auto res2 = convolve_cyclic(x2, y2, cyclen);
int nm = res0.size();
std::vector<Int> res;
res.reserve(nm);
for (int i = 0; i < nm; ++i) res.emplace_back(crt(u32(res0[i]), u32(res1[i]), u32(res2[i])));
return res;
}
} // namespace lib
#line 1 "math\\modulo\\factorial_modulo_prime.hpp"
/**
* @brief factorial modulo prime / 阶乘模素数
* @docs docs/math/modulo/factorial_modulo_prime.md
*/
#line 10 "math\\modulo\\factorial_modulo_prime.hpp"
#include <cmath>
#line 12 "math\\modulo\\factorial_modulo_prime.hpp"
#include <functional>
#include <iterator>
#include <numeric>
#line 16 "math\\modulo\\factorial_modulo_prime.hpp"
#line 1 "math\\formal_power_series\\sample_points_shift.hpp"
/**
* @brief sample points shift / 样本点平移
* @docs docs/math/formal_power_series/sample_points_shift.md
*/
#line 15 "math\\formal_power_series\\sample_points_shift.hpp"
#line 1 "math\\formal_power_series\\prime_binomial.hpp"
/**
* @brief prime binomial / 素数用二项式系数
*
*/
#line 11 "math\\formal_power_series\\prime_binomial.hpp"
namespace lib {
template <typename mod_t>
class PrimeBinomial {
public:
PrimeBinomial(int lim) {
static std::uint64_t mod = 0;
if (mod != mod_t::get_mod()) {
mod = mod_t::get_mod();
fac_table.clear();
ifac_table.clear();
}
if (fac_table.empty()) {
fac_table.emplace_back(1);
ifac_table.emplace_back(1);
}
init(lim);
}
~PrimeBinomial() = default;
/**
* @brief 预处理 [0, n) 的阶乘和其逆元
*/
static void init(int n) {
if (int(fac_table.size()) < n) {
int old_size = fac_table.size();
fac_table.resize(n);
ifac_table.resize(n);
for (int i = old_size; i < n; ++i) fac_table[i] = fac_table[i - 1] * mod_t(i);
ifac_table.back() = mod_t(1) / fac_table.back();
for (int i = n - 2; i >= old_size; --i) ifac_table[i] = ifac_table[i + 1] * mod_t(i + 1);
}
}
mod_t fac_unsafe(int n) const { return fac_table[n]; }
mod_t ifac_unsafe(int n) const { return ifac_table[n]; }
mod_t inv_unsafe(int n) const { return ifac_table[n] * fac_table[n - 1]; }
mod_t choose_unsafe(int n, int k) const {
// 返回 binom{n}{k} 注意上指标可以为负数但这里并未实现！
return n >= k ? fac_table[n] * ifac_table[k] * ifac_table[n - k] : mod_t(0);
}
private:
static inline std::vector<mod_t> fac_table, ifac_table;
};
} // namespace lib
#line 17 "math\\formal_power_series\\sample_points_shift.hpp"
namespace lib::internal {
/**
* @brief 样本点平移（通过拉格朗日插值公式）
* @note 不安全的算法
* @tparam mod_t 素数模数且点数不能超过模数！
* @tparam ConvolveCyclicFuncType
* @param n 返回值的点数
* @param pts f(0), f(1), …, f(k-1) 确定一个唯一的度数小于 k 的多项式
* @param m 平移距离 f(x) => f(x+m)
* @param f 循环卷积函数
* @return std::vector<mod_t> f(m), f(m+1), …, f(m+n-1)
*/
template <typename mod_t, typename ConvolveCyclicFuncType>
std::vector<mod_t> shift_sample_points_unsafe(int n, const std::vector<mod_t> &pts, mod_t m,
ConvolveCyclicFuncType f) {
if (n == 0) return {};
int s = pts.size(), deg_A = s - 1;
PrimeBinomial<mod_t> bi(s);
std::vector<mod_t> A(s), B(deg_A + n), p_sum(deg_A + n);
for (int i = 0; i < s; ++i) {
A[i] = pts[i] * bi.ifac_unsafe(i) * bi.ifac_unsafe(deg_A - i);
if ((deg_A - i) & 1) A[i] = -A[i];
}
for (int i = 0; i < deg_A + n; ++i) B[i] = m + mod_t(i - deg_A);
std::partial_sum(B.begin(), B.end(), p_sum.begin(), std::multiplies<>());
assert(p_sum.back() != mod_t(0));
mod_t p_inv = mod_t(1) / p_sum.back();
for (int i = deg_A + n - 1; i > 0; --i) {
mod_t t = p_inv * B[i];
B[i] = p_inv * p_sum[i - 1];
p_inv = t;
}
B[0] = p_inv;
A = f(A, B, get_ntt_len(s + s - 1 + n - (s < 2 ? 0 : deg_A - 1) - 1));
mod_t coeff(m);
for (int i = 1; i < s; ++i) coeff *= m - mod_t(i);
for (int i = 0; i < n; ++i) A[i] = A[deg_A + i] * coeff, coeff *= (m + mod_t(i + 1)) * B[i];
A.resize(n);
return A;
}
template <typename mod_t, typename ConvolveCyclicFuncType>
std::vector<mod_t> shift_sample_points_unsafe(const std::vector<mod_t> &pts, mod_t m,
ConvolveCyclicFuncType f) {
return shift_sample_points_unsafe(pts.size(), pts, m, f);
}
} // namespace lib::internal
namespace lib {
template <typename mod_t, typename ConvolveCyclicFuncType>
std::vector<mod_t> shift_sample_points(int n, const std::vector<mod_t> &pts, mod_t m,
ConvolveCyclicFuncType f) {
assert(n <= mod_t::get_mod());
assert(pts.size() <= mod_t::get_mod());
if (n == 0) return {};
using u64 = std::uint64_t;
u64 m_64 = u64(m), k = pts.size(), nm1 = u64(m + mod_t(n - 1));
if (m_64 < k) { // f(0), …, f(m), …, f(k-1)
if (m_64 + n <= k) { // f(0), …, f(m), …, f(n+m-1), …, f(k-1)
return std::vector<mod_t>(pts.begin() + m_64, pts.begin() + m_64 + n);
} else if (nm1 < k) { // f(0), …, f(n+m-1), …, f(m), …, f(k-1), …, f(mod-1)
std::vector<mod_t> res;
res.reserve(n);
std::copy_n(pts.begin() + m_64, k - m_64, std::back_inserter(res));
std::copy_n(
internal::shift_sample_points_unsafe(mod_t::get_mod() - k, pts, mod_t(k), f).begin(),
mod_t::get_mod() - k, std::back_inserter(res));
std::copy_n(pts.begin(), nm1 + 1, std::back_inserter(res));
return res;
} else { // f(0), …, f(m), …, f(k-1), …, f(n+m-1)
std::vector<mod_t> res;
res.reserve(n);
std::copy_n(pts.begin() + m_64, k - m_64, std::back_inserter(res));
std::copy_n(internal::shift_sample_points_unsafe(m_64 + n - k, pts, mod_t(k), f).begin(),
m_64 + n - k, std::back_inserter(res));
return res;
}
} else { // f(0), …, f(k-1), …, f(m)
if (nm1 >= m_64) { // f(0), …, f(k-1), …, f(m), …, f(n+m-1), …, f(mod-1)
return internal::shift_sample_points_unsafe(n, pts, m, f);
} else if (nm1 < k) { // f(0), …, f(n+m-1), …, f(k-1), …, f(m), …, f(mod-1)
std::vector<mod_t> res;
res.reserve(n);
std::copy_n(internal::shift_sample_points_unsafe(int(-m), pts, m, f).begin(), int(-m),
std::back_inserter(res));
std::copy_n(pts.begin(), nm1 + 1, std::back_inserter(res));
return res;
} else { // f(0), …, f(k-1), …, f(n+m-1), …, f(m), …, f(mod-1)
std::vector<mod_t> res;
res.reserve(n);
std::copy_n(internal::shift_sample_points_unsafe(int(-m), pts, m, f).begin(), int(-m),
std::back_inserter(res));
std::copy_n(pts.begin(), k, std::back_inserter(res));
std::copy_n(internal::shift_sample_points_unsafe(nm1 - k + 1, pts, mod_t(k), f).begin(),
nm1 - k + 1, std::back_inserter(res));
return res;
}
}
}
template <typename mod_t, typename ConvolveCyclicFuncType>
std::vector<mod_t> shift_sample_points(const std::vector<mod_t> &pts, mod_t m,
ConvolveCyclicFuncType f) {
return shift_sample_points(pts.size(), pts, m, f);
}
} // namespace lib
#line 19 "math\\modulo\\factorial_modulo_prime.hpp"
namespace lib {
/**
* @brief 素数模数的阶乘计算
* @note 在 O(p^{1/2}\log p) 时间预处理，单次询问 O(p^{1/2})
* @tparam mod_t 素数模数
*/
template <typename mod_t>
class FactorialModPrime {
public:
using u32 = std::uint32_t;
using u64 = std::uint64_t;
template <typename ConvolveCyclicFuncType>
FactorialModPrime(ConvolveCyclicFuncType conv)
: v_(static_cast<u64>(std::sqrt(modint_traits<mod_t>::get_mod()))) {
const mod_t ONE(1);
mod_t mv = mod_t(v_), iv = ONE / mv;
block_prod_ = std::vector<mod_t>{ONE, mv + ONE};
block_prod_.reserve(v_ + 1);
u64 mask = UINT64_C(1) << 63;
while ((mask & v_) == 0) mask >>= 1;
for (u32 d = 1; d != v_;) {
std::vector<mod_t> g0(shift_sample_points(d, block_prod_, mod_t(d + 1), conv));
std::vector<mod_t> g1(shift_sample_points(d << 1 | 1, block_prod_, mod_t(d) * iv, conv));
std::copy(g0.begin(), g0.end(), std::back_inserter(block_prod_));
d <<= 1;
for (int i = 0; i <= d; ++i) block_prod_[i] *= g1[i];
if ((mask >>= 1) & v_) {
mod_t k(d | 1), dpv(v_ * k), prod(ONE);
for (int i = 1; i <= d; ++i) prod *= (dpv += ONE);
block_prod_.emplace_back(prod);
d |= 1;
for (int i = 0; i <= d; ++i) block_prod_[i] *= k, k += mv;
}
}
std::partial_sum(block_prod_.begin(), block_prod_.end(), block_prod_.begin(),
std::multiplies<>());
}
~FactorialModPrime() = default;
mod_t fac(u64 n) const {
if (n >= modint_traits<mod_t>::get_mod()) return mod_t(0);
mod_t res(1);
u64 pass = n / v_;
if (pass != 0) res *= block_prod_[pass - 1];
for (mod_t ONE(1), mpass(pass * v_), mn(n); mpass != mn;) res *= (mpass += ONE);
return res;
}
private:
const u64 v_;
std::vector<mod_t> block_prod_;
};
} // namespace lib
#line 8 "remote_test\\yuki\\math\\502.0.test.cpp"
int main() {
#ifdef LOCAL
std::freopen("in", "r", stdin), std::freopen("out", "w", stdout);
#endif
std::ios::sync_with_stdio(false);
std::cin.tie(0);
using mint = lib::MontModInt<1000000007>;
long long v;
std::cin >> v;
std::cout << lib::FactorialModPrime<mint>([](const std::vector<mint> &x,
const std::vector<mint> &y, int c) {
return lib::convolve_cyclic_mod(x, y, mint::get_mod(), c);
}).fac(v)
<< '\n';
return 0;
}