結果

問題 No.8123 Calculated N !
ユーザー PNJ
提出日時 2025-04-01 22:46:59
言語 C++23
(gcc 13.3.0 + boost 1.87.0)
結果
AC  
実行時間 1,347 ms / 2,000 ms
コード長 12,986 bytes
コンパイル時間 4,194 ms
コンパイル使用メモリ 295,428 KB
実行使用メモリ 12,620 KB
最終ジャッジ日時 2025-04-01 22:47:19
合計ジャッジ時間 18,655 ms
ジャッジサーバーID
（参考情報） 		judge1 / judge3
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 6
other AC * 16
権限があれば一括ダウンロードができます

ソースコード

diff # 

#include <bits/stdc++.h>
using namespace std;

template <class T>
using vc = vector<T>;
template <class T>
using vvc = vector<vc<T>>;
template <class T>
using vvvc = vector<vvc<T>>;
template <class T>
using vvvvc = vector<vvvc<T>>;
template <class T>
using vvvvvc = vector<vvvvc<T>>;

#define vv(type, name, h, w) vector<vector<type>> name(h, vector<type>(w))
#define vvv(type, name, h, w, l) vector<vector<vector<type>>> name(h, vector<vector<type>>(w, vector<type>(l)))
#define vvvv(type, name, a, b, c, d) vector<vector<vector<vector<type>>>> name(a, vector<vector<vector<type>>>(b, vector<vector<type>>(c, vector<type>(d))))
#define vvvvv(type, name, a, b, c, d, e) vector<vector<vector<vector<vector<type>>>>> name(a, vector<vector<vector<vector<type>>>>(b, vector<vector<vector<type>>>(c, vector<vector<type>>(d, vector<type>(e)))))

#define elif else if

#define FOR1(a) for (long long _ = 0; _ < (long long)(a); _++)
#define FOR2(i, n) for (long long i = 0; i < (long long)(n); i++)
#define FOR3(i, l, r) for (long long i = l; i < (long long)(r); i++)
#define FOR4(i, l, r, c) for (long long i = l; i < (long long)(r); i += c)
#define FOR1_R(a) for (long long _ = (long long)(a) - 1; _ >= 0; _--)
#define FOR2_R(i, n) for (long long i = (long long)(n) - 1; i >= (long long)(0); i--)
#define FOR3_R(i, l, r) for (long long i = (long long)(r) - 1; i >= (long long)(l); i--)
#define FOR4_R(i, l, r, c) for (long long i = (long long)(r) - 1; i >= (long long)(l); i -= (c))
#define overload4(a, b, c, d, e, ...) e
#define FOR(...) overload4(__VA_ARGS__, FOR4, FOR3, FOR2, FOR1)(__VA_ARGS__)
#define FOR_R(...) overload4(__VA_ARGS__, FOR4_R, FOR3_R, FOR2_R, FOR1_R)(__VA_ARGS__)
#define FOR_in(a, A) for (auto a: A)
#define FOR_each(a, A) for (auto &&a: A)
#define FOR_subset(t, s) for(long long t = (s); t >= 0; t = (t == 0 ? -1 : (t - 1) & (s)))

#define all(x) x.begin(), x.end()
#define len(x) ll(x.size())

int popcount(int x) { return __builtin_popcount(x); }
int popcount(uint32_t x) { return __builtin_popcount(x); }
int popcount(long long x) { return __builtin_popcountll(x); }
int popcount(uint64_t x) { return __builtin_popcountll(x); }
// __builtin_clz(x)は最上位bitからいくつ0があるか.
int topbit(int x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); }
int topbit(uint32_t x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); }
int topbit(long long x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }
int topbit(uint64_t x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }

// 入力
void rd() {}
void rd(char &c) { cin >> c; }
void rd(string &s) { cin >> s; }
void rd(int &x) { cin >> x; }
void rd(uint32_t &x) { cin >> x; }
void rd(long long &x) { cin >> x; }
void rd(uint64_t &x) { cin >> x; }
template<class T>
void rd(vector<T> &v) {
  for (auto& x:v) rd(x);
}

void read() {}
template <class H, class... T>
void read(H &h, T &... t) {
  rd(h), read(t...);
}

#define CHAR(...) \
  char __VA_ARGS__; \
  read(__VA_ARGS__)

#define STRING(...) \
  string __VA_ARGS__; \
  read(__VA_ARGS__)

