#94409 (C++11) No.215 素数サイコロと合成数サイコロ (3-Hard)

提出ソース
結果

問題	No.215 素数サイコロと合成数サイコロ (3-Hard)
ユーザー	Min_25
提出日時	2016-05-28 09:35:36
言語	C++11 (gcc 11.4.0)
結果	AC
実行時間	470 ms / 4,000 ms
コード長	9,046 bytes
コンパイル時間	1,290 ms
コンパイル使用メモリ	89,440 KB
実行使用メモリ	7,968 KB
最終ジャッジ日時	2024-04-16 18:03:17
合計ジャッジ時間	3,005 ms
ジャッジサーバーID （参考情報）	judge1 / judge5
このコードへのチャレンジ
（要ログイン）
テストケース

テストケース表示
入力	結果	実行時間実行使用メモリ
testcase_00	AC	470 ms 7,808 KB
testcase_01	AC	470 ms 7,968 KB
権限があれば一括ダウンロードができます
ソースコード

raw source code
#include <cstdio>
#include <cassert>
#include <cmath>
#include <cstring>

#include <algorithm>
#include <iostream>
#include <vector>
#include <functional>

#define _fetch(_1, _2, _3, _4, name, ...) name
#define rep2(i, n) rep3(i, 0, n)
#define rep3(i, a, b) rep4(i, a, b, 1)
#define rep4(i, a, b, c) for (int i = int(a); i < int(b); i += int(c))
#define rep(...) _fetch(__VA_ARGS__, rep4, rep3, rep2, _)(__VA_ARGS__)

using namespace std;

using i64 = long long;
using u32 = unsigned;
using u64 = unsigned long long;
using f80 = long double;

namespace ntt {

template <u64 mod, u64 prim_root>
class Mod64 {
private:
  using u128 = __uint128_t;
  static constexpr u64 mul_inv(u64 n, int e=6, u64 x=1) {
    return e == 0 ? x : mul_inv(n, e-1, x*(2-x*n));
  }
public:
  static constexpr u64 inv = mul_inv(mod);
  static constexpr u64 r2 = -u128(mod) % mod;
  static constexpr int level = __builtin_ctzll(mod - 1);
  static_assert(inv * mod == 1, "invalid 1/M modulo 2^64.");

