ii. Gaussian Elimination, Rank and Inverse

Row reduction solves $Ax=b$, computes rank, inverse, null space and solution counts in $O(n^{2}m)$; over ${\mathbb{F}}_2$ the rows are bitsets and the cost drops by 64. mint from Modular Arithmetic; determinants are in Determinants.

Variant	What it gives
Gaussian elimination (row echelon form) + back substitution	solution of $Ax=b$; determinant = product of pivots, sign flipped per swap
Gauss–Jordan (reduced row echelon form)	solution read off directly; inverse from $[A\mid I]$; several right-hand sides at once; pivot and free columns expose a null-space basis
LU decomposition	records the multipliers: $PA=LU$, then each new $b$ costs $O(n^2)$
rank	number of pivots
over ${\mathbb{F}}_p$	exact; solution count $p^{m-\mathrm{rank}}$
over ${\mathbb{F}}_2$	XOR systems with bitset rows, $O(n^{2}m/64)$; the XOR basis is the same computation
over doubles	partial pivoting (largest $|a|$ in the column) for stability, compare with $\epsilon$
simplex	Gauss–Jordan pivots generalised to inequalities (Linear Programming)

Condition after reducing $[A\mid b]$ ($m$ unknowns)	Solutions
$\mathrm{rank}(A)<\mathrm{rank}([A\mid b])$ (a row $0=c\ne 0$)	none
$\mathrm{rank}(A)=\mathrm{rank}([A\mid b])=m$	unique
$\mathrm{rank}(A)=\mathrm{rank}([A\mid b])<m$	$m-\mathrm{rank}$ free variables: infinitely many over $\mathbb{R}$, exactly $p^{m-\mathrm{rank}}$ over ${\mathbb{F}}_p$ (usually what the problem asks)

int gaussDouble(vector<vector<double>> a, vector<double>& x) {   // a: n x (m+1) augmented; 0 none, 1 unique, 2 infinitely many
    const double EPS = 1e-9; int n = sz(a), m = sz(a[0]) - 1; vector<int> where(m, -1);
    for (int col = 0, row = 0; col < m && row < n; col++) {
        int sel = row; for (int i = row; i < n; i++) if (fabs(a[i][col]) > fabs(a[sel][col])) sel = i;
        if (fabs(a[sel][col]) < EPS) continue;
        swap(a[sel], a[row]); where[col] = row;
        for (int i = 0; i < n; i++) if (i != row) { double f = a[i][col] / a[row][col]; for (int j = col; j <= m; j++) a[i][j] -= a[row][j] * f; }
        row++;
    }
    x.assign(m, 0);
    for (int j = 0; j < m; j++) if (where[j] != -1) x[j] = a[where[j]][m] / a[where[j]][j];
    for (int i = 0; i < n; i++) { double s = 0; for (int j = 0; j < m; j++) s += x[j] * a[i][j]; if (fabs(s - a[i][m]) > EPS) return 0; }
    for (int j = 0; j < m; j++) if (where[j] == -1) return 2;
    return 1;
}
 
struct LinearSolution { bool ok = false; vector<mint> x; vector<vector<mint>> kernel; int rank = 0; };   // all solutions: x + span(kernel)
LinearSolution solveMod(vector<vector<mint>> a, vector<mint> b) {   // n equations, m unknowns, prime modulus
    int n = sz(a), m = n ? sz(a[0]) : 0; vector<int> where(m, -1); int row = 0;
    for (int i = 0; i < n; i++) a[i].push_back(b[i]);
    for (int col = 0; col < m && row < n; col++) {
        int sel = -1; for (int i = row; i < n; i++) if (a[i][col] != mint(0)) { sel = i; break; }
        if (sel < 0) continue;
        swap(a[sel], a[row]); where[col] = row;
        mint inv = a[row][col].inv(); for (int j = col; j <= m; j++) a[row][j] *= inv;
        for (int i = 0; i < n; i++) if (i != row && a[i][col] != mint(0)) { mint f = a[i][col]; for (int j = col; j <= m; j++) a[i][j] -= f * a[row][j]; }
        row++;
    }
    LinearSolution res; res.rank = row;
    for (int i = row; i < n; i++) if (a[i][m] != mint(0)) return res;   // 0 = c
    res.ok = true; res.x.assign(m, 0);
    for (int j = 0; j < m; j++) if (where[j] != -1) res.x[j] = a[where[j]][m];
    for (int f = 0; f < m; f++) if (where[f] == -1) {      // one kernel vector per free column
        vector<mint> v(m, 0); v[f] = 1;
        for (int j = 0; j < m; j++) if (where[j] != -1) v[j] = -a[where[j]][f];
        res.kernel.push_back(v);
    }
    return res;
}
 
