CUCUMBER - Editorial

[anton_lunyov]

PROBLEM LINK:

Practice
Contest

Author: Hiroto Sekido
Tester: Anton Lunyov
Editorialist: Anton Lunyov

DIFFICULTY:

HARD

PREREQUISITES:

Determinant, Cauchy-Binet formula, Gauss-Jordan elimination, Bitwise operations

PROBLEM:

Getting rid of the story line we have the following problem.
You are given N × N matrices Q₁, Q₂, …, Q_B.
For each a and b such that 1 ≤ a < b ≤ B we construct the matrix C_a,b by the formula
C_a,b[i][j] = Sum[Q_a[i][k] * Q_b[j][k] : 1 ≤ k ≤ N] for 1 ≤ i, j ≤ N.
For each permutation p = {p[1], p[2], …, p[N]} of integers from 1 to N, inclusive, we consider the sequence
Tr§ = {C_a,b[i][p[i]] : 1 ≤ i ≤ N}.
Denote by CNT(a, b) the number of permutations p, for which the sequence Tr§ has at least one odd element.
We call the pair (a, b) good if CNT(a, b) is odd.
Now the task is to calculate the total number of good pairs.

QUICK EXPLANATION:

In the case N = 1 the answer is M * (M − 1) / 2 where M be the number of odd values among {Q_a[1][1] : 1 ≤ a ≤ B}.

So let N be at least 2.

Since permanent of the matrix coincides with its determinant modulo 2 (see here) then the pair (a, b) is good if and only if det(D_a,b) mod 2 = 1, where D_a,b[i][j] = (C_a,b[i][j] + 1) mod 2.

Let’s add column composed of ones to each matrix Q_c and assume that each Q_c is actually N × (N+1) matrix having last column composed of ones. Then D_a,b = (Q_a)^T * Q_b mod 2, where A^T is the matrix transpose to A.

By Cauchy-Binet formula det(D_a,b) mod 2 = Sum[det(Q_a,c) * det(Q_b,c) : 1 ≤ c ≤ N+1] mod 2,
where Q_a,c is the matrix obtained from Q_a by deleting its c-th column.

Introducing numbers V_a = Sum[(det(Q_a,c) mod 2) * 2^c-1 : 1 ≤ c ≤ N+1] we can rewrite det(D_a,b) mod 2 as bitCount(V_a & V_b) mod 2, where X & Y is a bitwise “and” of the numbers X and Y, and bitCount(X) is the number of ones in binary of X. Using built-in function __builtin_parityll() of GNU C++ Compiler or doing some precalc we can find det(D_a,b) mod 2 in O(1) for each (a, b) after precomputing det(Q_a,c) mod 2 for all a and c.

Applying Gauss-Jordan elimination to Q_a modulo 2 we can assume that it is represented in reduced row echelon form since elementary row operations do not change any of the determinants det(Q_a,c) mod 2. If rank of Q_a is less than N then all determinants det(Q_a,c) mod 2 are zero. Otherwise Q_a could be obtained from identity matrix by inserting some column G = {G[1], …, G[N]} at some position (let it be k-th column of Q_a). Then
det(Q_a,c) mod 2 = G[c] for 1 ≤ c < k,
det(Q_a,c) mod 2 = 1 for c = k, and
det(Q_a,c) mod 2 = 0 for k < c ≤ N+1.

Since Gauss-Elimination could be performed in O(N²) time using bitwise operations we get the solution with complexity O(B * N² + B²).

EXPLANATION:

##The case N = 1
In this case the matrix C_a,b has only one element and it equals to Q_a[1][1] * Q_b[1][1].
We have only one permutation p = {1} in this case.
Hence, the pair (a, b) is good if and only if C_a,b[1][1] is odd, which is equivalent to the oddness of both numbers Q_a[1][1] and Q_b[1][1].
This observation allows to come up with O(B) solution for this special case. Let M be the number of odd values among {Q₁[1][1], Q₂[1][1], …, Q_B[1][1]}. Since every good pair corresponds to some pair of odd numbers in this list, then the number of good pairs is M * (M − 1) / 2 and we are done.

