結果
| 問題 | No.215 素数サイコロと合成数サイコロ (3-Hard) |
| ユーザー |
 Min_25
|
| 提出日時 | 2016-05-12 23:59:08 |
| 言語 | C++11 (gcc 11.4.0) |
| 結果 |
RE
|
| 実行時間 | - |
| コード長 | 12,932 bytes |
| コンパイル時間 | 1,016 ms |
| コンパイル使用メモリ | 91,324 KB |
| 実行使用メモリ | 7,616 KB |
| 最終ジャッジ日時 | 2024-10-05 14:07:48 |
| 合計ジャッジ時間 | 1,903 ms |
|
ジャッジサーバーID （参考情報） |
judge5 / judge1 |
（要ログイン）
ソースコード
#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 {
if (buff == nullptr) buff = vector<R>(8 * sa + 100).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(sa), buff(8 * sb + 100, 0);
copy(a, a + sa, tmp.data());
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;
}