template <class T> vector<vector<T>> matrixInverse(vector<vector<T>> a) {   // Gauss-Jordan on [A | I]; {} if singular
    int n = sz(a);
    for (int i = 0; i < n; i++) a[i].resize(2 * n), a[i][n + i] = T(1);
    for (int c = 0; c < n; c++) {
        int sel = c;
        for (int i = c + 1; i < n; i++) {
            if constexpr (is_floating_point_v<T>) { if (abs(a[i][c]) > abs(a[sel][c])) sel = i; }
            else { if (a[sel][c] == T(0)) sel = i; }
        }
        if constexpr (is_floating_point_v<T>) { if (abs(a[sel][c]) < 1e-9) return {}; }
        else { if (a[sel][c] == T(0)) return {}; }
        swap(a[sel], a[c]);
        T inv = T(1) / a[c][c]; for (int j = c; j < 2 * n; j++) a[c][j] *= inv;
        for (int i = 0; i < n; i++) if (i != c && a[i][c] != T(0)) { T f = a[i][c]; for (int j = c; j < 2 * n; j++) a[i][j] -= f * a[c][j]; }
    }
    for (auto& row : a) row.erase(row.begin(), row.begin() + n);
    return a;
}

Proof: elimination, the three cases, kernels and the inverse

Row operations (swap, scale by a non-zero factor, add a multiple of one row to another) are invertible and only recombine equations, so they preserve the solution set and the row space; hence they preserve the rank. Reduced form: by induction on the processed columns, each pivot column contains a single non-zero (normalised to 1 over ${\mathbb{F}}_p$) in its pivot row and zeros elsewhere: clearing column $c$ adds multiples of the pivot row, which is zero in every earlier pivot column. In the final form the non-zero rows are independent (each has a 1 where all others have 0), so the rank is the number of pivots. Three cases: a zero row of $A$ with non-zero right side is the contradiction $0=c$; otherwise every assignment of the free variables extends uniquely (each pivot variable is its right side minus the free part), so the solutions are $x_0+\mathrm{span}(\text{kernel})$, an affine space of dimension $m-\mathrm{rank}$, with $p^{m-\mathrm{rank}}$ points over ${\mathbb{F}}_p$. Kernel vectors: setting free variable $f$ to 1 and the other free variables to 0 forces pivot variable $j$ to $-a_{\mathrm{where}(j),f}$; these $m-\mathrm{rank}$ vectors are independent (each has its own free coordinate). Inverse: the row operations that reduce $A$ to $I$ multiply $[A\mid I]$ on the left by some $E$ with $EA=I$, so the right block is $E=A^{-1}$; a column with no pivot means $\mathrm{rank}<n$, i.e. singular. Doubles: partial pivoting divides by the largest available entry, bounding the multipliers by 1 in absolute value, which limits error growth; consistency is rechecked by substituting. LU: eliminating below the pivot multiplies by a unit lower-triangular matrix holding the multipliers; the product of these inverses is again unit lower triangular, so $PA=LU$. All functions were checked against brute force over ${\mathbb{F}}_5$ and ${\mathbb{F}}_7$ systems (exhaustive solution counts) and random real systems.