  Mod64() {}
  Mod64(u64 n) : x(init(n)) {};
  static u64 modulo() { return mod; }
  static u64 init(u64 w) { return reduce(u128(w) * r2); }
  static u64 reduce(const u128 w) { return u64(w >> 64) + mod - ((u128(u64(w) * inv) * mod) >> 64); }
  static Mod64 omega() { return Mod64(prim_root).pow((mod - 1) >> level); }
  Mod64& operator += (Mod64 rhs) { this->x += rhs.x; return *this; }
  Mod64& operator -= (Mod64 rhs) { this->x += 2 * mod - rhs.x; return *this; }
  Mod64& operator *= (Mod64 rhs) { this->x = reduce(u128(this->x) * rhs.x); return *this; }
  Mod64 operator + (Mod64 rhs) const { return Mod64(*this) += rhs; }
  Mod64 operator - (Mod64 rhs) const { return Mod64(*this) -= rhs; }
  Mod64 operator * (Mod64 rhs) const { return Mod64(*this) *= rhs; }
  u64 get() const { return reduce(this->x) % mod; }
  void set(u64 n) const { this->x = n; }
  Mod64 pow(u64 exp) const {
    Mod64 ret = Mod64(1);
    for (Mod64 base = *this; exp; exp >>= 1, base *= base) if (exp & 1) ret *= base;
    return ret;
  }
  Mod64 inverse() const { return pow(mod - 2); }
  friend ostream& operator << (ostream& os, const Mod64& m) { return os << m.get(); }
  static void debug() {
    printf("%llu %llu %llu %llu\n", mod, inv, r2, omega().get());
  }
  u64 x;
};

template <typename mod_t>
void convolute(mod_t* A, int s1, mod_t* B, int s2) {
  int s = s1 + s2 - 1;
  int size = 1;
  while (size < s) size <<= 1;

  mod_t roots[mod_t::level] = { mod_t::omega() };
  rep(i, 1, mod_t::level) roots[i] = roots[i - 1] * roots[i - 1];
  fill(A + s1, A + size, 0); ntt_dit4(A, size, 1, roots);
  fill(B + s2, B + size, 0); ntt_dit4(B, size, 1, roots);
  rep(i, size) A[i] *= B[i];
  ntt_dit4(A, size, -1, roots);
  mod_t inv = mod_t(size).inverse();
  rep(i, s) A[i] *= inv;
}

template <typename mod_t>
void rev_permute(mod_t* A, int n) {
  int r = 0, nh = n >> 1;
  rep(i, 1, n) {
    for (int h = nh; !((r ^= h) & h); h >>= 1);
    if (r > i) swap(A[i], A[r]);
  }
}

template <typename mod_t>
void ntt_dit4(mod_t* A, int n, int sign, mod_t* roots) {
  rev_permute(A, n);
  int logn = __builtin_ctz(n);
  if (logn & 1) rep(i, 0, n, 2) {
    mod_t a = A[i], b = A[i + 1];
    A[i] = a + b; A[i + 1] = a - b;
  }
  mod_t imag = roots[mod_t::level - 2];
  if (sign < 0) imag = imag.inverse();

  mod_t one = mod_t(1);
  rep(e, 2 + (logn & 1), logn + 1, 2) {
    const int m = 1 << e;
    const int m4 = m >> 2;
    mod_t dw = roots[mod_t::level - e];
    if (sign < 0) dw = dw.inverse();
    mod_t w = one, w2 = one, w3 = one;
    rep(j, m4) {
      rep(i, j, n, m) {
        mod_t a0 = A[i + m4 * 0] * one, a2 = A[i + m4 * 1] * w2;
        mod_t a1 = A[i + m4 * 2] * w,   a3 = A[i + m4 * 3] * w3;
        mod_t t02 = a0 + a2, t13 = a1 + a3;
        A[i + m4 * 0] = t02 + t13; A[i + m4 * 2] = t02 - t13;
        t02 = a0 - a2, t13 = (a1 - a3) * imag;
        A[i + m4 * 1] = t02 + t13; A[i + m4 * 3] = t02 - t13;
      }
      w *= dw; w2 = w * w; w3 = w2 * w;
    }
  }
}

const int size = 1 << 14;
using m64_1 = ntt::Mod64<34703335751681, 3>;
using m64_2 = ntt::Mod64<35012573396993, 3>;
m64_1 f1[size], g1[size];
m64_2 f2[size], g2[size];

} // namespace ntt

using R = u32;
class poly {
public:
  poly() {}
  poly(int n) : coefs(n) {}
  poly(int n, int c) : coefs(n, c % mod) {}
  poly(const vector<R>& v) : coefs(v) {}
  poly(const poly& f, int beg, int end=-1) {
    if (end < 0) end = beg, beg = 0;
    resize(end - beg);
    rep(i, beg, end) if (i < f.size()) coefs[i - beg] = f[i];
  }
  static u32 ilog2(u64 n) {
    return 63 - __builtin_clzll(n);
  }
  int size() const { return coefs.size(); }
  void resize(int s) { coefs.resize(s); }
  void push_back(R c) { coefs.push_back(c); }

  const R* data() const { return coefs.data(); }
  R* data() { return coefs.data(); }
  const R& operator [] (int i) const { return coefs[i]; }
  R& operator [] (int i) { return coefs[i]; }

  static void add(R& a, R b) { if ((a += b) >= mod) a -= mod; }
  static void sub(R& a, R b) { if (int(a -= b) < 0) a += mod; }

  poly operator - () {
    poly ret = *this;
    rep(i, ret.size()) ret[i] = (ret[i] == 0 ? 0 : mod - ret[i]);
    return ret;
  }
  poly& operator += (const poly& rhs) {
    if (size() < rhs.size()) resize(rhs.size());
    rep(i, rhs.size()) add(coefs[i], rhs[i]);
    return *this;
  }
  poly operator + (const poly& rhs) const {
    return poly(*this) += rhs;
  }
  poly& operator -= (const poly& rhs) {
    if (size() < rhs.size()) resize(rhs.size());
    rep(i, rhs.size()) sub(coefs[i], rhs[i]);
    return *this;
  }
  poly operator - (const poly& rhs) const {
    return poly(*this) -= rhs;
  }

  poly operator * (const poly& rhs) const {
    return this->mul(rhs);
  }
  poly& operator *= (const poly& rhs) {
    return *this = *this * rhs;
  }