#define INT(...) \
  int __VA_ARGS__; \
  read(__VA_ARGS__)

#define U32(...) \
  uint32_t __VA_ARGS__; \
  read(__VA_ARGS__)

#define LL(...) \
  long long __VA_ARGS__; \
  read(__VA_ARGS__)

#define U64(...) \
  uint64_t __VA_ARGS__; \
  read(__VA_ARGS__)

#define VC(t, a, n) \
  vector<t> a(n); \
  read(a)

#define VVC(t, a, h, w) \
  vector<vector<t>> a(h, vector<t>(w)); \
  read(a)

//出力
void wt() {}
void wt(const char c) { cout << c; }
void wt(const string s) { cout << s; }
void wt(int x) { cout << x; }
void wt(uint32_t x) { cout << x; }
void wt(long long x) { cout << x; }
void wt(uint64_t x) { cout << x; }
template<class T>
void wt(const vector<T> v) {
  int n = v.size();
  for (int i = 0; i < n; i++) {
    if (i) wt(' ');
    wt(v[i]);
  }
}

void print() { wt('\n'); }
template <class Head, class... Tail>
void print(Head &&head, Tail &&... tail) {
  wt(head);
  if (sizeof...(Tail)) wt(' ');
  print(forward<Tail>(tail)...);
}

/////////////////////////////////////////////////////////////////////////////////////////
u_int64_t random_u64(u_int64_t l, u_int64_t r) {
  static std::random_device rd;
  static std::mt19937_64 gen(rd());
  std::uniform_int_distribution<u_int64_t> dist(l, r);
  return dist(gen);
}

long long pow_mod(long long a, long long r, long long mod) {
  long long res = 1, p = a % mod;
  while (r) {
    if ((r % 2) == 1) res = res * p % mod;
    p = p * p % mod, r >>= 1;
  }
  return res;
}

long long mod_inv(long long a, long long mod) {
  if (mod == 1) return 0;
  a %= mod;
  long long b = mod, s = 1, t = 0;
  while (1) {
    if (a == 1) return s;
    t -= (b / a) * s;
    b %= a;
    if (b == 1) return t + mod;
    s -= (a / b) * t;
    a %= b;
  }
}

long long Garner(vector<long long> Rem, vector<long long> Mod, int MOD) {
  assert (Rem.size() == Mod.size());
  long long mod = MOD;
  Rem.push_back(0);
  Mod.push_back(mod);
  long long n = Mod.size();
  vector<long long> coffs(n, 1);
  vector<long long> constants(n, 0);
  for (int i = 0; i < n - 1; i++) {
    long long v = (Mod[i] + Rem[i] - constants[i]) % Mod[i];
    v *= mod_inv(coffs[i], Mod[i]);
    v %= Mod[i];
    for (int j = i + 1; j < n; j++) {
      constants[j] = (constants[j] + coffs[j] * v) % Mod[j];
      coffs[j] = (coffs[j] * Mod[i]) % Mod[j];
    }
  }
  return constants[n - 1];
}

long long Tonelli_Shanks(long long a, long long mod) {
  a %= mod;
  if (a < 2) return a;
  if (pow_mod(a, (mod - 1) / 2, mod) != 1) return -1;
  if (mod % 4 == 3) return pow_mod(a, (mod + 1) / 4, mod);

  long long b = 3;
  if (mod != 998244353) {
    while (pow_mod(b, (mod - 1) / 2, mod) == 1) {
      b = random_u64(2, mod - 1);
    }
  }

  long long q = mod - 1;
  long long Q = 0;
  while (q % 2 == 0) {
    Q++, q /= 2;
  }

  long long x = pow_mod(a, (q + 1) / 2, mod);
  b = pow_mod(b, q, mod);

  long long shift = 2;
  while ((x * x) % mod != a) {
    long long error = (((pow_mod(a, mod - 2, mod) * x) % mod) * x) % mod;
    if (pow_mod(error, 1 << (Q - shift), mod) != 1) {
      x = (x * b) % mod;
    }
    b = (b * b) % mod;
    shift++;
  }
  return x;
}

template <int mod>
struct modint {
  static constexpr uint32_t umod = uint32_t(mod);
  static_assert(umod < (uint32_t(1) << 31));
  uint32_t val;