Rank and its uses

Problem	Rank answer
distinct subset XORs of $n$ numbers	$2^{\mathrm{rank}}$
subsets with a given XOR	$2^{n-\mathrm{rank}}$ or 0
XOR equations (light switches, Lights Out)	${\mathbb{F}}_2$ elimination
perfect matching exists (general graph)	Tutte matrix with random values has rank $n$ (General Matching)
maximum matching size	$\frac{1}{2}$ rank of the random Tutte matrix (Lovász)
dimension of a span, independence of vectors	rank (= number of vectors iff independent)
solutions of a modular linear system	$p^{m-\mathrm{rank}}$
minimum generators of $({\mathbb{Z}}_2{)}^k$-subgroup	rank

Numerical rank. Over the reals “is this entry zero” is meaningless in floating point: use the singular value decomposition and count singular values above a tolerance. In contests, work over ${\mathbb{F}}_p$ instead (exact). Sparse systems: elimination causes fill-in; use Wiedemann or Lanczos, which only need matrix–vector products, $O(n\cdot \mathrm{nnz})$ time and $O(n)$ extra space.

Over GF(2)

int rankF2(vector<ull> a) {                                // rows as 64-bit masks
    int n = sz(a), rk = 0;
    for (int col = 63; col >= 0 && rk < n; col--) {
        int piv = -1; for (int r = rk; r < n; r++) if (a[r] >> col & 1) { piv = r; break; }
        if (piv < 0) continue;
        swap(a[rk], a[piv]);
        for (int r = 0; r < n; r++) if (r != rk && (a[r] >> col & 1)) a[r] ^= a[rk];
        rk++;
    }
    return rk;
}
 
template <int N> int solveXor(vector<bitset<N>> a, int m, bitset<N>& x) {   // column m holds b; returns -1 or #free vars
    int n = sz(a), row = 0; vector<int> where(m, -1);
    for (int col = 0; col < m && row < n; col++) {
        int piv = -1; for (int i = row; i < n; i++) if (a[i][col]) { piv = i; break; }
        if (piv < 0) continue;
        swap(a[piv], a[row]); where[col] = row;
        for (int i = 0; i < n; i++) if (i != row && a[i][col]) a[i] ^= a[row];
        row++;
    }
    for (int i = row; i < n; i++) if (a[i][m]) return -1;
    x.reset(); for (int j = 0; j < m; j++) if (where[j] != -1) x[j] = a[where[j]][m];
    return m - row;
}
 
template <int N> vector<bitset<N>> inverseMod2(vector<bitset<N>> a, int n) {   // {} if singular, O(n^3/64)
    vector<bitset<N>> inv(n); for (int i = 0; i < n; i++) inv[i][i] = 1;
    for (int i = 0; i < n; i++) {
        int piv = i; while (piv < n && !a[piv][i]) piv++;
        if (piv == n) return {};
        swap(a[i], a[piv]), swap(inv[i], inv[piv]);
        for (int j = 0; j < n; j++) if (j != i && a[j][i]) a[j] ^= a[i], inv[j] ^= inv[i];
    }
    return inv;
}
 
template <int N> int rankMod2(vector<bitset<N>> a, int m) {
    int n = sz(a), rank = 0;
    for (int col = 0; col < m && rank < n; col++) {
        int piv = rank; while (piv < n && !a[piv][col]) piv++;
        if (piv == n) continue;
        swap(a[rank], a[piv]);
        for (int j = rank + 1; j < n; j++) if (a[j][col]) a[j] ^= a[rank];
        rank++;
    }
    return rank;
}
 
