PROBLEM LINK:
Author: Jatin Yadav
Tester: Encho Mishinev
DIFFICULTY:
HARD.
PREREQUISITES:
Matrix Tree Theorem, Generating Functions, Linearity of expectation
PROBLEM:
Given a graph on n vertices. Each second, you choose a random unordered pair (a, b); If (a, b) is an edge in the graph, remove it, else add it. Given Q queries each containing a time t, find the expected number of spanning trees in the graph after time t.
QUICK EXPLANATION:
For each k, find the number of spanning trees(say T_k) of the complete graph that contain exactly k edges from the original graph G = (V, E) . Also, find for each k, the probability (say P_k(t)) that a spanning tree with k edges from the original graph is alive at time t. Using linearity of expectation, the answer is \displaystyle \sum_{k=0}^{n-1} P_k(t) T_k . We use a generalization of Kirchoff’s theorem to compute T_k and generating functions to compute P_k(t).
EXPLANATION:
Let H be the set of all trees with |V| vertices. Call a tree T alive at time t if all edges of T are present in the graph at time t. The number of spanning trees equals the number of alive trees in H. So, using linearity of expectation, we have to find:
where T_k is the number of trees sharing k edges with G, and P_k is the probability that a tree with k edges common with G is alive at time t.
So, now we have two subproblems : computing T_k and P_k(t).
Computing T_k :
The Matrix Tree Theorem states that the number of spanning trees in a graph G is equal to any cofactor of its laplacian matrix, defined as:
In general, we can put labels on each edge. For each edge (i, j), lets introduce an indeterminate x_{ij}. Consider the matrix L(\bold{x}), whose entries are given by:
Now, any cofactor of L(\bold{x}) is a polynomial in the variables x_{ij}, where the contribution of a spanning tree is the monomial equal to the product of the indeterminates corresponding to the edges of the tree. That is, each tree T contributes \displaystyle \prod_{(i, j) \in T} x_{ij}.
The proof of this generalization is very similar to the proof of the matrix tree theorem:
Proof outline
We’ll give a proof analogous to the proof for kirchoff’s theorem on wikipedia. Consider a modified incidence matrix E(\bold{x}), where, In each column corresponding to e = (u, v) with u < v , we have E(\bold{x})_{ue} = \sqrt{x_{(u, v)}}, and E(\bold{x})_{ve} = -\sqrt{x_{(u, v)}}. Clearly, L(\bold{x}) = E(\bold{x}) E(\bold{x})^T. Let F(\bold{x}) be the matrix formed on removing the first row of E(\bold{x}) .
Then, by cauchy binet formula, we have \det(L(\bold{x})_{11}) = \displaystyle \sum_{S} \det(F_S(\bold{x}))^2,
where S varies over all size n-1 subsets of edges, and F_S(\bold{x}) denotes the square matrix corresponding to the columns of F(\bold{x}) indexed by S.
We claim that \det(F_S(\bold{x}))^2 = 0 if there is a cycle formed by edges in S, and \displaystyle \prod _{(i, j) \in S} x_{ij} if S is a spanning tree.
If there is a cycle, we can find a linear combination of the columns corresponding to edges in the cycle with a sum 0.
Else, we prove using induction that the determinant is \pm \displaystyle \prod _{(i, j) \in S} \sqrt{x_{ij}}.
The base case is when we have two nodes and the determinant clearly is -\sqrt{x_{12}}
For the induction step, consider a tree T with \geq 3 vertices. Consider any leaf u \neq 1. There is only one edge (say e = (u, v)) incident to u, and in the determinant expansion, u must be mapped to e. So, we can erase the row u and the column e. The remaining matrix corresponds to a tree T' formed on removing the leaf u and edge e from T. Clearly, \det(F_T(\bold{x})) = \pm \sqrt{x_{u, v}} \times \pm \displaystyle \prod_{(i, j) \in T'} \sqrt{x_{ij}} = \pm \displaystyle \prod_{(i, j) \in T} \sqrt{x_{ij}}.
Notice that, if we replace x_{ij} by y for each edge (i, j) in G and by 1 for each edge not in G, any cofactor will be \displaystyle \sum_{k=0}^{n-1} T_k y^k
To find the cofactor, we evaluate it at n different values of y and then interpolate. This takes O(N^4).
Computing P_k(t):
Consider a tree T. Let X be the set of edges in T \cap E, Y be T \setminus E, and \bar{T} be the complement of the set T. Let |X| = k. We know |\bar{T}| = (n - 1)(n - 2) / 2 (say r).
For T to be alive, each edge in X must have been toggled an even number of times and each edge in Y must have been toggled an odd number of times.
Let 2x_1, 2x_2, \ldots 2x_k be the ordered number of times edges in X are toggled, 2y_1+1, 2y_2 + 1, \ldots 2y_{n-1-k} + 1 be the ordered number of times edges in Y are toggled, and z_1, z_2, \ldots z_r be the ordered number of times edges in \bar{T} are toggled. Then, we have:
And for each such triple of tuples, we have to add \frac{t!}{\prod (2x_i)! \prod (2y_i+1)! \prod z_i!}.
So, the total number of ways to toggle edges is given by t! times the coefficient of x^t in
Evaluate Q_k(x) as a polynomial in e^x. Say Q_k(x) = \sum q_{k,i} e^{ix}. Then, we have:
Notice that q_{k, i} is non zero only for i of the form 2j + r - n + 1, where j varies from 0 to n - 1. Let H_i = \sum_{k}T_k q_{k, i} . Clearly, there are only n non zero H_i. The answer E(t) at time t is given by:
The precomputation of q_{k,i} for a single k takes O(N^2), so for all k, this can be done in O(N^3). Therefore, H_i can also be computed in O(N^3).
Answering a query at time t takes O(N \log{t}).
The overall complexity is O(N^4 + Q N \log{t_{\text{max}}}).
SOLUTIONS:
Setter's Solution
#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define sd(x) scanf("%d", &(x))
const int N = 1e5 + 10, mod = 1e9 + 7;
int fact[N], invfact[N];
inline int add(int x, int y){ x += y; if(x >= mod) x -= mod; return x;}
inline int sub(int x, int y){ x -= y; if(x < 0) x += mod; return x;}
inline int mul(int x, int y){ return (((ll) x) * y) % mod;}
inline int powr(int a, ll b){
int x = 1 % mod;
while(b){
if(b & 1) x = mul(x, a);
a = mul(a, a);
b >>= 1;
}
return x;
}
inline int inv(int a){ return powr(a, mod - 2);}
void pre(){
fact[0] = invfact[0] = 1;
for(int i = 1;i < N; i++) fact[i] = mul(i, fact[i - 1]);
invfact[N - 1] = inv(fact[N - 1]);
for(int i = N - 2; i >= 1; i--) invfact[i] = mul(invfact[i + 1], i + 1);
assert(invfact[1] == 1);
}
inline int C(int n, int k){
if(n < k || k < 0) return 0;
return mul(fact[n], mul(invfact[k], invfact[n - k]));
}
/*
Find a degree n - 1 polynomial passing through n points (x[i], y[i]) in O(n^2)
*/
vector<int> interpolate(const vector<int> & x, const vector<int> & y){
int n = x.size();
vector<int> ret(n);
vector<int> numerator(n + 1);
vector<int> val(n + 1);
numerator[0] = 1;
for(int i = 0; i < n; i++){
for(int j = i + 1; j >= 1; j--)
numerator[j] = sub(numerator[j - 1], mul(x[i], numerator[j]));
numerator[0] = mul(numerator[0], mod - x[i]);
}
for(int i = 0; i < n; i++){
int den = 1;
for(int j = 0; j < n; j++) if(j != i) den = mul(den, sub(x[i], x[j]));
int coeff = mul(inv(den), y[i]);
copy(numerator.begin(), numerator.end(), val.begin());
int a = x[i];
for(int j = n - 1; j >= 0; j--){
int h = val[j + 1];
val[j] = add(val[j], mul(a, h));
ret[j] = add(ret[j], mul(h, coeff));
}
assert(val[0] == 0);
}
return ret;
}
int determinant(vector<vector<int>> & V){
int ret = 1;
int n = V.size();
vector<int> done(n, 0);
for(int j = 0; j < n; j++){
bool found = false;
for(int i = 0; i < n; i++) if(V[i][j] && !done[i]){
int coeff = inv(V[i][j]);
for(int k = 0; k < n; k++) if(!done[k] && i != k && V[k][j]){
// Rk -> Rk - V[k][j] / V[i][j] * Ri
int factor = mul(V[k][j], coeff);
for(int l = 0; l < n; l++){
V[k][l] = sub(V[k][l], mul(factor, V[i][l]));
}
}
found = true;
done[i] = 1;
ret = mul(ret, V[i][j]);
break;
}
if(!found) return 0;
}
return ret;
}
struct graph{
vector<vector<int>> col;
vector<int> deg;
int n;
graph(int n) : n(n), col(n, vector<int>(n, 0)), deg(n){}
void add_edge(int u, int v){
col[u][v] = col[v][u] = 1;
deg[u]++; deg[v]++;
}
vector<int> getSpanningTrees(){
// how many spanning trees of the complete graph with exactly i edges from this graph and the rest from its complement?
vector<int> x(n), y(n);
for(int i = 0; i < n; i++){
vector<vector<int>> arr(n - 1, vector<int>(n - 1, 0));
for(int j = 1; j < n; j++) for(int k = 1; k < n; k++){
if(j == k) arr[j - 1][k - 1] = deg[j] * i + n - 1 - deg[j];
else arr[j- 1][k - 1] = sub(0, col[j][k] * i + 1 - col[j][k]);
}
x[i] = i;
y[i] = determinant(arr);
}
return interpolate(x, y);
}
};
vector<int> naive_mul(const vector<int> & X, const vector<int> & Y){
vector<int> Z(X.size() + Y.size() - 1);
for(int i = 0; i < X.size(); i++)
for(int j = 0; j < Y.size(); j++)
Z[i + j] = add(Z[i + j], mul(X[i], Y[j]));
return Z;
}
const int L = 32005;
int powers[4][305][L + 1];
int main(){
pre();
int n = 100, m, q;
sd(n); sd(m); sd(q);
graph g(n);
for(int i = 0; i < m; i++){
int u, v; sd(u); sd(v);
u--; v--;
g.add_edge(u, v);
}
vector<int> T = g.getSpanningTrees();
int other = ((n - 2) * (n - 1)) / 2;
int coeff = powr((mod + 1) >> 1, n - 1);
vector<int> H(n);
for(int k = 0; k < n; k++){
/*
Given a spanning tree T with k original edges, what is the probability that it is alive?
Original edges in T have to be flipped an even number of times, 2*x_1, 2*x_2, ..., 2*x_k
Non-original edges in T have to be flipped an odd number of times, 2*y_1+1, 2*y_2+1, .., 2*y_{n-1-k}+1
Edges not in T can be flipped any number of times -> z_1, z_2, ..., z_other
Total number of ways -> t! / prod(2*x_i!) / prod((2 * y_i + 1)!) / prod(z_i!)
= t! * [x^t] ((e^x + e^-x) / 2)^k * ((e^x - e^-x) / 2) ^ (n - 1 - k) * e^(other * x)
Let y = e^x
(y^2 + 1) ^ k * (y^2 - 1) ^ (n - 1 - k) * y^(other - n + 1)
*/
vector<int> A(k + 1), B(n - k);
for(int i = 0; i <= k; i++) A[i] = C(k, i);
for(int i = 0; i <= n - 1 - k; i++){
B[i] = C(n - 1 - k, i);
if((n - 1 - k - i) & 1) B[i] = sub(0, B[i]);
}
vector<int> C = naive_mul(A, B);
for(int i = 0; i < n; i++) H[i] = add(H[i], mul(T[k], C[i]));
}
for(int i = 0; i < n; i++){
int v = 2 * i + other - n + 1;
if(v < 0) v = sub(mod, -v);
powers[0][i][0] = powers[1][i][0] = powers[2][i][0] = powers[3][i][0] = 1;
for(int j = 1; j <= L; j++){
powers[0][i][j] = mul(v, powers[0][i][j - 1]);
}
v = powers[0][i][L];
for(int j = 1; j <= L; j++){
powers[1][i][j] = mul(v, powers[1][i][j - 1]);
}
v = powers[1][i][L];
for(int j = 1; j <= L; j++){
powers[2][i][j] = mul(v, powers[2][i][j - 1]);
}
v = powers[2][i][L];
for(int j = 1; j <= L; j++){
powers[3][i][j] = mul(v, powers[3][i][j - 1]);
}
}
function<int(int, ll)> getPower = [&](int i, ll t){
int ret = 1;
for(int j = 0; j < 4; j++){
ll dv = t / L;
int md = t - L * dv;
ret = mul(ret, powers[j][i][md]);
t = dv;
}
assert(t == 0);
return ret;
};
// O(n) per query
function<int(ll)> query = [&](ll t){
int ret = 0;
for(int i = 0; i < n; i++){
// t! * [t]exp((2 * i + other - n + 1)x)
ret = add(ret, mul(H[i], getPower(i, t)));
}
int total = (n * (n - 1)) / 2;
return mul(coeff, mul(inv(powr(total, t)), ret));
};
while(q--){
ll t;
scanf("%lld", &t);
printf("%d\n", query(t));
}
}