  static modint raw(uint32_t v) {
    modint x;
    x.val = v % umod;
    return x;
  }

  constexpr modint() : val(0) {}
  constexpr modint(uint32_t x) : val(x % umod) {}
  constexpr modint(uint64_t x) : val(x % umod) {}
  constexpr modint(unsigned __int128 x) : val(x % umod) {}
  constexpr modint(int x) : val((x %= umod) < 0 ? x + umod : x) {};
  constexpr modint(long long x) : val((x %= umod) < 0 ? x + umod : x) {};
  constexpr modint(__int128 x) : val((x %= umod) < 0 ? x + umod : x) {};

  bool operator<(const modint &other) const { return val < other.val; }
  modint &operator+=(const modint &p) {
    if ((val += p.val) >= umod) val -= umod;
    return *this;
  }
  modint &operator-=(const modint &p) {
    if ((val += umod - p.val) >= umod) val -= umod;
    return *this;
  }
  modint &operator*=(const modint &p) {
    val = uint64_t(val) * p.val % umod;
    return *this;
  }
  modint &operator/=(const modint &p) {
    val = uint64_t(val) * p.inverse().val % umod;
    return *this;
  }
  modint operator-() const { return modint::raw(val ? umod - val : uint32_t(0)); }
  modint operator+(const modint &p) const { return modint(*this) += p; }
  modint operator-(const modint &p) const { return modint(*this) -= p; }
  modint operator*(const modint &p) const { return modint(*this) *= p; }
  modint operator/(const modint &p) const { return modint(*this) /= p; }
  bool operator==(const modint &p) const { return val == p.val; }
  bool operator!=(const modint &p) const { return val != p.val; }

  modint inverse() const {
    int a = val, b = umod, s = 1, t = 0;
    while (1) {
      if (a == 1) return modint(s);
      t -= (b / a) * s;
      b %= a;
      if (b == 1) return modint(t + umod);
      s -= (a / b) * t;
      a %= b;
    }

  }

  modint pow(long long n) const {
    assert(n >= 0);
    modint res(1), a(val);
    while (n > 0) {
      if (n & 1) res *= a;
      a *= a;
      n >>= 1;
    }
    return res;
  }

  uint32_t get() const { return val; }

  static constexpr int get_mod() { return umod; }
  
  static constexpr pair<int, int> ntt_info() {
    if (mod == 167772161) return {25, 17};
    if (mod == 469762049) return {26, 30};
    if (mod == 754974721) return {24, 362};
    if (mod == 880803841) return {23, 211};
    if (mod == 998244353) return {23, 31};
    return {-1, -1};
  }
};

template <int mod>
void rd(modint<mod> &x) {
  uint32_t y;
  cin >> y;
  x = y;
}

template <int mod>
void wt(modint<mod> x) {
  wt(x.val);
}

// https://nyaannyaan.github.io/library/prime/prime-enumerate.hpp
vector<int> prime_enumerate(int N) {
  vector<bool> sieve(N / 3 + 1, 1);
  for (int p = 5, d = 4, i = 1, sqn = sqrt(N); p <= sqn; p += d = 6 - d, i++) {
    if (!sieve[i]) continue;
    for (int q = p * p / 3, r = d * p / 3 + (d * p % 3 == 2), s = 2 * p,
             qe = sieve.size();
         q < qe; q += r = s - r)
      sieve[q] = 0;
  }
  vector<int> ret{2, 3};
  for (int p = 5, d = 4, i = 1; p <= N; p += d = 6 - d, i++)
    if (sieve[i]) ret.push_back(p);
  while (!ret.empty() && ret.back() > N) ret.pop_back();
  return ret;
}


vector<int> prime;

// https://nyaannyaan.github.io/library/multiplicative-function/prime-counting-o2d3.hpp
inline int64_t my_div(int64_t n, int64_t p) { return double(n) / p; };

