#include #include #include #include #include #include #include #include #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& 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 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 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 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({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(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& cs) { function< poly(int, int) > rec = [&](int beg, int end) { if (end - beg == 1) { return poly(vector({1, cs[beg] % mod})); } int mid = (beg + end) / 2; return rec(beg, mid) * rec(mid, end); }; return rec(0, cs.size()); } static vector multipoint_evaluation(const poly& f, vector& points) { int s = points.size(); int tree_size = 4 << ilog2(s - 1); vector tree(tree_size); function< void(int, int, int) > rec = [&](int beg, int end, int k) { if (end - beg == 1) { tree[k] = poly(vector({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 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 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 temp; if (buff == nullptr) { temp = vector(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 tmp(vector(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 coefs; static R mod; static u64 lmod; static vector facts, ifacts, invs; static vector tmp64; }; R poly::mod; u64 poly::lmod; vector poly::facts, poly::ifacts, poly::invs; vector 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; }