  struct fdiv {
    using u128 = __uint128_t;
    fdiv(u64 n) : m(n) {
      s = (n == 1) ? 0 : 127 - __builtin_clzll(n - 1);
      x = ((u128(1) << s) + n - 1) / n;
    }
    friend u64 operator / (u64 n, fdiv d) { return u128(n) * d.x >> d.s; }
    friend u64 operator % (u64 n, fdiv d) { return n - n / d * d.m; }
    u64 m, s, x;
  };

  poly mul(const poly& g) const {
    using namespace ntt;
    auto& f = *this;
    rep(i, f.size()) f1[i] = f[i];
    rep(i, g.size()) g1[i] = g[i];
    convolute(f1, f.size(), g1, g.size());
    rep(i, f.size()) f2[i] = f[i];
    rep(i, g.size()) g2[i] = g[i];
    convolute(f2, f.size(), g2, g.size());
    auto inv = m64_2(m64_1::modulo()).inverse();
    auto fmod = fdiv(mod);
    auto mod1 = m64_1::modulo() % fmod;
    poly ret(f.size() + g.size() - 1);
    rep(i, ret.size()) {
      u64 r1 = f1[i].get(), r2 = f2[i].get();
      ret[i] = (r1 + (m64_2(r2 + m64_2::modulo() - r1) * inv).get() % fmod * mod1) % fmod;
    }
    return ret;
  }

  poly mul_inv(int prec) const {
    poly ret(1, 1);
    for (int e = 1, ne; e < prec; e = ne) {
      ne = min(2 * e, prec);
      poly h = poly(ret, ne - e) * -poly(ret * poly(*this, ne), e, ne);
      rep(i, e, ne) ret.push_back(h[i - e]);
    }
    return ret;
  }

  pair<poly, poly> divmod_pre(const poly& b, const poly& inv) const {
    if (size() < b.size()) {
      return make_pair(poly(), poly(*this));
    }
    int sq = size() - b.size() + 1;
    assert(size() >= sq && inv.size() >= sq);
    poly q = poly(poly(*this, sq) * poly(inv, sq), sq);
    poly r = poly(*this - q * b, sq, size());
    return make_pair(q, r);
  }

  poly rem(const poly& f, const poly& inv) const {
    return divmod_pre(f, inv).second;
  }

  // x^e (mod f)
  static poly x_pow_mod(u64 e, const poly& f) {
    if (e == 0) return poly(1, 1);
    poly ret = poly(vector<R>({1, 0}));
    poly inv = f.mul_inv(f.size());
    ret = ret.rem(f, inv);
    u64 mask = (u64(1) << ilog2(e)) >> 1;
    while (mask) {
      ret *= ret;
      if (e & mask) ret.push_back(0);
      ret = ret.rem(f, inv);
      mask >>= 1;
    }
    return ret;
  }

  void print() const {
    printf("[");
    if (size()) {
      printf("%u", coefs[0]);
      rep(i, 1, size()) printf(" %u", coefs[i]);
    }
    puts("]");
  }

public:
  vector<R> coefs;
  static R mod;
};
R poly::mod;

u32 dp[301][3901];

poly init_poly(const u32* dice, u32 T) {
  rep(i, 0, T + 1) fill(dp[i], dp[i] + dice[5] * i + 1, 0);
  dp[0][0] = 1;
  rep(di, 6) rep(t, T) rep(i, t * dice[0], t * dice[di] + 1) if (dp[t][i]) {
    poly::add(dp[t + 1][i + dice[di]], dp[t][i]);
  }
  poly ret(dice[5] * T + 1);
  rep(i, dice[5] * T + 1) ret[i] = dp[T][i];
  return ret;
}

void solve() {
  poly::mod = 1e9 + 7;
  const u32 Ps[] = {2, 3, 5, 7, 11, 13};
  const u32 Cs[] = {4, 6, 8, 9, 10, 12};

  u64 N;
  u32 P, C;
  while (~scanf("%llu %u %u", &N, &P, &C)) {
    auto p1 = init_poly(Ps, P);
    auto p2 = init_poly(Cs, C);
    auto mod_f = -(p1 * p2);
    mod_f[0] = 1;
    auto r = poly::x_pow_mod(N + mod_f.size() - 2, mod_f);
    u64 ans = 0;
    rep(i, r.size()) ans += r[i];
    printf("%llu\n", ans % poly::mod);
  }
}

int main() {
  clock_t beg = clock();
  solve();
  clock_t end = clock();
  fprintf(stderr, "%.3f sec\n", double(end - beg) / CLOCKS_PER_SEC);
  return 0;
}
yukicoder

結果

テストケース

ソースコード
