PROBLEM LINK:
Practice
Contest: Division 1
Contest: Division 2
Contest: Division 3
Contest: Division 4
Author: aryansanghi
Tester: sushil2006
Editorialist: iceknight1093
DIFFICULTY:
TBD
PREREQUISITES:
Linearity of expectation
PROBLEM:
You’re given a graph G that’s not a path graph.
Let L(G) denote the line graph of G, and L^x(G) = L(L^{x-1}(G)).
You’re given x and k.
Let e^{(i)} denote the k-dimensional vector (0, 0, \ldots, 1024, \ldots, 0, 0) where every coordinate is 0 except for the i-th, with the i-th coordinate being 1024 instead.
To each vertex of G, assign a random vector among e^{(1)}, e^{(2)}, \ldots, e^{(k)}, each with equal probability.
Let t_G(v) denote the vector assigned to vertex v.
For the line graph L(G), define \displaystyle t_{L(G)}(v) = \frac{t_G(v_1) + t_G(v_2)}{2} assuming the vertex v in L(G) corresponds to edge (v_1, v_2) in G.
You’re also given a k-dimensional vector s.
Find the expected number of vertices in L^x(G) whose value equals s.
It’s guaranteed that 1 \le x \le 3.
EXPLANATION:
By linearity of expectation, computing the expected number of vertices whose value equals s can be done by computing the probability that each single vertex has a value of s, and summing all of these up.
Further, because the value of a vertex in a line graph is obtained by averaging the values of the corresponding vertices in the base graph (which is a linear combination), the value of any vertex in L^x(G) will be some linear combination of base values; since we’re just repeatedly taking linear combinations.
So, to compute the probability that some vertex v in L^x(G) has a final value of s, it would be quite useful to know exactly which vertices of G contribute to the values of v, along with their coefficients in the linear combination.
To find this, we use the fact that x is quite small; so we can solve for x = 1, 2, 3 separately.
Case 1: x = 1
This case is fairly easy.
Each vertex of L(G) corresponds to some edge of G.
So, if we compute p - the probability that a single vertex v has a value of s - then the answer will simply be p \cdot m by symmetry.
Thus, we focus on computing p.
The value p is probability that \displaystyle s = \frac{e^{(i)} + e^{(j)}}{2} where i and j are chosen independently and randomly in [1, k].
To find this, note that because each e^{(i)} has only a single non-zero coordinate, \frac{e^{(i)} + e^{(j)}}{2} can have at most two non-zero coordinates.
Now you can casework on whether it has one or two non-zero coordinates:
- If there’s only one non-zero coordinate in s, it must equal 1024, and the probability of obtaining it is \frac{1}{k^2}.
- If there are two non-zero coordinates in s, they must must both equal 512, and the probability of obtaining this is \frac{2}{k^2}.
- Any other configuration of s has a probability of 0.
This completely solves the x = 1 case.
Case 2: x = 2
The second case is not really all that different from the first.
Each vertex in L^2(G) corresponds to an edge in L(G).
So, say we’re looking at the edge (y_1, y_2) in L(G).
y_1 itself must correspond to some edge (v_1, v_2) in G, while y_2 must also correspond to some edge (v_3, v_4) in G.
However, vertices y_1 and y_2 are adjacent in L(G).
By definition, this is only possible if their corresponding edges in G share an endpoint.
So, if we let v_2 = v_4 be the shared endpoint, we have (v_1, v_2) and (v_2, v_3) as our edges.
This is just a length-2 path in the base graph!
(The length of a path is the number of edges in it.)
Each length-2 path in G thus corresponds to one vertex in L^2(G).
If the path is v_1 - v_2 - v_3, it can be verified that by taking averages, the corresponding weighting is:
- \frac{1}{4} for v_1 and v_3
- \frac{1}{2} for v_2
That is, we have
We now need to find the probability of this happening, and multiply it by the number of length 2 paths in G to obtain the expected value.
To find the probability: note that if s has more than 3 non-zero coordinates, it cannot be obtained as the weighted sum of three of the e^{(i)}'s.
So, s must have \le 3 positive coordinates.
While it is now possible to do some casework (as we did in the x = 1 case), there are many more valid configurations now.
Instead, it’s much simpler to simply brute force all possible configurations.
That is: without loss of generality, assume that the non-zero positions in s are 1, 2, 3.
(It’s possible that s has fewer than 3 non-zero positions, but that’s fine and won’t affect anything.)
Observe then that any e^{(i)} for i \gt 3 is useless, so each of t_G(v_1), t_G(v_2), t_G(v_3) must be some e^{(i)} for 1 \le i \le 3.
So, there are only 3^3 = 27 possibilities that need to be checked at all.
For each possibility, after fixing the values we can simply compute the weighted sum and check if it equals s or not.
This, divided by k^3, will give us the required probability for a single edge.
Finally, we need to multiply this by the number of length-2 paths in G.
Counting this is simple: if d(v) denotes the degree of v, the required count is just
because for each “central” vertex v, any two of its neighbors can be chosen.
Case 3: x = 3
Finally, we have x = 3.
Vertices in L^3(G) correspond to edges in L^2(G).
We saw that each vertex in L^2(G) was a length-2 path in G, so what do edges between them mean?
In fact, it means the two paths must share a common edge!
If we analyze what such configurations can look like in the base graph, there are three possibilities:
- A length-3 path, where the middle edge is common.
- A triangle, where any one of the three edges is common (and there’s exactly one way to choose two length-3 paths for each one).
- A “claw”, i.e. one vertex connected to three others, a.k.a K_{1, 3}.
Again any one of the three edges can be the common one, with exactly one way each.
We deal with each of these cases separately.
Case 3.1: Length-3 paths
Here, if the path is v_1 - v_2 - v_3 - v_4, it can be worked out that the coefficients in the linear combination are [\frac{1}{8}, \frac{3}{8}, \frac{3}{8}, \frac{1}{8}].
To count the number of assignments that result in s, we again simply brute force: if s has \gt 4 non-zero coordinates then the answer is 0, otherwise there are only 4^4 different assignments that need to be checked.
To count the number of such paths: fix the middle edge (v_2, v_3), then there are d(v_2)-1 choices for v_1 and d(v_3)-1 choices for v_4.
(Subtract one from each because we cannot choose v_2 or v_3 from the other neighbor.)
These choices are pretty much independent of each other - the only constraint is that we must ensure v_1 \ne v_4, so subtract out the number of common neighbors of v_2 and v_3 from the product (d(v_2)-1) \cdot (d(v_3)-1).
This can be easily done in \mathcal{O}(N^3) time which is fast enough for the given constraints (it’s however possible, though not needed, to achieve faster complexity.)
Case 3.2: Triangles
If the triangle vertices are (v_1, v_2, v_3), it can be worked out that the coefficients are [\frac{3}{8}, \frac{3}{8}, \frac{2}{8}], where the \frac{3}{8}'s correspond to the weights of the endpoints of whichever is the single common edge.
Brute force check all 3^3 possible assignments to find the probability of getting s.
Counting triangles is also trivial to do in \mathcal{O}(N^3) by just brute-forcing all possibilities (again, it’s possible to do this much faster, but not necessary here.)
Finally, note that each triangle in G actually corresponds to three vertices in L^3(G), depending on which edge is common.
So, there’s an additional multiplier of 3 here.
Case 3.3: Claws, a.k.a K_{1, 3}
The coefficients here are [\frac{4}{8}, \frac{2}{8}, \frac{1}{8}, \frac{1}{8}] where the \frac{4}{8} is for the central vertex, and the \frac{2}{8} goes to whichever vertex is the other endpoint of the common edge.
Again, there are only 4^4 configurations that need to be checked, so that’s not a problem.
As for counting: for each choice of central vertex v, there are \binom{d(v)}{3} choices of which three neighbors to pick.
Just as in the triangle case, this has an extra multiplier of 3, because after picking three edges, any one of them can be the “common” edge.
Add up the answers of all three cases to obtain the final answer.
TIME COMPLEXITY:
\mathcal{O}(N^3+(x+1)^{(x+1)}) per testcase.
CODE:
Editorialist's code (C++)
// #include <bits/allocator.h>
// #pragma GCC optimize("Ofast,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 fac = precalc_factorial<Zp>(200);
int t; cin >> t;
while (t--) {
int n, m, x, k; cin >> n >> m >> x >> k;
vector s(k, 0);
for (int &y : s) cin >> y;
vector adj(n, vector<int>());
for (int i = 0; i < m; ++i) {
int u, v; cin >> u >> v;
adj[--u].push_back(--v);
adj[v].push_back(u);
}
if (x == 1) {
// for each edge (u, v), average(t[u], t[v])
// either single 1024 or two 512
Zp ans = 0;
if (ranges::count(s, 1024) == 1 and ranges::count(s, 0) == k-1) {
// 1/k^2 per edge
ans = Zp(1) / (1ll*k*k);
ans *= m;
}
else {
if (ranges::count(s, 512) == 2 and ranges::count(s, 0) == k-2) {
// 2/k^2 per edge
ans = Zp(2) / (1ll*k*k);
ans *= m;
}
}
cout << ans << '\n';
}
else if (x == 2) {
// for each path (u, v, w), (t[u] + 2t[v] + t[w])/4
Zp paths = 0;
for (int i = 0; i < n; ++i) {
int deg = size(adj[i]);
paths += deg * (deg-1) / 2;
}
vector<int> reduced;
for (int y : s) if (y > 0) {
reduced.push_back(y);
}
int sz = size(reduced);
if (sz > 3) {
cout << 0 << '\n';
continue;
}
Zp ans = 0;
for (int i = 0; i < sz; ++i) for (int j = 0; j < sz; ++j) for (int ii = 0; ii < sz; ++ii) {
vector<int> cur(sz);
cur[i] += 1024 / 4;
cur[j] += 1024 / 2;
cur[ii] += 1024 / 4;
if (cur == reduced) ans += paths;
}
cout << ans / (1ll * k * k * k) << '\n';
}
else {
// for each path, (1, 3, 3, 1) / 8
// for each star, (1, 4, 2, 1) / 8
// for two sides of a triangle, (3, 3, 2) / 8
Zp prob = Zp(k) ^ 4;
prob = 1/prob;
auto common = [&] (int i, int j) {
vector mark(n, 0);
for (int u : adj[i]) mark[u] = 1;
int res = 0;
for (int u : adj[j]) res += mark[u];
return res;
};
vector<int> reduced;
for (int y : s) if (y > 0) {
reduced.push_back(y);
}
int sz = size(reduced);
if (sz > 4) {
cout << 0 << '\n';
continue;
}
Zp ans = 0;
for (int i = 0; i < n; ++i) {
for (int j : adj[i]) if (j > i) {
// fix this edge
Zp paths = (size(adj[i])-1) * (size(adj[j])-1);
Zp triangles = common(i, j);
paths -= triangles;
// path
for (int a = 0; a < sz; ++a) for (int b = 0; b < sz; ++b) for (int c = 0; c < sz; ++c) for (int d = 0; d < sz; ++d) {
vector<int> cur(sz);
cur[a] += 1024 / 8;
cur[b] += 3 * 1024 / 8;
cur[c] += 3 * 1024 / 8;
cur[d] += 1024 / 8;
if (cur == reduced) ans += paths * prob;
}
// triangle
for (int a = 0; a < sz; ++a) for (int b = 0; b < sz; ++b) for (int c = 0; c < sz; ++c) {
vector<int> cur(sz);
cur[a] += 2 * 1024 / 8;
cur[b] += 3 * 1024 / 8;
cur[c] += 3 * 1024 / 8;
if (cur == reduced) ans += triangles * prob * k;
}
}
int deg = size(adj[i]);
Zp stars = deg * (deg-1) * (deg-2) / 6;
for (int a = 0; a < sz; ++a) for (int b = 0; b < sz; ++b) for (int c = 0; c < sz; ++c) for (int d = 0; d < sz; ++d) {
vector<int> cur(sz);
cur[a] += 1024 / 8;
cur[b] += 4 * 1024 / 8;
cur[c] += 2 * 1024 / 8;
cur[d] += 1024 / 8;
if (cur == reduced) ans += 3 * stars * prob;
}
}
cout << ans << '\n';
}
}
}