int64_t prime_counting(long long N) {
  if(N < 2) return 0;

  using i64 = long long;

  i64 N6, N3, N2, N23, nsz;
  vector<i64> ns, h;
  vector<int> bit;

  // calculate N^{1/2}, N^{1/3}, N{1/6}, N{2/3}
  N2 = sqrt(N);
  N3 = pow(N, 1.0 / 3.0);
  while (N3 * N3 * N3 > N) N3--;
  while ((N3 + 1) * (N3 + 1) * (N3 + 1) <= N) N3++;
  N6 = sqrt(N3);
  N23 = N / N3;

  // precalc prime sieve below N ^ {1/2}
  // auto prime = prime_enumerate(N2 + 1000);
  // index of prime
  int pidx = 0;
  // restore pi(p - 1)
  i64 pi = 0;

  // initialize ns
  ns.reserve(N2 * 2 + 2);
  ns.push_back(0);
  for (int i = 1; i <= N2; i++) ns.push_back(my_div(N, i));
  for (int i = ns.back() - 1; i > 0; i--) ns.push_back(i);
  nsz = ns.size();

  // initialize h
  copy(begin(ns), end(ns), back_inserter(h));
  for (auto &i : h) --i;

  // p <= N ^ {1/6}
  while (prime[pidx] <= N6) {
    int p = prime[pidx];
    i64 p2 = i64(p) * p;
    for (i64 i = 1, n = ns[i]; i < nsz && n >= p2; n = ns[++i])
      h[i] -= h[i * p <= N2 ? i * p : nsz - my_div(n, p)] - pi;
    ++pidx;
    ++pi;
  }

  // fenwick tree, which restore [ h(p, 1), h(p, N ^ {2/3}) )
  // bit[i] corresponds to h[i + N3] (1 <= i)
  bit.resize(nsz - N3);

  auto dfs = [&](auto rec, i64 cur, int pid, int flag) -> void {
    if (flag) {
      // if cur > N^{1/2} :
      // cur <= floor(N / id)
      // -> cur * id + n = N, n < id < cur
      // -> id <= floor(N / cur)
      int id = cur <= N2 ? nsz - cur : my_div(N, cur);
      // decrease h[N3] ~ h[id]
      if (id > N3)
        for (id -= N3; id; id -= id & -id) --bit[id];
    }
    for (int dst = pid; cur * prime[dst] < N23; dst++)
      rec(rec, cur * prime[dst], dst, true);
  };

  // N ^ {1/6} < p <= N ^ {1/3}
  while (prime[pidx] <= N3) {
    int p = prime[pidx];
    i64 p2 = i64(p) * p;
    // update N ^ {2/3} <= n <= N
    for (i64 i = 1; i <= N3; i++) {
      // ns[i] >= p2 > N^{1/3}
      if (p2 > ns[i]) break;
      // id of floor(N/ip)
      int id = i * p <= N2 ? i * p : nsz - my_div(ns[i], p);
      // current value of h[id]

      i64 sm = h[id];
      // if floor(N/ip) < N^{2/3}, add sum of fenwick tree
      if (id > N3)
        for (id -= N3; id < (int)bit.size(); id += id & -id) sm += bit[id];
      h[i] -= sm - pi;
    }

    // update N ^ {1/3} <= n < N ^ {2/3}, using dfs
    dfs(dfs, p, pidx, false);

    ++pidx;
    ++pi;
  }

  // reflect fenwick tree
  for (int i = (int)bit.size() - 1; i; i--) {
    int j = i + (i & -i);
    if (j < (int)bit.size()) bit[i] += bit[j];
  }
  for (int i = 1; i < (int)bit.size(); i++) h[i + N3] += bit[i];

  // N ^ {1/3} < p <= N ^ {1/2}
  while (prime[pidx] <= N2) {
    int p = prime[pidx];
    i64 p2 = i64(p) * p;
    for (i64 i = 1, n = ns[i]; i < nsz && n >= p2; n = ns[++i])
      h[i] -= h[i * p <= N2 ? i * p : nsz - my_div(n, p)] - pi;
    ++pidx;
    ++pi;
  }

  return h[1];
}

using mint = modint<1000000007>;

int main() {
  LL(N);
  prime = prime_enumerate(10000000);
  long long n = prime_counting(N);
  mint x = mint(1);
  mint ans = mint(1);
  FOR(d, 2, N + 1) {
    long long m = prime_counting(N / d);
    x += mint(1);
    ans *= x.pow(n - m);
    n = m;
    if (n < prime.size()) break;
  }
  FOR(i, n) {
    long long p = (long long)prime[i];
    long long q = p * p;
    mint r = mint(1) + mint(N / p) + mint((N / q));
    while (q * p <= N) {
      q *= p;
      long long n = N / q;
      r += mint(N / q);
    }
    ans *= r;
  }
  print(ans);
}
0