From now on we assume that N ≥ 2.

##Reduction to the determinant
For the given pair (a, b) we construct another matrix such that its determinant is odd if and only if the pair (a, b) is good. Namely, let D_a,b[i][j] = (C_a,b[i][j] + 1) mod 2. Then for each permutation p = {p[1], p[2], …, p[N]} the product
Product[D_a,b[i][p[i]] : 1 ≤ i ≤ N] is zero if and only if we have at least one odd element among C_a,b[i][p[i]].
Indeed, when all they are even then D_a,b[i][p[i]] = 1 for all i and the product equals to 1.
On the other hand, when we have at least one odd element, say C_a,b[i][p[i]], then D_a,b[i][p[i]] = 0 and the product is zero.

Hence, the sum of these products over all permutations equals to the number of permutations p, for which the sequence Tr§ has all even elements, which is equal to N! − CNT(a, b).

Now the trick is that modulo 2 this sum coincide with the determinant of D_a,b which is denoted as det(D_a,b). Indeed, determinant also equals to the sum of these products over all permutations but each product is taken with some sign. But -1 and 1 has the same parity hence we can ignore signs when consider det(D_a,b) mod 2.

Since N! is even when N ≥ 2 (that is one of the reasons to handle the case N = 1 separetely) then we get det(D_a,b) mod 2 = CNT(a,b) mod 2.

Therefore, the pair (a, b) is good if and only if det(D_a,b) mod 2 = 1.

This observation allows to get a solution with complexity O(B² * N³) but this of course is too slow. We can optimize it to O(B² * N²) using some bitwise operations but this still will be too slow.

##Representing Da,b as a product of matrices
Note that D_a,b[i][j] = (C_a,b[i][j] + 1) mod 2
= (Sum[Q_a[i][k] * Q_b[j][k] : 1 ≤ k ≤ N] + 1) mod 2
= Sum[Q_a[i][k] * Q_b[j][k] : 1 ≤ k ≤ N+1] mod 2
if we put Q_c[i][N+1] = 1 for all c and i.
In other words, we add column composed of ones to each matrix Q_c.

So let’s assume that each Q_c is actually N × (N+1) matrix having last column composed of ones.

This observation allows us to look at D_a,b from another point of view.
Namely, D_a,b = (Q_a)^T * Q_b mod 2, where A^T denotes the matrix transpose to the matrix A and star here is the standard matrix multiplication.

##Cauchy and Binet will save us
Now using Cauchy-Binet formula and the property det(A) = det(A^T) we get
det(D_a,b) mod 2 = Sum[det(Q_a,c) * det(Q_b,c) : 1 ≤ c ≤ N+1] mod 2,
where Q_a,c is the matrix obtained from Q_a by deleting its c-th column.

This observation gives O(B * N³ + B²) solution in the following way. We preprocess each matrix Q_a by computing all the determinants det(Q_a,c) mod 2. Each determinant could be computed in O(N²) time using Gaussian elimination and bitwise operations (see below). Since we need N+1 determinants for each matrix this step takes O(B * N³) operations in all.

Introducing numbers
V_a = Sum[det(Q_a,c) mod 2 * 2^c-1 : 1 ≤ c ≤ N+1]
we can rewrite det(D_a,b) mod 2 as bitCount(V_a & V_b) mod 2, where X & Y is a bitwise “and” of the numbers X and Y, and bitCount(X) is the number of ones in binary of X. Since N ≤ 60 we can use 64-bit integer type to store V_a. If you are using GNU C++ compiler then to calculate bitCount(X) mod 2 you can use built-in functions __builtin_parityll() that returns exactly this number or more well-known __builtin_popcountll() that returns the number of ones in binary. They seem to be implemented very fast and are sufficient to keep the double loop over all pairs (a, b) fast enough. See author’s solution as a reference. Another way is to precompute bitCount values for numbers up to some small bound (2¹⁶ or 2²¹) and then divide number into several parts of almost equal length in binary and compute bitCount as a sum of bitCounts of these parts using precomputed tables.

But such solution also appears to be slow due to O(B * N³) part (only mugurelionut was managed to optimize this enough to get AC).

