結果

問題 No.215 素数サイコロと合成数サイコロ (3-Hard)
ユーザー k_zt0215k_zt0215
提出日時 2016-09-04 13:15:12
言語 C++11
(gcc 11.4.0)
結果
AC  
実行時間 1,776 ms / 4,000 ms
コード長 12,990 bytes
コンパイル時間 3,735 ms
コンパイル使用メモリ 149,128 KB
実行使用メモリ 8,108 KB
最終ジャッジ日時 2023-08-09 22:07:25
合計ジャッジ時間 8,333 ms
ジャッジサーバーID
(参考情報)
judge11 / judge13
このコードへのチャレンジ(β)

テストケース

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

ソースコード

diff #

#include <cstdio>
#include <cassert>
#include <cmath>
#include <ctime>

#include <iostream>
#include <vector>
#include <tuple>
#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 R = u32;
class poly {
  enum {
    KARATSUBA_CUTOFF = 64,
    DIV_CUTOFF = 128
  };
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);
  }
  static void init_mod(int s, R m) {
    mod = m;
    lmod = (u64(-1) / m - m) * m;

    facts.resize(s + 1, 1);
    ifacts.resize(s + 1, 1);
    invs.resize(s + 1, 1);
    rep(i, 2, s + 1) {
      invs[i] = u64(invs[mod % i]) * (mod - mod / i) % mod;
      facts[i] = u64(facts[i - 1]) * i % mod;
      ifacts[i] = u64(ifacts[i - 1]) * invs[i] % mod;
    }
  }

  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 add64(u64& a, u64 b) { if ((a += b) >= lmod) a -= lmod; }
  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;
  }

  // return a * b (mod x^prec)
  poly mul(const poly& b, int prec=-1) const {
    if (prec < 0) prec = max(0, size() + b.size() - 1);
    poly ret = poly(prec);
    amul(data(), size(), b.data(), b.size(), ret.data(), prec);
    return ret;
  }

  // return a / b (mod x^prec)
  poly rev_div(const poly& b, int prec) const {
    if (prec < 0) prec = size();
    poly q = poly(*this);
    q.resize(b.size() + prec - 1);
    return q.divmod(b).first;
  }

  // return 1 / b (mod x^prec)
  poly mul_inv(int prec) const {
    vector<int> precs;
    while (prec > 1) precs.push_back(prec), prec = (prec + 1) >> 1;

    poly ret(1, 1);
    for (int e = 1, ne; precs.size(); e = ne) {
      ne = precs.back(); precs.pop_back();
      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;
  }

  // return (q, r) such that a = q * b + r
  // - 2 * M(n/2) + 4 * M(n/4) + ...
  pair<poly, poly> divmod(const poly& b) const {
    if (size() < b.size()) {
      return make_pair(poly(), poly(*this));
    }
    poly q(size() - b.size() + 1);
    poly r(b.size() - 1);
    divmod_dc(data(), size(), b.data(), b.size(), q.data(), r.data());
    return make_pair(q, r);
  }

  // - 2 * M(n)
  pair<poly, poly> divmod_pre(const poly& b, const poly& inv) {
    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 {
    return divmod(f).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}));
    ret = ret.rem(f);
    u64 mask = (u64(1) << ilog2(e)) >> 1;
    while (mask) {
      ret *= ret;
      if (e & mask) ret.push_back(0);
      ret = ret.rem(f);
      mask >>= 1;
    }
    return ret;
  }

  // ----------------

  R evaluate(R x) const {
    R ret = 0;
    rep(i, size()) ret = (u64(ret) * x + coefs[i]) % mod;
    return ret;
  }

  static poly bernoullis(int N) {
    assert(int(ifacts.size()) > N + 1);
    poly ret = poly(vector<R>(ifacts.begin() + 1, ifacts.begin() + N + 2));
    ret = ret.mul_inv(ret.size());
    rep(i, ret.size()) ret[i] = u64(ret[i]) * facts[i] % mod;
    return ret;
  }

  static poly euler_numbers(int N) {
    assert(int(ifacts.size()) > N);
    poly ret = poly(N + 1);
    rep(i, 0, N + 1, 2) ret[i] = ifacts[i];
    ret = ret.mul_inv(ret.size());
    rep(i, ret.size()) ret[i] = u64(ret[i]) * facts[i] % mod;
    return ret;
  }

  static poly expand(vector<R>& cs) {
    function< poly(int, int) > rec = [&](int beg, int end) {
      if (end - beg == 1) {
        return poly(vector<R>({1, cs[beg] % mod}));
      } 
      int mid = (beg + end) / 2;
      return rec(beg, mid) * rec(mid, end);
    };
    return rec(0, cs.size());
  }

  static vector<R> multipoint_evaluation(const poly& f, vector<R>& points) {
    int s = points.size();
    int tree_size = 4 << ilog2(s - 1);

    vector<poly> tree(tree_size);
    function< void(int, int, int) > rec = [&](int beg, int end, int k) {
      if (end - beg == 1) {
        tree[k] = poly(vector<R>({1, (mod - points[beg] % mod) % mod}));
      } else {
        int mid = (beg + end) >> 1;
        rec(beg, mid, 2 * k + 1);
        rec(mid, end, 2 * k + 2);
        tree[k] = tree[2 * k + 1] * tree[2 * k + 2];
      }
    };
    rec(0, s, 0);

    vector<R> res(s);
    function< void(const poly&, int, int, int) > rec2 = [&](const poly& g, int beg, int end, int k) {
      auto r = g.rem(tree[k]);
      if (end - beg <= 32) {
        rep(i, beg, end) res[i] = r.evaluate(points[i]);
      } else {
        int mid = (beg + end) >> 1;
        rec2(r, beg, mid, 2 * k + 1);
        rec2(r, mid, end, 2 * k + 2);
      }
    };
    rec2(f, 0, s, 0);

    return res;
  }

  static R fact_mod(u32 N) {
    if (N >= mod) return 0;
    if (N <= 1) return 1 % mod;
    int v = sqrt(N);
    vector<R> cs(v);
    rep(i, v) cs[i] = (i * v + 1);
    auto f = expand(cs);

    rep(i, v) cs[i] = i;
    auto vs = multipoint_evaluation(f, cs);

    R ret = 1;
    rep(i, v) ret = u64(ret) * vs[i] % mod;
    rep(i, v * v + 1, N + 1) ret = u64(ret) * i % mod;
    return ret;
  }

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

