PROBLEM LINK:
Practice
Contest: Division 1
Contest: Division 2
Contest: Division 3
Contest: Division 4
Author: iceknight1093
Tester: sushil2006
Editorialist: iceknight1093
DIFFICULTY:
Medium
PREREQUISITES:
Dynamic programming, elementary combinatorics
PROBLEM:
Given an array, count the number of pairs of disjoint subsequences such that the MEX of one subsequence equals the \gcd of the other.
EXPLANATION:
Since we’re counting pairs of subsequences with \gcd = \text{mex}, a natural first step is to fix the value of the \gcd and count the number of pairs.
Let’s fix the value to be K.
Now, for each element of the array, we have to decide whether it goes into the first subsequence, the second subsequence, or neither.
First, let’s look at the mex.
For the mex of a subsequence to be K, the only condition is that the subsequence should contain all the elements 1, 2, 3, \ldots, K-1 at least once, and shouldn’t contain K.
Let f_x denote the number of occurrences of x in the array.
Simple combinatorics shows that there are
ways to achieve this; simply by taking the number of ways of choosing a non-empty subset of each element and multiplying them all.
What about the remaining copies of elements \lt K?
Well, they certainly can’t be in the gcd subsequence, since that would cause the gcd itself to be \lt K.
So, all the remaining elements \lt K must be in neither subsequence, i.e. there’s only one option for each of them.
We’re now left with only elements \geq K.
We want to select some of them such that their gcd equals K, and then the remaining ones can be either put into the mex subsequence, or left in neither one.
Since the gcd must equal K, all the chosen elements must themselves be multiples of K.
In particular, this means all elements \geq K that are not multiples of K have two choices each: either go into the mex subsequence, or be in neither one.
This is just another multiplier of a power of two and easy enough to compute.
This leaves us with only multiples of K.
Let M_K denote the number of multiples of K.
Suppose we choose some subset S of these multiples, such that \gcd(S) = K.
There are then M_K - |S| elements remaining, and each of them have two choices: go to the mex subsequence, or not.
So, that’s 2^{M_K - |S|} ways.
Our goal is to sum this up across all the different ways of choosing a valid subset S.
(You might note that this isn’t quite accurate – specifically, unused occurrences of K cannot go into the mex subsequence and so have only one option rather than two. Ignore this discrepancy for now, we’ll come back to it.)
This computation can be done using a fairly standard Möbius-style dynamic programming.
First, note that we can take out the common multiplier of 2^{M_K}, so that what we’re really trying to compute is the sum of 2^{-|S|} across all valid subsets.
This is helpful because the value contributed by a subset is now intrinsic to it, rather than being dependent on K.
Let’s define dp_K to be the sum of 2^{-|S|} across all subsets whose gcd equals K.
To compute this:
- There are M_K multiples of K, so only subsets of these M_K elements need to be considered.
- Each of these subsets will definitely have its gcd be a multiple of K.
- So, if we initialize dp_K with the sum across all subsets of the M_K elements, we can then subtract out the contributions corresponding to those gcd’s that are multiples of K.
That is, overestimate dp_K initially, and then fix it by subtracting out the values
dp_{2K}, dp_{3K}, dp_{4K}, \ldots
The only question remaining now is how to initialize dp_K.
Note that the initial value is the sum of 2^{-|S|} across all non-empty subsets of M_K elements.
Since sets of the same size have the same contribution, one way to write this summation is
This form is pretty much exactly what we obtain from the binomial expansion though, in particular it can be seen that
We want exactly the RHS minus the term for i = 0, so a closed form for what we’re looking for is simply
We now know how to compute all the dp_K values: initialize dp_K with \left(1 + \frac{1}{2}\right)^{M_K} - 1, and then subtract out the values of dp_{i\cdot K} for all i\gt 1.
Once this is done, for a fixed K we simply need to take the product of:
- \prod_{x=1}^{K-1} \left(2^{f_x} - 1\right)
- 2^y, where y is the number of elements \geq K that are not multiples of K.
- 2^{M_K} \cdot dp_K
The final answer is the sum of this across all 1 \leq K \leq \max(A).
However, there’s a couple more details.
The biggest issue is one already pointed out above: dp_K was computed by just assuming that every multiple of K not in the subset had two choices, but this isn’t always true – specifically for occurrences of K.
This, luckily, is not too hard to fix.
First, note that if there are no occurrences of K, then the dp is already correct and nothing needs to be done.
So, what we can do is split the counting into two parts: subsets that don’t include any occurrences of K, and those that do.
For subsets that don’t include any occurrences of K, observe that none of the occurrences of K matter at all anymore, so we functionally have only M_K - f_K elements to work with.
Simply running the dp with this count, rather than M_K, will give us the value we want.
As for subsets that do include at least one occurrence of K: note that as soon as K is included, the gcd of the subset is guaranteed to be K.
So, we can:
- Choose any non-empty subsequence of the K's.
There are 2^{f_K}-1 ways of doing this. - For every multiple of K that’s not K itself, it now has three options: go into one of the subsequences, or go into neither.
So, that’s 3^{M_K - f_K} possibilities.
This takes care of the issue.
Finally, there’s one last thing we haven’t considered: are the subsequences we choose non-empty?
For that, it can be seen that:
- The subsequence corresponding to gcd was never allowed to be empty.
- The subsequence corresponding to mex must include values 1, 2, 3\ldots, K-1
So, it can only potentially be empty if K = 1.
In particular, note that the mex subset will have been empty exactly once for each subsequence with gcd equal to 1.
So, we subtract the count of subsequences with gcd equal to 1 from the above computed answer to obtain the true answer.
The complexity of this is \mathcal{O}(N + \max(A)\log\max(A)) which is fast enough for our purposes.
TIME COMPLEXITY:
\mathcal{O}(N + \max(A)\log\max(A)) per testcase.
CODE:
Editorialist's code (C++)
// #include <bits/allocator.h>
// #pragma GCC optimize("O3,unroll-loops")
// #pragma GCC target("avx2,bmi,bmi2,lzcnt,popcnt")
#include "bits/stdc++.h"
using namespace std;
using ll = long long int;
mt19937_64 RNG(chrono::high_resolution_clock::now().time_since_epoch().count());
/**
* Integers modulo p, where p is a prime
* Source: Aeren (modified from tourist?)
* Modmul for 64-bit mod from kactl:ModMulLL
* Works with p < 7.2e18 with x87 80-bit long double, and p < 2^52 ~ 4.5e12 with 64-bit
*/
template<typename T>
struct Z_p{
using Type = typename decay<decltype(T::value)>::type;
static vector<Type> MOD_INV;
constexpr Z_p(): value(){ }
template<typename U> Z_p(const U &x){ value = normalize(x); }
template<typename U> static Type normalize(const U &x){
Type v;
if(-mod() <= x && x < mod()) v = static_cast<Type>(x);
else v = static_cast<Type>(x % mod());
if(v < 0) v += mod();
return v;
}
const Type& operator()() const{ return value; }
template<typename U> explicit operator U() const{ return static_cast<U>(value); }
constexpr static Type mod(){ return T::value; }
Z_p &operator+=(const Z_p &otr){ if((value += otr.value) >= mod()) value -= mod(); return *this; }
Z_p &operator-=(const Z_p &otr){ if((value -= otr.value) < 0) value += mod(); return *this; }
template<typename U> Z_p &operator+=(const U &otr){ return *this += Z_p(otr); }
template<typename U> Z_p &operator-=(const U &otr){ return *this -= Z_p(otr); }
Z_p &operator++(){ return *this += 1; }
Z_p &operator--(){ return *this -= 1; }
Z_p operator++(int){ Z_p result(*this); *this += 1; return result; }
Z_p operator--(int){ Z_p result(*this); *this -= 1; return result; }
Z_p operator-() const{ return Z_p(-value); }
template<typename U = T>
typename enable_if<is_same<typename Z_p<U>::Type, int>::value, Z_p>::type &operator*=(const Z_p& rhs){
#ifdef _WIN32
uint64_t x = static_cast<int64_t>(value) * static_cast<int64_t>(rhs.value);
uint32_t xh = static_cast<uint32_t>(x >> 32), xl = static_cast<uint32_t>(x), d, m;
asm(
"divl %4; \n\t"
: "=a" (d), "=d" (m)
: "d" (xh), "a" (xl), "r" (mod())
);
value = m;
#else
value = normalize(static_cast<int64_t>(value) * static_cast<int64_t>(rhs.value));
#endif
return *this;
}
template<typename U = T>
typename enable_if<is_same<typename Z_p<U>::Type, int64_t>::value, Z_p>::type &operator*=(const Z_p &rhs){
uint64_t ret = static_cast<uint64_t>(value) * static_cast<uint64_t>(rhs.value) - static_cast<uint64_t>(mod()) * static_cast<uint64_t>(1.L / static_cast<uint64_t>(mod()) * static_cast<uint64_t>(value) * static_cast<uint64_t>(rhs.value));
value = normalize(static_cast<int64_t>(ret + static_cast<uint64_t>(mod()) * (ret < 0) - static_cast<uint64_t>(mod()) * (ret >= static_cast<uint64_t>(mod()))));
return *this;
}
template<typename U = T>
typename enable_if<!is_integral<typename Z_p<U>::Type>::value, Z_p>::type &operator*=(const Z_p &rhs){
value = normalize(value * rhs.value);
return *this;
}
template<typename U>
Z_p &operator^=(U e){
if(e < 0) *this = 1 / *this, e = -e;
Z_p res = 1;
for(; e; *this *= *this, e >>= 1) if(e & 1) res *= *this;
return *this = res;
}
template<typename U>
Z_p operator^(U e) const{
return Z_p(*this) ^= e;
}
Z_p &operator/=(const Z_p &otr){
Type a = otr.value, m = mod(), u = 0, v = 1;
if(a < (int)MOD_INV.size()) return *this *= MOD_INV[a];
while(a){
Type t = m / a;
m -= t * a; swap(a, m);
u -= t * v; swap(u, v);
}
assert(m == 1);
return *this *= u;
}
template<typename U> friend const Z_p<U> &abs(const Z_p<U> &v){ return v; }
Type value;
};
template<typename T> bool operator==(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value == rhs.value; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator==(const Z_p<T>& lhs, U rhs){ return lhs == Z_p<T>(rhs); }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator==(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) == rhs; }
template<typename T> bool operator!=(const Z_p<T> &lhs, const Z_p<T> &rhs){ return !(lhs == rhs); }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator!=(const Z_p<T> &lhs, U rhs){ return !(lhs == rhs); }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> bool operator!=(U lhs, const Z_p<T> &rhs){ return !(lhs == rhs); }
template<typename T> bool operator<(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value < rhs.value; }
template<typename T> bool operator>(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value > rhs.value; }
template<typename T> bool operator<=(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value <= rhs.value; }
template<typename T> bool operator>=(const Z_p<T> &lhs, const Z_p<T> &rhs){ return lhs.value >= rhs.value; }
template<typename T> Z_p<T> operator+(const Z_p<T> &lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) += rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator+(const Z_p<T> &lhs, U rhs){ return Z_p<T>(lhs) += rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator+(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) += rhs; }
template<typename T> Z_p<T> operator-(const Z_p<T> &lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) -= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator-(const Z_p<T>& lhs, U rhs){ return Z_p<T>(lhs) -= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator-(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) -= rhs; }
template<typename T> Z_p<T> operator*(const Z_p<T> &lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) *= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator*(const Z_p<T>& lhs, U rhs){ return Z_p<T>(lhs) *= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator*(U lhs, const Z_p<T> &rhs){ return Z_p<T>(lhs) *= rhs; }
template<typename T> Z_p<T> operator/(const Z_p<T> &lhs, const Z_p<T> &rhs) { return Z_p<T>(lhs) /= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator/(const Z_p<T>& lhs, U rhs) { return Z_p<T>(lhs) /= rhs; }
template<typename T, typename U, typename enable_if<is_integral<U>::value>::type* = nullptr> Z_p<T> operator/(U lhs, const Z_p<T> &rhs) { return Z_p<T>(lhs) /= rhs; }
template<typename T> istream &operator>>(istream &in, Z_p<T> &number){
typename common_type<typename Z_p<T>::Type, int64_t>::type x;
in >> x;
number.value = Z_p<T>::normalize(x);
return in;
}
template<typename T> ostream &operator<<(ostream &out, const Z_p<T> &number){ return out << number(); }
/*
using ModType = int;
struct VarMod{ static ModType value; };
ModType VarMod::value;
ModType &mod = VarMod::value;
using Zp = Z_p<VarMod>;
*/
// constexpr int mod = 1e9 + 7; // 1000000007
constexpr int mod = (119 << 23) + 1; // 998244353
// constexpr int mod = 1e9 + 9; // 1000000009
using Zp = Z_p<integral_constant<decay<decltype(mod)>::type, mod>>;
template<typename T> vector<typename Z_p<T>::Type> Z_p<T>::MOD_INV;
template<typename T = integral_constant<decay<decltype(mod)>::type, mod>>
void precalc_inverse(int SZ){
auto &inv = Z_p<T>::MOD_INV;
if(inv.empty()) inv.assign(2, 1);
for(; inv.size() <= SZ; ) inv.push_back((mod - 1LL * mod / (int)inv.size() * inv[mod % (int)inv.size()]) % mod);
}
template<typename T>
vector<T> precalc_power(T base, int SZ){
vector<T> res(SZ + 1, 1);
for(auto i = 1; i <= SZ; ++ i) res[i] = res[i - 1] * base;
return res;
}
template<typename T>
vector<T> precalc_factorial(int SZ){
vector<T> res(SZ + 1, 1); res[0] = 1;
for(auto i = 1; i <= SZ; ++ i) res[i] = res[i - 1] * i;
return res;
}
int main()
{
ios::sync_with_stdio(false); cin.tie(0);
auto p2 = precalc_power<Zp>(2, 200005);
auto p3 = precalc_power<Zp>(3, 200005);
int t; cin >> t;
while (t--) {
int n; cin >> n;
vector a(n, 0);
for (int &x : a) cin >> x;
int m = ranges::max(a);
vector freq(m+1, 0), mulct(m+1, 0);
for (int x : a) ++freq[x];
for (int i = 1; i <= m; ++i) for (int j = i; j <= m; j += i)
mulct[i] += freq[j];
Zp ans = 0, pref = 1;
int rem = n;
vector dp(m+1, Zp(0));
vector ways(m+1, Zp(0));
for (int i = m; i >= 1; --i) {
dp[i] = p3[mulct[i]] / p2[mulct[i]];
dp[i] -= 1;
ways[i] = p2[mulct[i]] - 1;
for (int j = 2*i; j <= m; j += i) {
dp[i] -= dp[j];
ways[i] -= ways[j];
}
}
for (int i = 1; i <= m; ++i) {
if (pref == 0) break;
// no occurrence of i
Zp cur = p3[mulct[i] - freq[i]] / p2[mulct[i] - freq[i]] - 1;
for (int j = 2*i; j <= m; j += i) cur -= dp[j];
ans += cur * pref * p2[rem - mulct[i]] * p2[mulct[i] - freq[i]];
// yes occurrence of i
ans += pref * p2[rem - mulct[i]] * (p2[freq[i]] - 1) * p3[mulct[i] - freq[i]];
pref *= p2[freq[i]] - 1;
rem -= freq[i];
}
ans -= ways[1];
cout << ans << '\n';
}
}