##Gauss-Jordan elimination
The final step is to figure it out how to get the sequence {det(Q_a,c) mod 2 : 1 ≤ c ≤ N+1} in O(N²) time.
For this we apply Gauss-Jordan elimination to the matrix Q_a modulo 2.
That is, we reduce this matrix to reduced row echelon form using only elementary row operations performing all operations in the field Z₂.
Note that each elementary row operation does not change any of the determinants det(Q_a,c) mod 2.
Hence let’s assume that matrix Q_a is already in its reduced row echelon form modulo 2.

At first we note that when rank of the matrix Q_a mod 2 is less than N than all these determinants are zero.
When rank is N then its reduced row echelon form could be like this

100100
010000
001100
000010
000001

So it is obtained from identity matrix by inserting exactly one column. Moreover, if we insert the column G = {G[1], …, G[N]} to be the k-th column then we should have G[i] = 0 for i ≥ k so that it will not ruin row echelon form, while values G[1], …, G[k-1] could be arbitrary. The inserted column is bold in the example above.

Now the funny thing is that
det(Q_a,c) mod 2 = G[c] for 1 ≤ c < k,
det(Q_a,c) mod 2 = 1 for c = k, and
det(Q_a,c) mod 2 = 0 for k < c ≤ N+1.

To see this we at first note that if 1 ≤ i₁ < … < i_L < k are positions of all ones in G then G is the sum of columns of the matrix Q_a with indices i₁, … i_L. Hence if c is not equal to any of the i₁, … i_L and k then in the matrix Q_a,c we have the column G and it is the linear combination of other columns of Q_a,c. Hence det(Q_a,c) mod 2 = 0 in this case. When we delete k-th column from Q_a we simply get the identity matrix; therefore det(Q_a,k) mod 2 = 1. Finally if c equals to one of the i₁, … i_L then all columns in Q_a,c are linearly independent since column G can’t represented as sum of other columns; therefore det(Q_a,c) mod 2 = 1 in this case. So we prove all formulas above.

The Gauss-Jordan elimination modulo 2 could be performed in O(N²) time using bitwise operation. For this we need represent each row of the matrix as 64-bit integer with binary form equal to this row. Then each pivoting could be made using “xor” operation in O(1) time.

Thus we get the solution with complexity O(B * N² + B²) which is fast enough for this problem.

AUTHOR’S AND TESTER’S SOLUTIONS:

Author’s solution can be found here.
Tester’s solution can be found here.

RELATED PROBLEMS:

Alkhwarizm 2012 - Ignite Flower - IGNITE

[triplem]

It’s going to take me a while to think about the problem from a different viewpoint to see exactly how our algorithms match up, but in the meantime - fantastic problem. I learned many unexpected and cool facts about matrices mod 2 in the process.

[anton_lunyov]

I encourage you to repost your brief explanation here now,
since comments to problem are always poorly formatted.

I also find the problem fantastic.

BTW, the last step was discovered in 2 days before the contest.

I told Hiroto:
“I feel we can make O(N^3) to be O(N^2) but I have troubles with it”
and he then realized how to do this

[shimanshu]

@anton_lunyov how do you actually test a solution?

[anton_lunyov]

What do you mean? What solution?

[shimanshu]

i mean the code which we submit during contest.

[triplem]

OK, I’m going to attempt to write up my solution. Hopefully someone may understand some of it, even though I expect it will get a little complicated towards the end.

N=1 is an easy special case, so assume N>=2. We are given B matrices of size NxN and want to calculate how many pairs (X,Y) have det(XY^T+U)=1, where U is a matrix filled with 1s (all work is done modulo 2).

A couple of definitions just to make things easier to understand:

1. Given a vector v, v is odd if it has an odd number of 1s, and even if it has an even (> 0) number of 1s.

2. Given a vector v and matrix, v is odd in M if ∃ an odd w : Mw = v, and v is even in M if ? an even w : Mw = v. Note that a vector may be both odd an even in M.

We know that the determinant of a matrix is 0 iff it can be multiplied to a nonzero vector to get 0.