private:
  // f * g
  static void amul(const R* a, int sa, const R* b, int sb, R* res, int prec, R* buff=nullptr) {
    if (sa < sb) return amul(b, sb, a, sa, res, prec, buff);
    if (prec < 0) prec = max(0, sa + sb - 1);
    if (sb < KARATSUBA_CUTOFF) {
      mul_basecase(a, sa, b, sb, res, prec);
    } else {
      // ...
      vector<R> temp;
      if (buff == nullptr) {
        temp = vector<R>(8 * sa + 100);
        buff = temp.data();
      }
      int q = sa / sb, r = sa % sb;
      if (r > 0 && q * pow(sa / float(q), 1.59) > (q + 1) * pow(sb, 1.59)) q += 1;
      int s = (sa + q - 1) / q;
      if (sb * q < sa) { 
        copy(b, b + sb, buff); fill(buff + sb, buff + s, 0); b = buff; buff += s; sb = s; 
      }
      if (sb * q > sa) {
        copy(a, a + sa, buff); fill(buff + sa, buff + sb * q, 0); a = buff; buff += sb * q;
      }
      fill(res, res + prec, 0);
      rep(i, q) {
        mul_karatsuba(a + i * sb, b, sb, buff, buff + 2 * sb - 1);
        rep(j, i * sb, min((i + 2) * sb - 1, prec)) add(res[j], buff[j - i * sb]);
      }
    }
  }

  static void mul_karatsuba(const R* a, const R* b, int s, R* res, R* buff) {
    if (s <= KARATSUBA_CUTOFF) {
      return mul_basecase(a, s, b, s, res, 2 * s - 1);
    }
    int sh = s / 2, sl = s - s / 2;
    
    mul_karatsuba(a, b, sl, res, buff);
    res[2 * sl - 1] = 0;
    mul_karatsuba(a + sl, b + sl, sh, res + 2 * sl, buff);

    auto* q1 = buff; copy(a, a + sl, q1); buff += sl;
    auto* q2 = buff; copy(b, b + sl, q2); buff += sl;
    auto* r1 = buff; buff += 2 * sl;

    rep(i, sh) add(q1[i], a[i + sl]);
    if (a != b) {
      rep(i, sh) add(q2[i], b[i + sl]);
    } else {
      q2 = q1;
    }
    mul_karatsuba(q1, q2, sl, r1, buff);

    rep(i, 2 * sl - 1) sub(r1[i], res[i]);
    rep(i, 2 * sh - 1) sub(r1[i], res[i + 2 * sl]);
    rep(i, 2 * sl - 1) add(res[i + sl], r1[i]);

    buff -= 4 * sl;
  }

  static void square_basecase(const R* a, int s, R* res, int prec=-1) {
    if (prec < 0) prec = max(0, 2 * s - 1);
    tmp64.assign(prec, 0);
    rep(i, s) tmp64[2 * i] = u64(a[i]) * a[i];
    rep(i, s) if (a[i]) {
      u32 c = (a[i] << 1) % mod;
      rep(j, i + 1, min(prec - i, s)) add64(tmp64[i + j], u64(c) * a[j]);
    }
    rep(i, prec) res[i] = tmp64[i] % mod;
  }

  static void mul_basecase(const R* a, int sa, const R* b, int sb, R* res, int prec=-1) {
    if (a == b) return square_basecase(a, sa, res, prec);
    if (prec < 0) prec = max(0, sa + sb - 1);
    tmp64.assign(prec, 0);
    rep(i, sb) if (b[i]) rep(j, min(prec - i, sa)) add64(tmp64[i + j], u64(b[i]) * a[j]);
    rep(i, prec) res[i] = tmp64[i] % mod;
  }

  // f / g (mod x^prec)
  static void rev_div_basecase(const R* a, int sa, const R* b, int sb, R* res, int prec) {
    assert(b[0] == 1);
    tmp64.assign(prec, 0);
    rep(i, min(prec, sa)) tmp64[i] = a[i];
    rep(i, prec) {
      R c = tmp64[i] % mod;
      if (c) rep(j, 1, min(prec - i, sb)) add64(tmp64[i + j], u64(mod - c) * b[j]);
      res[i] = c;
    }
  }

  // f % g
  static void divmod_dc32(R* a, int sa, const R* b, int sb, R* buff) {
    if (sa < sb) return;
    int d = sa - sb;
    divmod_dc21(a, 2 * d + 1, b, d + 1, buff);
    amul(a, d + 1, b + d + 1, sb - (d + 1), buff, sb - 1, buff + sb - 1);
    rep(i, sb - 1) sub(a[sa - 1 - i], buff[sb - 2 - i]);
  }

  static void divmod_dc21(R* a, int sa, const R* b, int sb, R* buff) {
    if (sb < DIV_CUTOFF || sa - sb < DIV_CUTOFF) {
      return divmod_basecase(a, sa, b, sb, a, a + sa - sb + 1);
    }
    int h = sb >> 1;
    if (sa - h >= sb) {
      divmod_dc32(a, sa - h, b, sb, buff);
      divmod_dc32(a + (sa - h) - (sb - 1), h + (sb - 1), b, sb, buff);
    } else {
      divmod_dc32(a, sa, b, sb, buff);
    }
  }

  static void divmod_dc(const R* a, int sa, const R* b, int sb, R* q, R* r) {
    assert(sa >= sb);
    int dq = sa / sb, dr = sa % sb;
    vector<R> tmp(vector<R>(a, a + sa)), buff(8 * sb + 100, 0);
    auto* t = tmp.data();
    rep(i, dq) {
      int end = dr + sb * (i + 1);
      int beg = max(0, end - (2 * sb - 1));
      divmod_dc21(t + beg, end - beg, b, sb, buff.data());
    }
    rep(i, sa - sb + 1) q[i] = t[i];
    rep(i, sb - 1) r[i] = t[i + sa - sb + 1];
  }

  static void divmod_basecase(const R* a, int sa, const R* b, int sb, R* q, R* r) {
    assert(sb >= 1 && b[0] == 1);
    tmp64.resize(sa);
    rep(i, sa) tmp64[i] = a[i];
    int d = sa - sb + 1;
    rep(i, d) {
      R c = tmp64[i] % mod;
      if (c) rep(j, 1, sb) add64(tmp64[i + j], u64(mod - c) * b[j]);
      q[i] = c;
    }
    rep(i, d, sa) r[i - d] = tmp64[i] % mod;
  }

public:
  vector<R> coefs;
  static R mod;
  static u64 lmod;
  static vector<R> facts, ifacts, invs;
  static vector<u64> tmp64;
};
R poly::mod;
u64 poly::lmod;
vector<R> poly::facts, poly::ifacts, poly::invs;
vector<u64> poly::tmp64;

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() {
  u32 mod = 1e9 + 7;
  poly::init_mod(0, mod);
  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;
}
0