Search for a command to run...
Progressive hints first, then the full explanation and implementation when you're ready to cash out.
Review status
AI-generated and still unreviewed. Double-check the details before internalizing them.
Hints
Open only as much as you need to keep the solve alive.
Forget prime powers. For this problem, , , and only care about their set of distinct prime factors. So replace each by its squarefree prime set .
You need . Since , each has at most primes. That tiny bound is the whole door into the problem.
Use the identity
where is a Stirling number of the second kind. Here , so only matters, even though is huge.
For a fixed , counts the number of -element prime sets such that . So swap the summation: count, for every small prime set , how many pairs have .
For a fixed , the bad pairs are pairs where some prime from is missing from both numbers. Inclusion-exclusion over subsets : count elements whose prime set avoids . Because every actual number has at most primes and , all needed subset states can be generated sparsely; no giant bitmask over all primes. Big math, tiny implementation. Classic Codeforces nonsense.
We need compute
The first trap is thinking the product matters as an integer. It does not. Prime powers are dead weight.
For each value , define as the set of distinct primes dividing it. Then
So the task is really:
Useful bound
Since , a number cannot have many distinct prime factors. In fact,
while multiplying by already exceeds . So every has size at most .
Therefore every union has size at most .
That turns the scary into something manageable, because we only ever evaluate for .
The transform we use
There is a standard identity:
where is a Stirling number of the second kind.
Why is this useful? Because has a clean counting meaning: it is the number of -prime sets contained inside .
So:
where
Now sum over all pairs:
And now swap the inner counting:
So the problem becomes: for many small prime sets , count pairs whose union covers .
Counting pairs covering a fixed
For a fixed set , a pair fails to cover if at least one prime in is absent from both numbers.
For , let be the number of array elements whose prime set is disjoint from :
By inclusion-exclusion, the number of pairs whose union covers all of is:
This is the key formula. No pair-by-pair garbage. That would explode.
But aren't there too many prime sets?
If we tried all subsets of all primes up to , yes, instant tragedy.
But we only need sets that can appear inside . Since each has at most primes, every relevant is a subset of the union of two actual element-prime-sets.
The implementation builds a sparse map over relevant prime sets only.
Represent a prime set as a sorted vector of compressed prime ids. To store it in hash maps, encode it into a compact integer key using fixed-width chunks. There are fewer than relevant prime ids, so bits per prime is enough. With at most primes, this fits in __int128.
Precomputing sparse data
For every array element set , enumerate all subsets and add freqSub[U]++.
Then for any query set , the number of elements intersecting can be found by inclusion-exclusion over subsets of :
So
Because , this costs at most map lookups. Totally fine.
Which do we iterate?
We generate all relevant by taking every pair of distinct element prime-set types, forming their union, and enumerating its subsets.
There is one important optimization: group equal prime sets first. The number of different squarefree kernels up to is manageable, and each set is tiny.
For every pair of types , enumerate all subsets of and insert them into a hash set of relevant keys. Since , each pair contributes at most subsets.
That sounds chunky, but with the actual limit and tiny prime-factor sets, it fits. The math is doing the heavy lifting; the implementation just needs to not be goofy.
Computing the coefficients
We need
For small , use the explicit formula:
Compute by fast exponentiation modulo .
Note that because and .
Final algorithm
For each test case:
freqSub[X]: how many elements have all primes of .cover(T) using inclusion-exclusion over subsets ;The empty set contributes nothing because , so it can be ignored.
Complexity
Each number has at most distinct prime factors, and every union has at most .
The algorithm is sparse over actually reachable prime sets. The heaviest parts are subset enumeration over sets of size at most and hash-map lookups. With total , this is designed for the constraints.
The main conceptual trick is the Stirling transform. It converts the annoying nonlinear value into counting covered prime subsets, where inclusion-exclusion can slap the problem into shape.
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
void setIO() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
}
const int MOD = 998244353;
const int SHIFT = 18;
const int MAXR = 12;
using u128 = __uint128_t;
struct U128Hash {
size_t operator()(u128 x) const {
unsigned long long lo = (unsigned long long)x;
unsigned long long hi = (unsigned long long)(x >> 64);
lo ^= lo >> 33;
lo *= 0xff51afd7ed558ccdULL;
lo ^= lo >> 33;
hi ^= hi >> 33;
hi *= 0xc4ceb9fe1a85ec53ULL;
hi ^= hi >> 33;
return (size_t)(lo ^ (hi + 0x9e3779b97f4a7c15ULL + (lo << 6) + (lo >> 2)));
}
};
ll modpow(ll a, ll e) {
ll r = 1;
while (e) {
if (e & 1) r = r * a % MOD;
a = a * a % MOD;
e >>= 1;
}
return r;
}
int addmod(int a, int b) {
int s = a + b;
if (s >= MOD) s -= MOD;
return s;
}
int submod(int a, int b) {
int s = a - b;
if (s < 0) s += MOD;
return s;
}
u128 encode_vec(const vector<int>& v) {
u128 key = 0;
for (int i = 0; i < (int)v.size(); i++) {
key |= (u128)(v[i] + 1) << (SHIFT * i);
}
return key;
}
int key_size(u128 key) {
int r = 0;
while (key) {
r++;
key >>= SHIFT;
}
return r;
}
vector<int> decode_key(u128 key) {
vector<int> v;
while (key) {
v.push_back((int)(key & ((1u << SHIFT) - 1)) - 1);
key >>= SHIFT;
}
return v;
}
void gen_subsets_from_vec(const vector<int>& v, vector<u128>& out) {
int m = (int)v.size();
out.clear();
out.reserve(1 << m);
for (int mask = 0; mask < (1 << m); mask++) {
u128 key = 0;
int pos = 0;
for (int i = 0; i < m; i++) {
if (mask >> i & 1) {
key |= (u128)(v[i] + 1) << (SHIFT * pos++);
}
}
out.push_back(key);
}
}
vector<int> unite_vecs(const vector<int>& a, const vector<int>& b) {
vector<int> c;
c.reserve(a.size() + b.size());
int i = 0, j = 0;
while (i < (int)a.size() || j < (int)b.size()) {
if (j == (int)b.size() || (i < (int)a.size() && a[i] < b[j])) c.push_back(a[i++]);
else if (i == (int)a.size() || b[j] < a[i]) c.push_back(b[j++]);
else {
c.push_back(a[i]);
i++, j++;
}
}
return c;
}
int main() {
setIO();
const int LIM = 200000;
vector<int> spf(LIM + 1);
for (int i = 2; i <= LIM; i++) {
if (!spf[i]) {
spf[i] = i;
if ((ll)i * i <= LIM) {
for (ll j = (ll)i * i; j <= LIM; j += i) {
if (!spf[j]) spf[j] = i;
}
}
}
}
spf[1] = 1;
int C[13][13]{};
for (int i = 0; i <= MAXR; i++) {
C[i][0] = C[i][i] = 1;
for (int j = 1; j < i; j++) C[i][j] = addmod(C[i - 1][j - 1], C[i - 1][j]);
}
int tc;
cin >> tc;
while (tc--) {
int n;
ll k;
cin >> n >> k;
vector<int> coef(MAXR + 1);
for (int r = 1; r <= MAXR; r++) {
int cur = 0;
for (int q = 0; q <= r; q++) {
int term = (ll)C[r][q] * modpow(q, k) % MOD;
if ((r - q) & 1) cur = submod(cur, term);
else cur = addmod(cur, term);
}
coef[r] = cur;
}
unordered_map<int, int> prime_id;
prime_id.reserve(n * 2 + 10);
int pcnt = 0;
unordered_map<u128, int, U128Hash> type_count;
type_count.reserve(n * 2 + 10);
vector<int> tmp;
vector<u128> subs;
for (int i = 0; i < n; i++) {
int x;
cin >> x;
tmp.clear();
while (x > 1) {
int p = spf[x];
auto it = prime_id.find(p);
if (it == prime_id.end()) it = prime_id.emplace(p, pcnt++).first;
tmp.push_back(it->second);
while (x % p == 0) x /= p;
}
sort(tmp.begin(), tmp.end());
type_count[encode_vec(tmp)]++;
}
vector<pair<u128, int>> types;
types.reserve(type_count.size());
for (auto &p : type_count) types.push_back(p);
unordered_map<u128, int, U128Hash> sup_count;
sup_count.reserve(n * 16 + 10);
vector<vector<int>> type_vecs;
type_vecs.reserve(types.size());
for (auto &[key, cnt] : types) {
vector<int> v = decode_key(key);
type_vecs.push_back(v);
gen_subsets_from_vec(v, subs);
for (u128 s : subs) sup_count[s] += cnt;
}
unordered_set<u128, U128Hash> relevant;
relevant.reserve(types.size() * 64 + 10);
vector<int> uni;
vector<u128> usubs;
int mtypes = (int)type_vecs.size();
for (int i = 0; i < mtypes; i++) {
for (int j = i; j < mtypes; j++) {
uni = unite_vecs(type_vecs[i], type_vecs[j]);
gen_subsets_from_vec(uni, usubs);
for (u128 s : usubs) {
if (s) relevant.insert(s);
}
}
}
unordered_map<u128, int, U128Hash> avoid_memo;
avoid_memo.reserve(relevant.size() * 2 + 10);
auto avoid_count = [&](u128 key) -> int {
auto found = avoid_memo.find(key);
if (found != avoid_memo.end()) return found->second;
vector<int> v = decode_key(key);
gen_subsets_from_vec(v, subs);
int hit = 0;
for (u128 s : subs) {
if (!s) continue;
auto it = sup_count.find(s);
if (it == sup_count.end()) continue;
if (key_size(s) & 1) hit += it->second;
else hit -= it->second;
}
int res = n - hit;
avoid_memo[key] = res;
return res;
};
ll ans = 0;
for (u128 tkey : relevant) {
vector<int> tv = decode_key(tkey);
int r = (int)tv.size();
if (coef[r] == 0) continue;
gen_subsets_from_vec(tv, usubs);
ll cover = 0;
for (u128 u : usubs) {
int c = avoid_count(u);
ll pairs = (ll)c * (c - 1) / 2 % MOD;
if (key_size(u) & 1) cover -= pairs;
else cover += pairs;
}
cover %= MOD;
if (cover < 0) cover += MOD;
ans = (ans + cover * coef[r]) % MOD;
}
cout << ans << '\n';
}
}#include <bits/stdc++.h>
using namespace std;
using ll = long long;
void setIO() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
}
const int MOD = 998244353;
const int SHIFT = 18;
const int MAXR = 12;
using u128 = __uint128_t;
struct U128Hash {
size_t operator()(u128 x) const {
unsigned long long lo = (unsigned long long)x;
unsigned long long hi = (unsigned long long)(x >> 64);
lo ^= lo >> 33;
lo *= 0xff51afd7ed558ccdULL;
lo ^= lo >> 33;
hi ^= hi >> 33;
hi *= 0xc4ceb9fe1a85ec53ULL;
hi ^= hi >> 33;
return (size_t)(lo ^ (hi + 0x9e3779b97f4a7c15ULL + (lo << 6) + (lo >> 2)));
}
};
ll modpow(ll a, ll e) {
ll r = 1;
while (e) {
if (e & 1) r = r * a % MOD;
a = a * a % MOD;
e >>= 1;
}
return r;
}
int addmod(int a, int b) {
int s = a + b;
if (s >= MOD) s -= MOD;
return s;
}
int submod(int a, int b) {
int s = a - b;
if (s < 0) s += MOD;
return s;
}
u128 encode_vec(const vector<int>& v) {
u128 key = 0;
for (int i = 0; i < (int)v.size(); i++) {
key |= (u128)(v[i] + 1) << (SHIFT * i);
}
return key;
}
int key_size(u128 key) {
int r = 0;
while (key) {
r++;
key >>= SHIFT;
}
return r;
}
vector<int> decode_key(u128 key) {
vector<int> v;
while (key) {
v.push_back((int)(key & ((1u << SHIFT) - 1)) - 1);
key >>= SHIFT;
}
return v;
}
void gen_subsets_from_vec(const vector<int>& v, vector<u128>& out) {
int m = (int)v.size();
out.clear();
out.reserve(1 << m);
for (int mask = 0; mask < (1 << m); mask++) {
u128 key = 0;
int pos = 0;
for (int i = 0; i < m; i++) {
if (mask >> i & 1) {
key |= (u128)(v[i] + 1) << (SHIFT * pos++);
}
}
out.push_back(key);
}
}
vector<int> unite_vecs(const vector<int>& a, const vector<int>& b) {
vector<int> c;
c.reserve(a.size() + b.size());
int i = 0, j = 0;
while (i < (int)a.size() || j < (int)b.size()) {
if (j == (int)b.size() || (i < (int)a.size() && a[i] < b[j])) c.push_back(a[i++]);
else if (i == (int)a.size() || b[j] < a[i]) c.push_back(b[j++]);
else {
c.push_back(a[i]);
i++, j++;
}
}
return c;
}
int main() {
setIO();
const int LIM = 200000;
vector<int> spf(LIM + 1);
for (int i = 2; i <= LIM; i++) {
if (!spf[i]) {
spf[i] = i;
if ((ll)i * i <= LIM) {
for (ll j = (ll)i * i; j <= LIM; j += i) {
if (!spf[j]) spf[j] = i;
}
}
}
}
spf[1] = 1;
int C[13][13]{};
for (int i = 0; i <= MAXR; i++) {
C[i][0] = C[i][i] = 1;
for (int j = 1; j < i; j++) C[i][j] = addmod(C[i - 1][j - 1], C[i - 1][j]);
}
int tc;
cin >> tc;
while (tc--) {
int n;
ll k;
cin >> n >> k;
vector<int> coef(MAXR + 1);
for (int r = 1; r <= MAXR; r++) {
int cur = 0;
for (int q = 0; q <= r; q++) {
int term = (ll)C[r][q] * modpow(q, k) % MOD;
if ((r - q) & 1) cur = submod(cur, term);
else cur = addmod(cur, term);
}
coef[r] = cur;
}
unordered_map<int, int> prime_id;
prime_id.reserve(n * 2 + 10);
int pcnt = 0;
unordered_map<u128, int, U128Hash> type_count;
type_count.reserve(n * 2 + 10);
vector<int> tmp;
vector<u128> subs;
for (int i = 0; i < n; i++) {
int x;
cin >> x;
tmp.clear();
while (x > 1) {
int p = spf[x];
auto it = prime_id.find(p);
if (it == prime_id.end()) it = prime_id.emplace(p, pcnt++).first;
tmp.push_back(it->second);
while (x % p == 0) x /= p;
}
sort(tmp.begin(), tmp.end());
type_count[encode_vec(tmp)]++;
}
vector<pair<u128, int>> types;
types.reserve(type_count.size());
for (auto &p : type_count) types.push_back(p);
unordered_map<u128, int, U128Hash> sup_count;
sup_count.reserve(n * 16 + 10);
vector<vector<int>> type_vecs;
type_vecs.reserve(types.size());
for (auto &[key, cnt] : types) {
vector<int> v = decode_key(key);
type_vecs.push_back(v);
gen_subsets_from_vec(v, subs);
for (u128 s : subs) sup_count[s] += cnt;
}
unordered_set<u128, U128Hash> relevant;
relevant.reserve(types.size() * 64 + 10);
vector<int> uni;
vector<u128> usubs;
int mtypes = (int)type_vecs.size();
for (int i = 0; i < mtypes; i++) {
for (int j = i; j < mtypes; j++) {
uni = unite_vecs(type_vecs[i], type_vecs[j]);
gen_subsets_from_vec(uni, usubs);
for (u128 s : usubs) {
if (s) relevant.insert(s);
}
}
}
unordered_map<u128, int, U128Hash> avoid_memo;
avoid_memo.reserve(relevant.size() * 2 + 10);
auto avoid_count = [&](u128 key) -> int {
auto found = avoid_memo.find(key);
if (found != avoid_memo.end()) return found->second;
vector<int> v = decode_key(key);
gen_subsets_from_vec(v, subs);
int hit = 0;
for (u128 s : subs) {
if (!s) continue;
auto it = sup_count.find(s);
if (it == sup_count.end()) continue;
if (key_size(s) & 1) hit += it->second;
else hit -= it->second;
}
int res = n - hit;
avoid_memo[key] = res;
return res;
};
ll ans = 0;
for (u128 tkey : relevant) {
vector<int> tv = decode_key(tkey);
int r = (int)tv.size();
if (coef[r] == 0) continue;
gen_subsets_from_vec(tv, usubs);
ll cover = 0;
for (u128 u : usubs) {
int c = avoid_count(u);
ll pairs = (ll)c * (c - 1) / 2 % MOD;
if (key_size(u) & 1) cover -= pairs;
else cover += pairs;
}
cover %= MOD;
if (cover < 0) cover += MOD;
ans = (ans + cover * coef[r]) % MOD;
}
cout << ans << '\n';
}
}