Note that Uv = 0 if v has an even number of 1s, and 1 (a vector of 1s) otherwise. Let’s call v an odd vector if it has an odd number of 1s, and an even vector if it has an even (>0) number of ones.

det(XY^T+U)=0
⇔ ∃v : v!=0 and (XY^T+U)v = 0
⇔ ∃v : v!=0 and XY^Tv + Uv = 0
⇔ ∃v : (v is even and XY^Tv = 0) or (v is odd and XY^Tv = 1)
⇔ ∃w : (Xw = 0 and w is 0 or even in Y^T) or (Xw = 1 and w is odd in Y^T)

So we’re going to need to perform two tasks:

1. Quickly tell whether a vector is even or odd (or both, or neither) in each given (transposed) matrix.
2. Find vectors with Xw = 0 and Xw = 1 in each given matrix.

Now, let’s first look at a special case - if ?v : v is even and Y^Tv = 0, then clearly XY^Tv = 0 for any X. So if we can find such buildings, we can ignore them entirely for the rest of the problem (and therefore also ignore the ‘w is 0’ case in the last step above.

Let’s start by looking at some particular matrix Y^T, augment it with the identity matrix on the right, and perform Gauss-Jordan elimination, just as we would do to calculate the inverse of Y^T. (Indeed, if Y^T is invertible, the right half will end up as its inverse).

If this results in two or more zero rows on the left, that means there are two independent vectors v and w with Y^Tv = Y^Tw = 0. One of v, w, and v+w will be even - so we can immediately throw out this matrix. From now on, we can assume that there is at most one zero row with everything else looking like the identity matrix (ie, at most one set of dependent vectors).

Note also that this therefore will be true about matrix Y too. In particular, if Y^T has no zero rows, there is precisely one vector with Yv = 1, and if Y^T has one zero row, there is precisely one nonzero vector with with Yv = 0.

The left half of the reduced matrix will either be the identity matrix, or have one zero row with the corresponding column containing some additional ones. For example:

100100|101111
010100|000110
001000|011001
000000|001100
000010|010100
000001|100011

OK, now let’s solve our two problems.

1. Here we’re trying to find, given v, the parity of solutions for w to Y^Tw = v. The matrix we have just calculated gives us solutions for w in terms of sums of elements of v, plus possibly a ‘free’ variable occurring in some positions - ie in this case:

w1 = v1 + v3 + v4 + v5 + v6 + X
w2 = v4 + v5 + X
w3 = v2 + v3 + v6
w4 = X
w5 = v2 + v4
w6 = v1 + v5 + v6

and 0 = v3 + v4

In this case, we have found an even number of vectors adding to 0, so can throw out the matrix like before.

Otherwise, if there is a zero row, then we have found a way of determining if w does not exist. In this case, let z = (0 0 1 1 0 0) - then w exists iff v.z = 0. We store z for later use.

If w does exist, then we can add together all of the formulae above to calculate the parity of w in terms of elements of v, and possibly a free vector.

In this case, w1 + … + w6 = v4 + v5 + v6 + X.

If there is an X term, we can set X to either 0 or 1, and therefore w can be both even and odd, like in this case.

If there is no X term, then we have found a vector which, when dotted with v, tells the parity of w, which we also store.

2. Rather than doing another row reduction on Y, we can solve Yv = 0 and Yv = 1 from this matrix itself.

If there are no zero rows, the right half of this matrix is the inverse of Y^T, which is the transpose of the inverse of Y, so we can multiply the transpose by a vector of 1s to find v with Yv = 1.

If there is a zero row, it’s not too hard to see (but I’ve run out of stamina to explain :P) that the right half of this row is a vector v with Yv = 0.

FINALLY - process all matrices in this manner. For each pair X,Y, find the 0 or 1 vector corresponding to X; and check to see whether this vector is even, odd, neither, or both in Y^T. If we get the result we were looking for in that last step mentioned earlier, add one to the final answer. This just involves a couple of dot products of already calculated vectors which can be processed in virtually constant time.

[triplem]

Have only skimmed the editorial so far but the steps sound similar generally, just with the way things are stored in matrices shuffled around a bit and looked at from a slightly different perspective.

[damians]

Probably one of the best problems I’ve ever solved. At least when we consider “mathy” problems. Thanks to Hiroto, Anton and all the staff!

[iscsi]

I didn’t solve this problem, but I think I find a solution, but I very afraid to implement it.

So my solution is little bit similar with triplem solution. I’m using his notations so we would like to decide det(XY^T+U)=1, where U is a matrix filled with 1s, but instead of calculate determinant I calculate the rank, because det(XY^T+U)=1 is equivalent rank(XY^T T+U)=N. If this matrix has a full rank than rank(XY^T)=N-1 or rank(XY^T)=N. So I had to check 3 cases (because rank(AB)<=min(rank(A),rank(B))):

1)rank(X)=N, rank(Y)=N

