PROBLEM LINK:
Practice
Div-2 Contest
Div-1 Contest
Author: Sayantan Jana
Tester: Radoslav Dimitrov
Editorialist: Sayantan Jana
DIFFICULTY:
HARD
PREREQUISITES:
Disjoint Set Union, Dynamic Connectivity
PROBLEM:
Initially given is a prime p and a matrix M of size N * N .
The matrix is required to be curious, i.e., once it is completely filled,
- Each cell should contain an integer between 1 and p-1 inclusive .
- For each non trivial square submatrix A of M (a submatrix containing more than 1 cell) , its determinant |A| should be a multiple of p.
You are required to handle Q queries . After each addition or delete query in the matrix M, report the number of possible ways to fill up the empty cells in M to make it curious modulo (1e9+7) .
QUICK EXPLANATION:
- It can be shown that the 2nd curious condition is equivalent to the matrix having rank 1 with respect to modulo p. In a rank 1 matrix, each row is a multiple of every other row and each column is a multiple of every other column . Hence there is a cyclic dependency between the cells, i.e., 1 \leq a,b,c,d \leq N, M_{a,b}*M_{c,d} \equiv M_{a,d}*M_{c,b} mod p, thus if any 3 of these cells have values in them, the 4th cell necesaarily needs to have a particular value to obey the rank 1 matrix condition .
- There are actually a total of 2N-1 individual cells whose values independently can determine the values to be filled in all other cells for the matrix to have rank 1 . Intuitively just consider the 1st row and 1st column filled, a total 2N-1 cells . These are actually the independent constants in a graph, the count of free constants available decide the count of ways to form a curious matrix, (p-1)^f , f being the current count of free constants.
- The problem can be modelled into a graph problem. Consider 2N nodes in the graph, each node representative of a row or a column, nodes 1 to N for the rows and N+1 to 2N for the columns . For every filled cell M_{x, y} containing value v, consider a directed edge from node x to node y+N with weight v and another directed edge from node y+N to node x with weight 1/v (inverse of v w.r.t prime p) . Now to obey the rank 1 matrix condition, the products of all edge weights in a cycle in the graphymust be 1 . The count of free constants is decided by the number of disconnected components in graph subtracted by 1.
- The addition query alone could have been handled in O(log N) by a DSU. However to handle deletion queries, dynamic connectivity comes into play .
EXPLANATION:
Counting ways to get a curious matrix from an initially empty matrix
Let’s count the number of ways to fill an empty matrix to make it curious .
If we were asked to just satisfy the first condition only, the number of ways would have been (p-1)^{N*N} . Does the count remain same if we were asked to incorporate the 2nd condition as well ? No, it doesn’t since by the second condition, there is certain dependency imposed on certain cells once some other cells are filled . This dependency is what is enforced by the rank 1 property of the matrix .
Proving the 2nd Curious condition is equivalent to rank 1 requirement
Consider a,b such that 1 \leq a,b \leq N-1, consider the 2*2 submatrix formed by rows a and a+1 and columns b and b+1 as,
The curious condition requires the determinant of the submatrix, |A| mod p = 0, i.e., M_{a,b}*M_{a+1,b+1} \equiv M_{a+1,b}*M_{a,b+1} mod p . This implies, \frac{M_{a,b}}{M_{a+1,b}} \equiv \frac{M_{a,b+1}}{M_{a+1,b+1}} mod p . Thus the bottom row of submatrix A is a multiple of the top row, now varying b from 1 to N-1, it can be found that (a+1)-th row is a multiple of a-th row in the matrix M. It can also be said the other way round, i.e., a-th row is a multiple of (a+1)-th row This can be generalised to show that each row is a multiple of every other row in a curious matrix . Similarly, this can also be shown for columns too . Hence, considering modulo prime p, the rank of any curious matrix is 1 .
Now we understand the cyclic dependency that is created among the cells . Taking cue from the rank 1 property we just proved, we use the fact there are N + (M - 1) free constants in a Rank 1 matrix of dimension N*M, the N-components of the first column vector (let’s call them x_0, x_1 … x_{N-1}), and the multiplier of each of the other column vectors with respect to the first (let’s call them y_1, y_2 … y_{m-1}). Hence the number of ways to fill up an empty matrix to make it curious is (p-1)^{2N-1}, since the answer is the number of values each constant can take to the power of number of such constants.
To better understand the notion of rank 1 matrices and free constants, look watch this video.
Transformation to a Graph problem
Since, we have already established that the rows are multiples of each other and columns are multiples of each other too, we have to somehow effecitively store the ratios between them to handle the queries .
Consider the matirx in this state ,
At this stage, we know c-th row is a multiple of a-th row by a factor of M_{c,b}/M_{a,b}, also d-th row is a multiple of b-th row by a factor of M_{a,d}/M_{a,b} . We attempt to maintain this ratios between these rows and columns together in a graph such that :
- The graph contains 2N nodes, each node representative of a row or a column, hence 2N of them, x-th row denoted by node x, y-th row denoted by node y+N .
- For every filled cell M_{x, y} containing value v, consider a directed edge from node x to node y+N with weight v and another directed edge from node y+N to node x with weight 1/v (inverse of v w.r.t prime p) .
On careful observation on the construction of the graph, we can infer that the ratio between a-th row and c-th can be obtained as the product of the edge weights in the path from node c to a . Notice that there can be many directed paths from c to a, but the product of each of those paths should be same . By construction the reverse path would have a product that is an inverse of the product of edge weights in original path . Hence, we can say, it is needed to satisfy that product of edge weights in each directed cycle in the graph should be 1. Another crucial observation is that the count of free constants in the matrix being constructed is equal to the number of components in the graph subtracted by 1.
Let’s now look at how to handle the queries.
Now on an addition query x y v,
- If there was earlier a path between x and y+N, it needs to obey that the product of all the edges in the path from node x to node y+N should be v .
- If there wasn’t a path earlier, we add an edge as described above. In notion of the rank 1 matrix count of the free constants decrease by 1 since number of components decrease and hence after this query the number of ways to obtain curious matrix decreases .
Answering the query
If the query is valid, i.e., there exists a way to fill up the remaining cells of the matrices , the answer is calculated as, (p-1)^f , f being the current count of free constants. Count of free constants can be found by getting the current number of components in the graph. Report the answer modulo 1e9+7.
Efficiently handling addition queries through Disjoint Set Union
Both the above checks can be handled faster by maintaing a disjoint set union structure . Note usually the DSU’s are associated with undirected edges, we too can work with undirected edges only since between any two nodes, the weight in one of the directed edge is just inverse of the the directed edge in opposite direction , this can be maintained using undirected edges only, just maintain the edge that is directed from a node to its parent in the DSU structure .
For an addition query x y v,
- To check if there is a path between x and y+N, we just need to check if the two nodes are in same subset of DSU or not .
- If the nodes are in the same component, find p_x, the product of all edges from x leading up to r_x, the root of the subset containg the node x and p_y, the product of all edges from y+N leading up to r_y,the root of the subset containg the node y+N (All products are found modulo prime p) . Note that we have actually simulated the reverse path from y+N to r_y . Since for obeying rank 1 matrix rule, the product of the path from x to y+N was required to be v and hence it must hold that,
v \equiv \frac{p_x}{p_y} mod p . - If the nodes are not in the same component, we need to connect r_x and r_y, the roots of the subsets containing x and y+N . We use smaller to bigger merging and accordingly take care of the edge weight we put between them . If r_y is made a child of r_x, the edge connecting r_y to r_x (we only maintain the edge in upward direction) is assigned an edge weight, w = \frac{p_x}{v*p_y} mod p , while if r_x is made a child of r_y, it is assigned an edge weight, w = \frac{v*p_y}{p_x} mod p . The number of components get decreased and in sense the free constants as well.
But note that deletion queries cannot be handled efficiently just by DSU . Hence so far we have solved subtask 1 and 2 . For subtask 1, it is not needed to handle deletions efficiently .
Efficiently Handling Deletion Queries through Dynamic Connectivity
Finally, Handling delete queries is a simple albeit programmatically tedious modification, what follows is a standard technique described here: Deletion in O(log N). We need to add delete operations to the DSU, for which we use a Segment Tree and a regular DSU as described above which does and undoes operations on a stack (by storing all changes made after each operation and undoing them later). We loop over the query list and pair up the set and unset queries, and the ranges in which the value is set, we add it to the Segment Tree. Then we perform a DFS on the Segment Tree, passing the DSU around, and this allows all the operation to be applied as deleted in the stack-order. Since path compression is not possible with deletion (roll backs), merging smaller into bigger is why we used in DSU as well .
For subtask 3, it is actually not needed to maintain edge weights, just maintaing the edges is fine . Complexity wise it’s same as in subtask 4, just the careful handling of edge weights is what is not needed in the subtask .
For subtask 4, we use everything discussed so far to efficiently answer the queries in O(Q log ^{2}Q) .
If you are new to dynamic connectivity, it’s recommended to try the problem Extending Set of Points and look at the attached editorial.
Also give a try at Dynamic connectivity contest.
SOLUTIONS:
Setter's Solution
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef vector<long long> vll;
typedef pair<long long, long long> pll;
long long MOD1;
const long long MOD2 = 1000000007;
ll mod_power(ll a, ll b, ll MOD) {
ll cumulative = a, result = 1;
for (; b > 0; b /= 2) {
if (b % 2 == 1)
result = (result * cumulative) % MOD;
cumulative = (cumulative * cumulative) % MOD;
}
return result;
}
class DynamicConnectivity {
void __dfs(int v, int l, int r, vll& res) {
ll last_ans = answer;
int state = save_ptr;
for (auto x : tree[v])
merge(x.u, x.v, x.ratio);
if (l == r - 1)
res[l] = answer;
else {
int m = (l + r) / 2;
__dfs(v * 2 + 1, l, m, res);
__dfs(v * 2 + 2, m, r, res);
}
while (save_ptr != state)
rollback();
answer = last_ans;
};
public:
int size_nodes;
int size_query;
vector<ll> parent, comp_size;
ll answer = 0;
vector<ll*> saved_object;
vector<ll> saved_value;
int save_ptr = 0;
vector<ll> factor;
ll comp_count;
struct Query {
int u, v;
ll ratio;
Query(pair<int, int> p, ll r) {
u = p.first, v = p.second;
ratio = r;
}
};
vector<vector<Query>> tree;
DynamicConnectivity(int n = 600000, int q = 300000) {
size_nodes = n;
size_query = q;
parent = vector<ll>(n);
comp_size = vector<ll>(n, 1);
ll tree_size = 1;
while (tree_size < q)
tree_size <<= 1;
tree = vector<vector<Query>>(2 * tree_size);
iota(parent.begin(), parent.end(), 0);
saved_object = vector<ll*>(max(3 * n, 1000000));
saved_value = vector<ll>(max(3 * n, 1000000));
factor = vector<ll>(n, 1);
comp_count = n;
}
void change(ll& object, ll value) {
saved_object[save_ptr] = &object;
saved_value[save_ptr] = object;
object = value;
save_ptr++;
}
void rollback() {
save_ptr--;
(*saved_object[save_ptr]) = saved_value[save_ptr];
}
int find(int x) {
if (parent[x] == x)
return x;
return find(parent[x]);
}
ll find_factor(int x) {
if (parent[x] == x)
return 1;
return 1ll*factor[x]*find_factor(parent[x])%MOD1;
}
void merge(int x, int y, ll ratio) {
ll factor_x = find_factor(x);
ll factor_y = find_factor(y);
x = find(x);
y = find(y);
if (x == y) {
if (!(factor_x == (ratio * factor_y) % MOD1))
change(comp_count, 0);
answer = comp_count;
return;
}
ll tmp_var = 1ll*ratio*factor_y%MOD1;
if (comp_size[x] > comp_size[y]) {
change(parent[y], x);
change(comp_size[x], comp_size[x] + comp_size[y]);
change(factor[y], (factor_x * mod_power(tmp_var, MOD1 - 2, MOD1)) % MOD1);
change(comp_count, comp_count - 1);
}
else
{
change(parent[x], y);
change(comp_size[y], comp_size[x] + comp_size[y]);
change(factor[x], (tmp_var * mod_power(factor_x, MOD1 - 2, MOD1)) % MOD1);
change(comp_count, comp_count - 1);
}
answer = comp_count;
}
void add(int l, int r, Query edge, int node = 0, int x = 0, int y = -1) {
if (y == -1)
y = size_query;
if (l >= r)
return;
if (l == x && r == y)
tree[node].emplace_back(edge);
else {
int m = (x + y) / 2;
add(l, min(r, m), edge, node * 2 + 1, x, m);
add(max(m, l), r, edge, node * 2 + 2, m, y);
}
}
vll results(int v = 0, int l = 0, int r = -1) {
if (r == -1)
r = size_query;
vll vec(size_query);
__dfs(v, l, r, vec);
return vec;
}
};
int main() {
ios::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
ll n, q;
cin >> n >> q >> MOD1;
map<pll, ll> last;
DynamicConnectivity dsu(n + n, q);
vector<DynamicConnectivity::Query> queries;
queries.reserve(q);
for (int i = 0; i < q; i++) {
ll x, y, val;
cin >> x >> y >> val;
x--, y--;
pll p(x, y + n);
queries.emplace_back(p, val);
if (last.count(p)) {
dsu.add(last[p], i, queries.back());
last.erase(p);
} else {
last[p] = i;
}
}
for (auto x : last)
dsu.add(x.second, q, queries[x.second]);
vll res = dsu.results();
for (int i = 0; i < q; i++)
cout << ((res[i] <= 0) ? 0 : mod_power(MOD1-1, res[i] - 1, MOD2)) << "\n";
}
Tester's Solution
#include <bits/stdc++.h>
#define endl '\n'
#define SZ(x) ((int)x.size())
#define ALL(V) V.begin(), V.end()
#define L_B lower_bound
#define U_B upper_bound
#define pb push_back
using namespace std;
template<class T, class T1> int chkmin(T &x, const T1 &y) { return x > y ? x = y, 1 : 0; }
template<class T, class T1> int chkmax(T &x, const T1 &y) { return x < y ? x = y, 1 : 0; }
const int MAXN = (1 << 18);
const int mod = (int)1e9 + 7;
int pw(int x, int p, int m) {
int r = 1;
while(p) {
if(p & 1) r = r * 1ll * x % m;
x = x * 1ll * x % m;
p >>= 1;
}
return r;
}
int inv(int x, int m) {
return pw(x, m - 2, m);
}
int n, q, p;
struct Fraction {
int a, b;
Fraction(int u = 0, int d = 1) {
a = u;
b = d;
}
Fraction operator*(Fraction oth) {
return Fraction(a * 1ll * oth.a % p, b * 1ll * oth.b % p);
}
int val() {
return a * 1ll * ::inv(b, p) % p;
}
Fraction inv() {
return Fraction(b, a);
}
};
struct expr {
int main_pw, rev;
Fraction coef;
expr() { coef = Fraction(0); main_pw = 0; rev = 0; }
expr(Fraction _coef, int _main_pw, int _rev = 0) {
coef = _coef;
main_pw = _main_pw;
rev = _rev;
}
expr operator*(expr other) {
if(!this->rev) return expr(this->coef * other.coef, this->main_pw + other.main_pw, other.rev);
else {
expr tmp = other.inv();
return expr(this->coef * tmp.coef, this->main_pw + tmp.main_pw, tmp.rev);
}
}
expr inv_f() {
if(this->rev) return expr(this->coef, this->main_pw, this->rev);
else return expr(this->coef.inv(), -this->main_pw, this->rev);
}
expr inv() {
return expr(this->coef.inv(), -this->main_pw, this->rev ^ 1);
}
};
struct persistent_dsu {
int par[MAXN], sz[MAXN];
expr e[MAXN];
bool failed;
int main_val, cnt_comps;
void init(int n) {
for(int i = 0; i < n; i++) {
par[i] = i;
sz[i] = 1;
e[i] = expr(1, 0);
}
failed = false;
cnt_comps = n;
main_val = -1;
}
pair<int, expr> root(int x) {
if(x == par[x]) {
return {x, e[x]};
}
pair<int, expr> abv = root(par[x]);
abv.second = e[x] * abv.second;
return abv;
}
vector<pair<int, int>> snapshots;
void unite(int x, int y, int v) {
if(failed) return;
auto p = root(x);
auto q = root(y);
if(q.first == 0 || (p.first != 0 && sz[p.first] < sz[q.first])) {
swap(p, q);
}
expr ex = p.second, ey = q.second;
x = p.first, y = q.first;
if(x == y) {
// cycle
// ex(x) * ey(y) = main * v
int cf = ex.coef.val() * 1ll * ey.coef.val() % ::p;
int pw_0 = ex.main_pw + ey.main_pw - 1;
if(ex.rev) pw_0--;
else pw_0++;
if(ey.rev) pw_0--;
else pw_0++;
cf = cf * 1ll * inv(v, ::p) % ::p;
if(pw_0 == 0) {
failed = (cf != 1);
} else {
if(pw_0 > 0) cf = inv(cf, ::p), pw_0 *= -1;
assert(pw_0 == -1);
if(main_val == -1) {
main_val = cf;
} else {
failed = (cf != main_val);
}
}
} else {
// make x the root
cnt_comps--;
par[y] = x;
sz[x] += sz[y];
// ex(x) * ey(y) = main * v
// ey(y) = main * v * ex(x).inv
// y = apply |ey.inv_f()| (main * v * ex(x).inv)
e[y] = ey.inv_f() * (expr(Fraction(v), 1) * ex.inv());
snapshots.push_back({x, y});
}
}
void rollback(bool was_failed, int snapshots_sz, int main_val_prv) {
main_val = main_val_prv;
failed = was_failed;
while(SZ(snapshots) > snapshots_sz) {
int x = snapshots.back().first, y = snapshots.back().second;
sz[x] -= sz[y];
par[x] = x;
par[y] = y;
e[x] = expr(1, 0);
e[y] = expr(1, 0);
cnt_comps++;
snapshots.pop_back();
}
}
}d;
void read() {
cin >> n >> q >> p;
}
int ans[MAXN];
vector<pair<int, pair<int, int>>> li[MAXN << 2];
void add(int ql, int qr, pair<int, pair<int, int>> to_add, int l, int r, int idx) {
if(ql <= l && r <= qr) {
li[idx].pb(to_add);
return;
}
int mid = (l + r) >> 1;
if(ql <= mid) add(ql, qr, to_add, l, mid, 2 * idx + 1);
if(mid < qr) add(ql, qr, to_add, mid + 1, r, 2 * idx + 2);
}
void solve(int l, int r, int idx) {
bool was_failed = d.failed;
int sz = SZ(d.snapshots), main_val = d.main_val;
for(auto upd: li[idx]) {
int x = upd.second.first, y = upd.second.second, v = upd.first;
if(y > 0) y += n - 1;
d.unite(x, y, v);
}
if(l == r) {
if(d.failed) ans[l] = 0;
else ans[l] = pw(p - 1, d.cnt_comps - (d.main_val != -1), mod);
} else if(!d.failed) {
int mid = (l + r) >> 1;
solve(l, mid, 2 * idx + 1);
solve(mid + 1, r, 2 * idx + 2);
}
d.rollback(was_failed, sz, main_val);
}
void solve() {
d.init(2 * n - 1);
// We can notice that if the determinants for all 2x2 matrices are 0, the whole matrix will be curious
// Then we can easily convert the constraints to the type a[0][0] * v = a[x][0] * a[0][y]
map<pair<int, pair<int, int>>, int> last;
for(int i = 0; i < q; i++) {
pair<int, pair<int, int>> curr;
cin >> curr.second.first >> curr.second.second >> curr.first;
curr.second.first--;
curr.second.second--;
if(!last.count(curr)) {
last[curr] = i;
} else {
add(last[curr], i - 1, curr, 0, q - 1, 0);
last.erase(curr);
}
}
for(auto it: last) {
add(it.second, q - 1, it.first, 0, q - 1, 0);
}
solve(0, q - 1, 0);
for(int i = 0; i < q; i++) {
cout << ans[i] << endl;
}
}
int main() {
ios_base::sync_with_stdio(false);
cin.tie(nullptr);
read();
solve();
return 0;
}
This is the first time I have handled the full responsibilty of setting a problem on Codechef (earlier my problems were used but I didn’t take part in curation) . I hope the participants liked the ideas involved in the problem, pardon me for the tedious implementation involved . Suggestions are most welcome 
.