問題 (difficulty = 2169)
atcoder.jp
解説
時間内に解けませんでした…
$M \le 10^{16}$ なので、$O(\sqrt{M})$ の素因数分解がギリギリできます。また、一般に素因数は少ないです。小さい順から $14$ 個の素数の積 $2 \times 3 \times 5 \times 7 \times 11 \times 13 \times 17 \times 19 \times 23 \times 29 \times 31 \times 37 \times 41 \times 43 = 1.3 \times 10^{16} > 10^{16} $ なので、$ M $ の素因数 $ p_k $ の個数 $K$ は高々 $13$ 個です。従ってそれぞれの $ p_k $ について、$ M $ が $ p_k$ を持つ個数 $ e_k $ *1 に、部分列の最小公倍数が、既になっているか否かの flag を持った bit DP が考えられます。
まず $A_i$ の内、$ M $ の約数でないものを部分列に入れた瞬間明らかに条件を満たさないので、それらは無視していいです。
それぞれの $A_i$ で $ p_k$ ごとに分解して考えるのがよさそうです。$ p_k$ について、$A_i$ の持つ個数が $ e_k $ 未満であるならば、それは $0$ 個であるのと変わりません。なぜならそうであるならば、$ e_k $ 個を持つ数をあとでどうせ部分列にいれなければならないからです。
例えば $ M = 36$ のとき、$ A = \{2, 9, 4, 12 \} $ でも $ A = \{1, 9, 4, 12 \} $ でも状況が変わらないことが分かります。$A_1 = 2$ の場合、素因数 $2$ を $1$ 個持っていますが、部分列の最小公倍数を $12$ (素因数 $2$ を $2$ 個持つ) にするには、どのみち後で素因数 $2$ を $2$ 個持つ $4$ か $12$ を部分列に入れなければなりません。従って、$A_1 = 1$ と状況が変わりません。すなわち、素因数が何個か、ではなく、素因数の個数が $ e_k $ 未満か、$ e_k $ に等しいか、だけが大事で*2それを $0, \ 1$ で管理すればよいです。
以上より、全ての $A_i$ を $2^K$ 通りの状態にまとめることができます。部分列の積の状態も $2^K$ 通りにまとめることができます。これで bit DP をすれば、この DP 部分の計算量は $O(2^K \times 2^K)$ つまり $O(4^K)$ です。$4^{13} = 6.7 \times 10^7$ なので間に合います。
ただ、このように $A_i$ をまとめて計算量を削減しなくても、 $O(N 2^K)$ で bit DP できますが、速い言語ならそれも意図せず通ってしまうようです。
遷移の詳細については、以下の実装例を見てください。
(1) $O(N 2^K)$ の bit DP (C++). 448 ms. これでもまあまあ余裕があります。
#include "bits/stdc++.h" #include <atcoder/all> using namespace atcoder; using namespace std; #define rep(i, n) for(long long i = 0; i < (n); ++i) #define rrep(i, n) for(long long i = (n) - 1; i >= 0; --i) #define PII pair<int, int> typedef long long ll; constexpr int MOD = 998244353; // 7 * 17 * 2^23 + 1 // 繰り返し2乗法 (非再帰) // T = modint でも動く。 template<typename T = long long> constexpr T my_pow(T N, long long a) { assert(0 <= a); T x = N, res = (T)1; while (a) { if (a & 1) res *= x; x *= x; // x は *this の (2のべき乗) 乗を管理する。 a >>= 1; } return res; } // 素因分解アルゴリズム (O(sqrt(N)). // T = long long (defalt) template<typename T = long long> vector<pair<T, T>> prime_factor(T N) { vector<pair<T, T>> res; T i = 2; while (i * i <= N) { if (N % i == 0) res.push_back({ i, 0 }); while (N % i == 0) { res.back().second++; N /= i; } i += 1 + (i % 2); //i == 2 の場合だけ +1, その他の場合は +2 } if (N > 1) { //sqrt((元の N)) より大きな素因数は高々1つしかない。 res.push_back({ N, 0 }); res.back().second++; } return res; } signed main() { ll N, M; cin >> N >> M; vector<ll> A(N); rep(i, N) cin >> A[i]; //A[i] を個数で管理。M を割り切らないものは無視。 map<ll, ll> cnt; rep(i, N) if (M % A[i] == 0) cnt[A[i]]++; //素因数が 0 個で壊れそうなので、M = 1 は別で処理。 if (M == 1) { cout << (my_pow((modint998244353)2, cnt[1]) - 1).val() << endl; return 0; } auto mp = prime_factor<ll>(M); int K = mp.size(); //(p_k)^(e_k) をあらかじめ計算しておく。 vector<ll> beki(K); rep(k, K) beki[k] = my_pow(mp[k].first, mp[k].second); vector<modint998244353> dp(1LL << K, 0); dp[0] = 1; for (auto&& iter = cnt.begin(); iter != cnt.end(); iter++) { modint998244353 T = my_pow((modint998244353)2, iter->second); int flag = 0; rep(k, K) { if (iter->first % beki[k] == 0) flag |= 1LL << k; } //ナップサック dp おなじみの、1行しか状態を持たないで、逆から遷移する工夫。 rrep(bit, 1LL << K) { //1個も選ばない場合は除くので、T - 1 dp[bit | flag] += dp[bit] * (T - 1); } } cout << dp[(1LL << K) - 1].val() << endl; }
(2) $O(4^K)$ の bit DP (C++). 193 ms. 確かに (1) より大分早くなっています。
#include "bits/stdc++.h" #include <atcoder/all> using namespace atcoder; using namespace std; #define rep(i, n) for(long long i = 0; i < (n); ++i) #define rrep(i, n) for(long long i = (n) - 1; i >= 0; --i) #define PII pair<int, int> typedef long long ll; constexpr int MOD = 998244353; // 7 * 17 * 2^23 + 1 // 繰り返し2乗法 (非再帰) // T = modint でも動く。 template<typename T = long long> constexpr T my_pow(T N, long long a) { assert(0 <= a); T x = N, res = (T)1; while (a) { if (a & 1) res *= x; x *= x; // x は *this の (2のべき乗) 乗を管理する。 a >>= 1; } return res; } // 素因分解アルゴリズム (O(sqrt(N)). // T = long long (defalt) template<typename T = long long> vector<pair<T, T>> prime_factor(T N) { vector<pair<T, T>> res; T i = 2; while (i * i <= N) { if (N % i == 0) res.push_back({ i, 0 }); while (N % i == 0) { res.back().second++; N /= i; } i += 1 + (i % 2); //i == 2 の場合だけ +1, その他の場合は +2 } if (N > 1) { //sqrt((元の N)) より大きな素因数は高々1つしかない。 res.push_back({ N, 0 }); res.back().second++; } return res; } signed main() { ll N, M; cin >> N >> M; vector<ll> A(N); rep(i, N) cin >> A[i]; //A[i] を個数で管理。M を割り切らないものは無視。 map<ll, ll> cnt; rep(i, N) if (M % A[i] == 0) cnt[A[i]]++; //素因数が 0 個で壊れそうなので、M = 1 は別で処理。 if (M == 1) { cout << (my_pow((modint998244353)2, cnt[1]) - 1).val() << endl; return 0; } auto mp = prime_factor<ll>(M); int K = mp.size(); //(p_k)^(e_k) をあらかじめ計算しておく。 vector<ll> beki(K); rep(k, K) beki[k] = my_pow(mp[k].first, mp[k].second); //2^K の状態の個数を数える。 vector<ll> state(1LL << K, 0); for (auto&& iter = cnt.begin(); iter != cnt.end(); iter++) { int flag = 0; rep(k, K) { if (iter->first % beki[k] == 0) flag |= 1LL << k; } state[flag] += iter->second; } vector<modint998244353> dp(1LL << K, 0); dp[0] = 1; rep(bit1, 1LL << K) { modint998244353 T = my_pow((modint998244353)2, state[bit1]); //ナップサック dp おなじみの、1行しか状態を持たないで、逆から遷移する工夫。 rrep(bit2, 1LL << K) { //1個も選ばない場合は除くので、T - 1 dp[bit2 | bit1] += dp[bit2] * (T - 1); } } cout << dp[(1LL << K) - 1].val() << endl; }