So we would like to detect that exists an v!=0 v(XY^T+U)=0. If this happen than v parity is odd, and v(XY^T)=(1,…,1).

Let u=(1…1)(Y^T)^-1 and v=u(X^T)^-1. But we don’t need to calculate v we need just the v parity.

Let x^p=(x_1,…,x_N) where x_i=(X^T)^-1 matrix i-th column parity, so v=u*x^p.

Rank(XY^T+U)=N if u*x^p parity even.

2)rank(X)=N-1, rank(Y)=N

image(X)={a: where exists z zX=a}

kernel(X)={a: where aX=0} (this set has exactly one item…) x⁺=kernel(X).

x⁺=kernel(XY^T)=kernel(X) if x parity even than x⁺(XY^T+U)=0 else x⁺(XY^T+U)=(1,…,1). We need check there exists y⁺!=x⁺ where y⁺(XY^T+U)=(1,…,1). this means y⁺(^)=(1,…,1) because just x⁺ in kernel(XY^T).
Let u=(1…1)*(Y^T)^-1. If u in image(X) this is equivalent u is orthogonal kernel(X^T).

Rank(XY^T+U)=N if kernel(X) item parity odd, and (1…1)*(Y^T)^-1 is orthogonal kernel(X^T).

3)rank(X)=N-1, rank(Y)=N-1

similarly the previous case: x⁺=kernel(XY^T)=kernel(X) parity odd.

y⁺=kernel(Y^T) need to be orhogonal to image(X) because than would be rank(XY^T) < N-1 than rank(XY^T+U)< N, but image(X) is orthogonal kernel(X^T) so kernel(X^T)=kernel(Y^T) .

doesn’t exist u(XY^T)=(1,…,1) this means image(Y^T) doesn’t contain (1,…,1) this can be solve with gauss elimination.

Rank(XY^T+U)=N if kernel(X) item parity odd, kernel(X^T)=kernel(Y^T) and image(Y^T) doesn’t contain (1,…,1).

I will try to implement this for practice… it seems to me, this algorithm is O(BN²+B²).

[mugurelionut]

I will explain how I optimized my O(B * N^3 + B^2) solution to get AC. So… I won’t get into the other details of the problem (they are well explained in the editorial) and I’ll skip to the fact that for each building a we need to compute N+1 determinants. We have the matrix Qa which has N rows and N+1 columns (column N+1 is all 1s) and we need to compute the determinants of the matrices Qa,c, obtained by removing column c from Qa (1<=c<=N+1).

Computing each determinant from scratch is too slow and will get TLE (the time complexity of this approach is O(N^3)). Let’s define a value K (which would ideally be around N/2, but I used K=32 for N>32 and K=0 for N<=32). Let’s consider the process of computing the determinants of the matrices Qa,c (K+1<=c<=N+1). For each such matrix, the first K steps of the Gaussian elimination algorithm (which will bring the matrix to the row echelon form) will be exactly the same (step i consists of bringing a 1 to row i and column i and reducing all the rows below row i). Thus, we start running the Gaussian elimination algorithm on the matrix Qa. Once we reach step c (K+1<=c<=N) we will “take a snapshot” of the rows c,c+1,…,N (i.e. we copy them separately) - let’s call this snapshot©. Within snapshot© we will also swap column N+1 with column c (that’s the same as moving column c to position N+1, which won’t contribute to the determinant when using Gaussian elimination without column swapping).