int lightsOut(int R, int C, const vector<string>& g) {     // min presses to turn all '1' off (press toggles cell + 4 neighbours), -1 if impossible
    const int N = 401; int n = R * C; vector<bitset<N>> a(n);
    for (int r = 0; r < R; r++) for (int c = 0; c < C; c++) {
        int i = r * C + c; a[i][n] = g[r][c] == '1';
        int dr[] = {0, 1, -1, 0, 0}, dc[] = {0, 0, 0, 1, -1};
        for (int d = 0; d < 5; d++) { int rr = r + dr[d], cc = c + dc[d]; if (rr >= 0 && rr < R && cc >= 0 && cc < C) a[i][rr * C + cc] = 1; }
    }
    int row = 0; vector<int> where(n, -1);                 // reduce, then enumerate the free variables (2^free)
    for (int col = 0; col < n && row < n; col++) {
        int piv = -1; for (int i = row; i < n; i++) if (a[i][col]) { piv = i; break; }
        if (piv < 0) continue;
        swap(a[piv], a[row]); where[col] = row;
        for (int i = 0; i < n; i++) if (i != row && a[i][col]) a[i] ^= a[row];
        row++;
    }
    for (int i = row; i < n; i++) if (a[i][n]) return -1;
    vector<int> fr; for (int j = 0; j < n; j++) if (where[j] == -1) fr.push_back(j);
    int best = INT_MAX;
    for (int mask = 0; mask < (1 << sz(fr)); mask++) {
        bitset<N> x; for (int t = 0; t < sz(fr); t++) x[fr[t]] = mask >> t & 1;
        for (int j = 0; j < n; j++) if (where[j] != -1) { bitset<N> rowv = a[where[j]]; rowv[n] = 0; rowv[j] = 0; x[j] = a[where[j]][n] ^ ((rowv & x).count() & 1); }
        best = min(best, (int)x.count());
    }
    return best;
}

Lights Out. Each cell is an ${\mathbb{F}}_2$ variable (pressed or not; order and repetition don’t matter), each light gives one equation. Rank $=n$: unique solution; rank $<n$: $2^{n-\mathrm{rank}}$ solutions, enumerate the free variables for the fewest presses; inconsistent: impossible. This models the whole family of toggle puzzles. GF(2) determinant and matchings. $\det \equiv \mathrm{per}(\mathrm{mod}2)$, so $\det$ of a bipartite adjacency matrix mod 2 is the parity of the number of perfect matchings; CF 736D: if the count is odd, removing edge $(u,v)$ keeps it odd iff inv[v][u] == 0.

Proof: XOR systems, Lights Out and the matching-parity trick

XOR systems are linear systems over ${\mathbb{F}}_2$, where addition is XOR and the only non-zero scalar is 1, so elimination is “XOR the pivot row into every other row with a 1 in the column”; solveXor returns a particular solution (free variables 0) and the number of free variables, so there are $2^{\text{free}}$ solutions. Lights Out: pressing a cell twice is the identity and presses commute, so a strategy is a subset of cells, and light $i$ ends off iff the number of presses in its closed neighbourhood has the parity of its initial state: one linear equation per light. Every solution is obtained by choosing the free variables and back-substituting each pivot variable from its reduced row (the row contains its pivot and free columns only), so minimising over all $2^{\text{free}}$ choices is exact. Determinant mod 2: $\det A={\sum }_{\sigma }\mathrm{sgn}\sigma {\prod }_i{A}_{i\sigma (i)}$ and $\mathrm{sgn}\sigma \equiv 1(\mathrm{mod}2)$, so $\det \equiv \mathrm{per}$. CF 736D: matchings avoiding the edge $(u,v)$ number $\mathrm{per}(A)-\mathrm{per}(A_{\hat{u}\hat{v}})$ (matchings using it are matchings of the minor). Mod 2 this is $1-C_{uv}$ with $C_{uv}$ the cofactor, and $(A^{-1}{)}_{vu}=C_{uv}/\det A=C_{uv}$. Verified against brute force (exhaustive press sets for grids up to $4\times 4$, matching counts for random bipartite graphs).

Related notes

• XOR Basis
• Determinants
• Linear Programming (simplex)