結果
問題 | No.215 素数サイコロと合成数サイコロ (3-Hard) |
ユーザー |
![]() |
提出日時 | 2015-05-23 15:03:43 |
言語 | C++11 (gcc 13.3.0) |
結果 |
AC
|
実行時間 | 1,233 ms / 4,000 ms |
コード長 | 12,938 bytes |
コンパイル時間 | 949 ms |
コンパイル使用メモリ | 96,456 KB |
実行使用メモリ | 8,836 KB |
最終ジャッジ日時 | 2024-07-06 06:20:11 |
合計ジャッジ時間 | 5,081 ms |
ジャッジサーバーID (参考情報) |
judge3 / judge5 |
(要ログイン)
ファイルパターン | 結果 |
---|---|
other | AC * 2 |
ソースコード
#include <cstdio>#include <cmath>#include <cstring>#include <cstdlib>#include <ctime>#include <cassert>#include <iostream>#include <utility>#include <algorithm>#include <queue>#include <functional>#include <vector>#include <map>#include <set>#include <complex>#define getchar getchar_unlocked#define putchar putchar_unlockedusing namespace std;typedef long long int64;typedef long long unsigned uint64;typedef long double float80;typedef unsigned short uint16;typedef unsigned uint;typedef unsigned char uint8;uint get_uint() {uint n;int c;while( (c = getchar()) < '0') {;}n = c - '0';while( (c = getchar()) >= '0') {n = n * 10 + (c - '0');}return n;}void put_uint(uint n) {uint8 stack[30];int top = 0;do {stack[top++] = n % 10 + '0';n /= 10;} while(n != 0);while(top > 0) {putchar(stack[--top]);}putchar('\n');}// mod < 2^63template <uint64 mod, uint64 z, uint period, uint64 n_prime, uint64 r2>class Mod64 {private:typedef __uint128_t uint128;public:Mod64() {}Mod64(uint64 v) {n_ = montgomery_init(v);}static Mod64 root2pow(uint i) {return Mod64(z).pow(1ull << (period - i));}Mod64 operator+ (Mod64 rhs) const {Mod64 ret;uint64 t = this->n_ + rhs.n_;ret.set_value(t >= mod ? t - mod : t);return ret;}Mod64 operator- (Mod64 rhs) const {Mod64 ret;uint64 t = this->n_ - rhs.n_;ret.set_value(int64(t) < 0 ? t + mod : t);return ret;}Mod64 operator* (Mod64 rhs) const {Mod64 ret;uint64 t = montgomery_reduction(uint128(this->n_) * rhs.n_);ret.set_value(t);return ret;}Mod64 operator+= (Mod64 rhs) {return *this = *this + rhs;}Mod64 operator-= (Mod64 rhs) {return *this = *this - rhs;}Mod64 operator*= (Mod64 rhs) {return *this = *this * rhs;}uint64 get_raw_value() const {return this->n_;}uint64 get_value() const {return montgomery_reduction(this->n_);}void set_value(uint64 val) {this->n_ = val;}Mod64 pow(uint64 exp) const {Mod64 base = *this;Mod64 ret = Mod64(1);while (exp) {if (exp & 1) {ret *= base;}exp >>= 1;base *= base;}return ret;}Mod64 inverse() const {return pow(mod - 2);}private:uint64 n_;static uint64 montgomery_init(uint64 w) {return montgomery_reduction(uint128(w) * r2);}static uint64 montgomery_reduction(const uint128 w) {uint64 x = uint64(w) * n_prime;uint128 y = uint128(x) * mod + w;uint64 ret = y >> 64;if (ret >= mod) {ret -= mod;}return ret;}};template <typename mod_t>class NTT {public:static void auto_convolute(mod_t* poly_a, uint size_a) {uint size = size_a;uint ntt_size = 1;uint ldn = 0;while (ntt_size < size) {ntt_size <<= 1;ldn++;}ntt_size <<= 1;++ldn;fill(poly_a + size_a, poly_a + ntt_size, 0);ntt_dit4(poly_a, ldn, 1);for (uint i = 0; i < ntt_size; ++i) {poly_a[i] *= poly_a[i];}ntt_dit4(poly_a, ldn, -1);mod_t inv = mod_t(ntt_size).inverse();for (uint i = 0; i < 2 * size_a - 1; ++i) {poly_a[i] *= inv;}}static void convolute(mod_t* poly_a, uint size_a, mod_t* poly_b, uint size_b) {uint size = max(size_a, size_b);uint ntt_size = 1;uint ldn = 0;while (ntt_size < size) {ntt_size <<= 1;ldn++;}ntt_size <<= 1;++ldn;fill(poly_a + size_a, poly_a + ntt_size, 0);fill(poly_b + size_b, poly_b + ntt_size, 0);ntt_dit4(poly_a, ldn, 1);ntt_dit4(poly_b, ldn, 1);for (uint i = 0; i < ntt_size; ++i) {poly_a[i] *= poly_b[i];}ntt_dit4(poly_a, ldn, -1);mod_t inv = mod_t(ntt_size).inverse();for (uint i = 0; i < size_a + size_b - 1; ++i) {poly_a[i] *= inv;}}static void ntt_dit4(mod_t* f, uint ldn, int sign) {revbin_permute(f, 1u << ldn);ntt_dit4_core(f, ldn, sign);}private:static inline void sumdiff(mod_t& a, mod_t& b) {mod_t t = a - b;a += b;b = t;}static void revbin_permute(mod_t* A, uint n) {if (n <= 2) {return;}uint r = 0;uint nh = n >> 1;for (uint x = 1; x < n; ++x) {uint h = nh;while (! ((r ^= h) & h)) {h >>= 1;}if (r > x) {swap(A[x], A[r]);}}}static void ntt_dit4_core(mod_t *f, uint ldn, int sign) {const uint LX = 2;const uint n = 1u << ldn;if (ldn & 1) {for (uint i = 0; i < n; i += 2) {sumdiff(f[i], f[i+1]);}}mod_t imag = mod_t::root2pow(2);if (sign < 0) {imag = imag.inverse();}uint ldm = LX + (ldn & 1);mod_t one = mod_t(1);for (; ldm <= ldn; ldm += LX) {const uint m = 1u << ldm;const uint m4 = m >> LX;mod_t dw = mod_t::root2pow(ldm);if (sign < 0) {dw = dw.inverse();}mod_t w = one;mod_t w2 = w;mod_t w3 = w;for (uint j = 0; j < m4; ++j) {for (uint r = 0, i = j + r; r < n; r += m, i += m) {mod_t a0 = f[i + m4 * 0];mod_t a2 = f[i + m4 * 1] * w2;mod_t a1 = f[i + m4 * 2] * w;mod_t a3 = f[i + m4 * 3] * w3;mod_t t02 = a0 + a2;mod_t t13 = a1 + a3;f[i + m4 * 0] = t02 + t13;f[i + m4 * 2] = t02 - t13;t02 = a0 - a2;t13 = a1 - a3;t13 *= imag;f[i + m4 * 1] = t02 + t13;f[i + m4 * 3] = t02 - t13;}w *= dw;w2 = w * w;w3 = w * w2;}}}};const uint64 P1 = 0x3f91300000000001ull;const uint64 P2 = 0x3f93300000000001ull;const uint64 P1_INV_MOD_P2 = 0x19D3CB8000001FCAull;typedef Mod64<P1, 0x1941B388165C78EBull, 44,0x3f912fffffffffffull, 0x0298CD3E4612D42Aull> mod64_t1;typedef Mod64<P2, 0x394ba9c52fec1825ull, 44,0x3f932fffffffffffull, 0x30da2c2b74e27a1full> mod64_t2;// ------------------------------------------------------------------------------const uint mod = 1e9 + 7;const uint64 lim = (0xFFFFFFFFFFFFFFFFull / mod - mod) * mod;const uint MAX_PC = 300;const uint NTT_BUFF_SIZE = (1 << 14) + 1;uint dp[MAX_PC + 1][MAX_PC * 13 + 1];uint poly_P[MAX_PC * 13 + 1];uint poly_C[MAX_PC * 12 + 1];uint poly[MAX_PC * 25 + 1];uint poly_inv[MAX_PC * 25 + 100];uint sums[MAX_PC * 25 + 1];uint counts[MAX_PC * 25 + 1];uint rem[2 * MAX_PC * 25 + 100];uint64 temp[2 * MAX_PC * 25 + 100];mod64_t1 A1[NTT_BUFF_SIZE], B1[NTT_BUFF_SIZE];mod64_t2 A2[NTT_BUFF_SIZE], B2[NTT_BUFF_SIZE];void _poly_mul_restore(const uint size, mod64_t1* a1, mod64_t2* a2, uint* res) {for(uint i = 0; i < size; ++i) {uint64 x1 = a1[i].get_value();uint64 x2 = a2[i].get_value();if(x1 != x2) {uint64 dx = x2 - x1;if(int64(dx) < 0) {dx += P2; // P1 < P2}dx = (mod64_t2(dx) * mod64_t2(P1_INV_MOD_P2)).get_value();x1 = (x1 + dx % mod * (P1 % mod));}res[i] = x1 % mod;}}uint poly_square(const uint* p1, uint s1, uint* res) {for (uint i = 0; i < s1; ++i) {A1[i] = mod64_t1(p1[i]);A2[i] = mod64_t2(p1[i]);}NTT<mod64_t1>::auto_convolute(A1, s1);NTT<mod64_t2>::auto_convolute(A2, s1);_poly_mul_restore(s1 * 2 - 1, A1, A2, res);return s1 * 2 - 1;}uint poly_mul(const uint* p1, uint s1, const uint* p2, uint s2, uint* res) {for (uint i = 0; i < s1; ++i) {A1[i] = mod64_t1(p1[i]);A2[i] = mod64_t2(p1[i]);}for (uint i = 0; i < s2; ++i) {B1[i] = mod64_t1(p2[i]);B2[i] = mod64_t2(p2[i]);}NTT<mod64_t1>::convolute(A1, s1, B1, s2);NTT<mod64_t2>::convolute(A2, s1, B2, s2);_poly_mul_restore(s1 + s2 - 1, A1, A2, res);return s1 + s2 - 1;}uint _poly_rem(const uint* p1, uint s1, const uint* p2, uint s2, uint ofs, uint* ret) {if (s1 < s2) {copy(p1, p1 + s1, ret);return s1;}copy(p1, p1 + s1, temp);uint dif = s1 - s2 + 1;for (uint i = 0; i < dif; ++i) {uint64 coef = temp[i] % mod;if (!coef) {continue;}coef = mod - coef;for (uint j = ofs; j < s2; ++j) {temp[i + j] += coef * p2[j];if (temp[i + j] >= lim) {temp[i + j] -= lim;}}}for (uint i = dif; i < s1; ++i) {ret[i - dif] = temp[i] % mod;}return s2 - 1;}uint poly_tmp[2 * MAX_PC * 25 + 100];uint poly_tmp2[2 * MAX_PC * 25 + 100];uint poly_q[2 * MAX_PC * 25 + 100];uint poly_div(const uint* poly1, uint s1, const uint* , uint s2,const uint *inv, uint inv_size, uint *res) {assert(s1 >= s2);uint needed_size = s1 - s2 + 1;assert(inv_size >= needed_size);poly_mul(poly1, needed_size, inv, needed_size, res);return needed_size;}uint poly_rem(const uint* poly1, uint s1, const uint* poly2, uint s2, uint* inv, uint inv_size, uint* res) {if (s1 < s2) {copy(poly1, poly1 + s1, res);return s1;}uint q_size = poly_div(poly1, s1, poly2, s2, inv, inv_size, poly_q);poly_mul(poly_q, q_size, poly2, s2, poly_q);uint ofs = s1 - s2 + 1;for (uint i = 0; i < s2 - 1; ++i) {res[i] = (poly1[ofs + i] + mod - poly_q[ofs + i]) % mod;}return s2 - 1;}uint poly_inverse(const uint *poly, uint size, uint* inv, uint needed_size) {uint degs[100];uint deg = needed_size - 1;uint deg_pos = 0;while (deg) {degs[deg_pos++] = deg;deg >>= 1;}if (size < needed_size) {copy(poly, poly + size, poly_tmp);fill(poly_tmp + size, poly_tmp + needed_size, 0);} else {copy(poly, poly + needed_size, poly_tmp);}uint inv_size = 1;inv[0] = 1;while (deg_pos) {uint t = degs[--deg_pos];uint added = t + 1 - inv_size;poly_mul(poly_tmp, t + 1, inv, inv_size, poly_tmp2);poly_mul(poly_tmp2 + inv_size, added, inv, added, poly_tmp2);for (uint i = inv_size; i < inv_size + added; ++i) {inv[i] = (mod - poly_tmp2[i - inv_size] % mod) % mod;}inv_size += added;}return inv_size;}uint poly_power_rem(uint64 e, const uint* divisor, uint div_size, uint* ret) {uint64 mask = 1;while (mask <= e) {mask <<= 1;}mask >>= 1;uint ret_size = 0;ret[ret_size++] = 1;const uint inv_size = poly_inverse(divisor, div_size, poly_inv, div_size);while (mask) {if (e & mask) {ret[ret_size++] = 0;}mask >>= 1;if (mask == 0) {break;}ret_size = poly_square(ret, ret_size, ret);ret_size = poly_rem(ret, ret_size, divisor, div_size, poly_inv, inv_size, ret);}if (ret_size >= div_size) {ret_size = poly_rem(ret, ret_size, divisor, div_size, poly_inv, inv_size, ret);}return ret_size;}void conv(uint* poly, uint size) {poly[0] = 1;for (uint i = 1; i < size; ++i) {if (poly[i]) {poly[i] = mod - poly[i];}}}inline void add_mod(uint& a, uint b) {a += b;if (a >= mod) {a -= mod;}}void init_poly(const uint* dice, const uint T, uint* result) {for (uint i = 1; i <= T; ++i) {fill(dp[i], dp[i] + dice[5] * i + 1, 0);}dp[0][0] = 1;for (uint di = 0; di < 6; ++di) {const uint d = dice[di];for (uint t = 0; t < T; ++t) {for (uint i = t * dice[0]; i <= t * dice[di]; ++i) {if (dp[t][i]) {add_mod(dp[t + 1][i + d], dp[t][i]);}}}}for (uint i = 0; i <= dice[5] * T; ++i) {result[i] = dp[T][i];}}void solve() {const uint Ps[] = {2, 3, 5, 7, 11, 13};const uint Cs[] = {4, 6, 8, 9, 10, 12};uint64 N;uint P, C;while (~scanf("%llu %u %u", &N, &P, &C)) {init_poly(Ps, P, poly_P);init_poly(Cs, C, poly_C);uint size_p = 13 * P + 1;uint size_c = 12 * C + 1;poly_mul(poly_P, size_p, poly_C, size_c, poly);const uint poly_size = 13 * P + 12 * C + 1;sums[0] = 0;for (uint i = 1; i < poly_size; ++i) {sums[i] = (sums[i - 1] + poly[i]) % mod;}fill(counts, counts + poly_size, 0);counts[0] = 1;for (uint i = 0; i < poly_size; ++i) {if (counts[i]) {for (uint j = 0; i + j < poly_size; ++j) {add_mod(counts[i + j], uint64(counts[i]) * poly[j] % mod);}}}uint64 E = max(0ll, int64(N - poly_size + 1));conv(poly, poly_size);uint rem_size = poly_power_rem(E, poly, poly_size, rem);uint64 ans = 0;for (uint64 e = E; e < N; ++e) {uint total = 0;for (uint i = 0; i < rem_size; ++i) {total = (total + uint64(rem[rem_size - 1 - i]) * counts[i]) % mod;}ans = (ans + uint64(total) * (sums[poly_size - 1] - sums[N - e - 1] + mod)) % mod;rem[rem_size++] = 0;rem_size = _poly_rem(rem, rem_size, poly, poly_size, 2 * P + 4 * C, rem);}printf("%llu\n", ans);}}int main() {solve();return 0;}