At the end of this determinant computation algorithm we will get the determinant of Qa,N+1 (because column N+1 is always considered to be the “removed” column which doesn’t contribute to the determinant’s value). Then, we take each snapshot© (K+1<=c<=N) and run Gaussian elimination on it starting from step c. Thus, the time complexity for processing snapshot© is only O((N-c+1)^2).

Let’s see what this means for N=60 and K=32. Computing the determinant of Qa,N+1 takes O(N^2) time, but computing the determinants of Qa,c for c=33,34,…,60 will take 28^2 + 27^2 + 26^2 + … + 3^2 + 2^2 + 1^2.

In order to make things fast enough, note that each snapshot© (K+1<=c<=N) contains up to 32 columns - thus, in this case, we do not need to use 64-bit types and operations, and instead we can use 32-bit types (e.g. unsigned int) and operations (e.g. and, xor).

This part showed how to compute the determinants for the matrices Qa,c for K+1<=c<=N+1 in O(N^2 + approximately (N/2)^2 + (N/2-1)^2 + (N/2-2)^2 + … + 1^2). If we compute the sum in the big-O expression we still get O(N^3), but the constant factor before N^3 is very small.

In order to compute the determinants for the matrices Qa,c (1<=c<=K) we will proceed similarly, but starting from the last column. This time we will use Gaussian elimination which will either try to make the secondary diagonal filled with 1s (and 0s below it) or will try to fill the main diagonal as before, but starting from row/column N towards row/column 1. Once we reach a column c<=K (keep in mind that columns are now processed in descending order), we will again take a snapshot, of the rows above c this time (i.e. rows 1,2,…,c). To keep things short, this part will compute the determinants of the matrices Qa,c (1<=c<=K) in O(N^2 + K^2 + (K-1)^2 + (K-2)^2 + … + 2^2 + 1^2) which is again O(N^3) (since K should be considered to be about N/2), with a very small constant factor. Each snapshot contains at most 32 columns, so 32-bit types and operations can be used in this case, too.

[balakrishnan_v]

(1) The price of the pizza as a result of the collaboration between shop j1 of building i1 and shop j2 of building i2 is given by the dot product between the vectors Q[i1,j1] and Q[i2,j2].

(2) The number of ways in which you have at-least one match to be odd is given by N! - Number of ways in which all matches are even.

(3) In order to compute the number of ways in which all matches are even, we first construct a graph of possible pairings between the shops of buildings i1 and i2 such that the pairs of shops result in even pizza price. This condition (even pizza price) is same as saying Q[i1,j1]’ * Q[i2,j2] mod 2 = 0. OR 1- Q[i1,j1]’ * Q[i2,j2] mod 2 = 1.

(4) Thus, the adjacency matrix of the graph of possible pairings has its (j1,j2) element as 1- Q[i1,j1]’ * Q[i2,j2] mod 2. This may be compactly written as O - (Q[i1] mod 2) * (Q[i2] mod 2) where O is a matrix of all ones.

(5) The number of ways in which Cucumber boy pays more than cucumber girl is higher than vice-versa for a given building pair (i1,i2) if the number of perfect matchings where each matching has atleast one odd-priced pizza is odd. The number of perfect matchings for a bipartite graph is given as the permanent of its adjacency matrix. Thus N!-permanent(O - (Q[i1] mod 2) * (Q[i2] mod 2)’) must be odd. Since permanent of a matrix is same as its determinant in GF[2], we want that N!-det(O + Q[i1] * Q[i2] ) must be odd for every compatible pairs of buildings i1,i2.

(6) Next note that O = o * o’ where o is a vector of all ones. This means that det(O+Q[i1] * Q[i2]’) = det(o * o’ + Q[i1] * Q[i2]’). By the Matrix determinant lemma,
det(o*o’ + Q[i1]*Q[i2]’) = det(Q[i1]*Q[i2]’) + < Adj(Q[i1])*o, Adj(Q[i2])*o>.
where < x,y> is the dot product between 2 vectors x and y and Adj(X) is the adjugate of a matrix given by X^(-1)/det(X) for invertible matrices. However, we are not necessarily dealing with invertible matrices, Hence each element of Adj(X) is the determinant of the (N-1 x N-1) matrix obtained by crossing out the row and column corresponding to the element, and then taking a transpose.

(7) The first term det(Q[i1]*Q[i2]’) = det(Q[i1]) * det(Q[i2]) and hence can be computed for each i. However, computing Adj(X) can take O(N^4) which is not very efficient.

(8) In order to compute more efficiently, we need to transform Q[i] appropriately as follows:
Let us do a gaussian row-elimination on Q[i] for every i such that P[i]*Q[i] = R[i] where R[i] is either identity or has the last row to be completely zero (which means rank(Q[i]) <= N-1). Since, P[i] is simply a row-transformation on Q[i], det(P[i])=1.

(9) Thus det(o o’+Q[i1]Q[i2]’) = det(P[i1]) det(o o’+Q[i1]*Q[i2]’) * det(P[i2]’) = det( (P[i1] o) (P[i2] o)’ + R[i1] R[i2]’)

(10) By the determinant lemma, det( (P[i1]*o) * (P[i2]*o)’ + R[i1]*R[i2]’) = det(R[i1]*R[i2]) + < Adj(R[i1]) * P[i1]*o, Adj(R[i2]) * P[i2]*o> = det(R[i1])*det(R[i2]) + < Adj(R[i1]) * P[i1]*o, Adj(R[i2]) * P[i2]*o>.

(11) So for every i, we need to compute det(R[i]) and Adj(R[i]) * P[i]*o. det(R[i]) is already an outcome of the Gaussian elimination (if it is identity, it is one. Otherwise, the determinant is zero).

(12) Adj(R[i]) is identity if R[i] is identity. If R[i] is not identity, R[i] has the the last row to be zero. Thus the first (N-1) rows of Adj(R[i]) are all zeros.

(13) If the N-1^th row of R[i] is also zero (ie rank(R[i])<=N-2), Adj(R[i]) is all zeros (trivially). Otherwise (rank(R[i])=N-1), for each j<=N, we check if removing column j from R[i] renders some row i<=N-1 of the matrix to be all zeros. If so, Adj(R[i])[N,j]=0, otherwise it Adj(R[i])[N,j]=1.

(14) By matrix elimination, we already have P[i] and thus P[i]*o (if P[i] is stored as a 64-bit integer, this is simply its bitcount).

(15) Thus O[i]=Adj(R[i]) * P[i]*o is computed for every i. Since this is a boolean array, we store O[i] are a 64-bit integer.

(16) Then, in O(N^2) time, we may compute the pairity of N! - det(R[i1])*det(R[i2]) - < Adj(R[i1]) * P[i1]*o, Adj(R[i2]) * P[i2]*o> for every pair(i1,i2). Note that the inner product can be written as bit_counts(O[i1]&O[i2]).

[anton_lunyov]

I was thinking on the same approach while trying to move from O(N^3) to O(N^2).
I should implement it even after we figure it out O(N^2) solution.
Maybe we could prevent this solution to get AC
But at least you made a lot of submits and optimizations to accept it.

[mugurelionut]

Well, I was rather sure during the contest that an O(N^2) solution existed for computing all the N+1 determinants. However, rather than trying to find it, I attempted to optimize as much as possible the O(N^3) solution (which eventually worked). If I had run out of optimizations (without getting AC), there would have still been enough time to think about finding an O(N^2) solution (but after getting AC I had no incentive to think about it anymore).

[mugurelionut]

By the way, one other thing I forgot to mention in my comment is that if the rank of the first N columns of Qa is N-2 or less then all the N+1 determinants are 0 (and so I do not use the snapshots at all) - the same thing is mentioned in the editorial (if the rank of Qa is less than N all the determinants are 0). However, from my measurements, this was just a minor optimization (it did not reduce the running time of my solution too much).

[cptcharisma]

The editorials are really fun… It reeks of the hard work put in by the problem setters and esp our veteran @anton_lunyov. A big thanks to you Anton!!!