Skip to content

Gaussian Elimination

Quick Revisit

  • USE WHEN: solving linear systems, computing matrix rank, counting solutions; for max XOR subset see Linear Basis
  • INVARIANT: row operations preserve the solution set; rank equals number of pivot columns
  • TIME: O(n^2 m) over reals, O(n m) over GF(2) | SPACE: O(n m)
  • CLASSIC BUG: not using partial pivoting over reals -- tiny pivots amplify floating-point error
  • PRACTICE: ../12-practice-ladders/08-ladder-mixed.md

Systematic row reduction to solve systems of linear equations -- over reals, modular arithmetic, and GF(2) (XOR). The same framework handles classical algebra problems, matrix rank computation, counting solutions, and the "maximum XOR subset" family that appears constantly at CF 1800+.

Difficulty: [Advanced]Prerequisites: Math Fundamentals, Bit ManipulationCF Rating Range: 1800-2100


Contents

Three implementations — pick before you code

The same row-reduction skeleton has three meaningfully different implementations, with different failure modes. Decide which one applies before writing any code:

VariantCoefficient fieldDivisionFailure modes
Real-valuedR via doubledouble / doubleTiny pivots blow up rounding; needs partial pivoting and an epsilon for "is this zero?"
ModularZ/pZ for prime pModular inverse via FermatPivot must be coprime to p; non-prime modulus breaks the field assumption
GF(2) / XOR{0,1} with XOR as additionNo division needed"Add a multiple of row i to row j" collapses to a[j] ^= a[i]; the only failure is misidentifying the pivot column

The real-valued and modular variants share most of their code but have opposite failure modes — one fails on numerical instability, the other on arithmetic that does not exist. The GF(2) variant uses bitsets and is typically 32-64x faster than either. The dedicated lesson on the GF(2) variant for max-XOR-subset queries is 11-linear-basis-xor.md.


Intuition

The Pain — three equations, three unknowns

You have three equations in three unknowns:

 2x + 1y - 1z =   8
-3x - 1y + 2z = -11
-2x + 1y + 2z =  -3

The textbook says "use substitution." Fine -- isolate x from the first equation, plug into the second, simplify, isolate y, plug into the third... By the third substitution you are drowning in fractions and sign errors. With 4 unknowns it is practically hopeless by hand.

Now imagine the contest version: n up to 500 unknowns, or 60 bitmask equations over XOR. You need a mechanical, automatable procedure -- not cleverness, just a reliable algorithm.

The Key Insight

Systematically eliminate variables one by one using row operations -- the same idea works over reals, modular arithmetic, AND over GF(2) (XOR) for bitmask problems.

Three legal row operations preserve the solution set:

  1. Swap two rows.
  2. Scale a row by a nonzero constant.
  3. Add a multiple of one row to another.

By choosing these operations carefully, we zero out entries below (and optionally above) each pivot, converting the matrix to upper-triangular (or reduced row-echelon) form. Back-substitution then reads off the answer.

Over GF(2) (the field {0,1} with XOR as addition), "add a multiple of row i to row j" simplifies to just XOR-ing two bitmasks. The entire elimination becomes a loop of a[j] ^= a[i] -- blazingly fast.

Where it breaks: the coefficient matrix must be over a field (so division is possible). Over plain integers without a modulus you get rationals -- use double with care, or better, work modulo a prime.

Visual Walkthrough

Example over reals -- 3x3 system:

Augmented matrix [Ab]:

[  2   1  -1 |   8 ]
[ -3  -1   2 | -11 ]
[ -2   1   2 |  -3 ]

Step 1 -- Pivot on column 0 (eliminate x from rows 1, 2):

R1R1+32R0, R2R2+R0:

[  2    1   -1  |   8  ]          pivot = 2
[  0   1/2  1/2 |   1  ]     <--  R1 += (3/2)*R0
[  0    2    1  |   5  ]     <--  R2 += (2/2)*R0

Step 2 -- Pivot on column 1 (eliminate y from row 2):

R2R24R1:

[  2    1   -1  |   8  ]
[  0   1/2  1/2 |   1  ]          pivot = 1/2
[  0    0   -1  |   1  ]     <--  R2 -= 4*R1

Upper-triangular! Back-substitute:

z = 1 / (-1)       = -1
y = (1 - 1/2*z) / (1/2) = 3
x = (8 - 1*y + 1*z) / 2 = 2

Solution: (x,y,z)=(2,3,1). Three steps, no creativity needed.

Step-by-step elimination flow diagram:

  AUGMENTED MATRIX [A | b]         ROW OPERATIONS          RESULT
  +---+---+---+-----+
  | 2 | 1 |-1 |   8 |
  |-3 |-1 | 2 | -11 |        Step 1: Eliminate column 0
  |-2 | 1 | 2 |  -3 |                (kill x in R1, R2)
  +---+---+---+-----+
         |
         | R1 += (3/2)*R0
         | R2 += (2/2)*R0
         v
  +---+-----+-----+-----+
  | 2 |  1  | -1  |   8 |    <-- pivot row (col 0)
  | 0 | 1/2 | 1/2 |   1 |
  | 0 |  2  |  1  |   5 |
  +---+-----+-----+-----+
         |
         | R2 -= 4*R1          Step 2: Eliminate column 1
         v                            (kill y in R2)
  +---+-----+-----+-----+
  | 2 |  1  | -1  |   8 |
  | 0 | 1/2 | 1/2 |   1 |    <-- pivot row (col 1)
  | 0 |  0  | -1  |   1 |
  +---+-----+-----+-----+
         |
         | Back-substitute       Step 3: Read off solution
         v                              (bottom to top)
  z = 1/(-1) = -1
  y = (1 - 1/2*z) / (1/2) = 3
  x = (8 - y + z) / 2 = 2
         |
         v
  SOLUTION: (x, y, z) = (2, 3, -1)

Geometric interpretation -- 3 planes in 3D:

  Each equation defines a plane in 3D space.
  A unique solution = all three planes meet at one point.

             z ^
               |    Plane 1: 2x + y - z = 8
              /|\        _______________
             / | \      /              /
            /  |  \    /   Plane 2    /
           /   *---\--/---- (-3x - y + 2z = -11)
          / * |     \/              /
    -----/*---+-----/\------------/---------> y
        / *   |    /   \         /
       /  *   |   / Plane 3    /
      /   *   |  / (-2x + y + 2z = -3)
     /    *   | /             /
    /     *   |/             /
   /------*---+-------------/
  /       *  /
 x        *          * = intersection point (2, 3, -1)
          *             all three planes pass through here

Example over GF(2) -- XOR basis construction:

Given values: {6,5,3,2} -- find max XOR of any subset.

Binary representations (4 bits):

6 = 0110
5 = 0101
3 = 0011
2 = 0010

Build the basis by inserting each value:

Insert 6 = 0110:
  Highest set bit = bit 2.  basis[2] is empty --> basis[2] = 0110.
  Basis: [ _ , _ , 0110, _ ]

Insert 5 = 0101:
  Highest bit = bit 2.  basis[2] = 0110 --> XOR: 0101 ^ 0110 = 0011.
  Highest bit = bit 1.  basis[1] is empty --> basis[1] = 0011.
  Basis: [ _ , 0011, 0110, _ ]

Insert 3 = 0011:
  Highest bit = bit 1.  basis[1] = 0011 --> XOR: 0011 ^ 0011 = 0000.
  Value reduced to 0 --> 3 is linearly dependent, skip.

Insert 2 = 0010:
  Highest bit = bit 1.  basis[1] = 0011 --> XOR: 0010 ^ 0011 = 0001.
  Highest bit = bit 0.  basis[0] is empty --> basis[0] = 0001.
  Basis: [ 0001, 0011, 0110, _ ]

Max XOR: greedily XOR from highest basis element down:

ans = 0
ans ^= 0110 = 0110  (bit 2 was 0, now 1 -- improves)
ans ^= 0011 = 0101  (bit 1 was 1, now 0 -- skip? no, check: 0110^0011 = 0101 < 0110)

Wait -- greedy rule: XOR with basis[i] only if it increases the current answer (i.e., the answer's bit i is 0):

ans = 0
ans ^= basis[2] = 0110   (bit 2 of ans was 0 --> XOR improves) ans = 0110
ans ^= basis[1] = 0011   (bit 1 of ans is 1   --> skip)        ans = 0110
ans ^= basis[0] = 0001   (bit 0 of ans is 0   --> XOR improves) ans = 0111

Maximum XOR = 01112=7 (achieved by 652=7, or equivalently 652). The basis has rank 3, so 23=8 distinct XOR values are reachable from these 4 numbers.

Worked Example -- Solving a 3x3 system in GF(2) (mod 2 arithmetic):

Solve the following system over GF(2), where + means XOR:

  x1 ^ x2 ^ x3 = 1
  x1      ^ x3 = 0
       x2 ^ x3 = 1

Represent as augmented bitmask matrix (each row is a 4-bit value: x1 x2 x3 | rhs):

  Initial:          Binary bitmasks:
  +----+----+----+-----+       eq[0] = 1 1 1 | 1  = 0b1111
  | x1 | x2 | x3 | rhs |       eq[1] = 1 0 1 | 0  = 0b1010
  +----+----+----+-----+       eq[2] = 0 1 1 | 1  = 0b0111
  |  1 |  1 |  1 |  1  |
  |  1 |  0 |  1 |  0  |
  |  0 |  1 |  1 |  1  |
  +----+----+----+-----+

  Step 1: Pivot on column 0 (x1).
          Pivot row = eq[0]. Eliminate x1 from eq[1]:
          eq[1] ^= eq[0]:  1010 ^ 1111 = 0101

  +----+----+----+-----+
  |  1 |  1 |  1 |  1  |   <-- pivot
  |  0 |  1 |  0 |  1  |   <-- eq[1] ^ eq[0]
  |  0 |  1 |  1 |  1  |
  +----+----+----+-----+

  Step 2: Pivot on column 1 (x2).
          Pivot row = eq[1]. Eliminate x2 from eq[0] and eq[2]:
          eq[0] ^= eq[1]:  1111 ^ 0101 = 1010
          eq[2] ^= eq[1]:  0111 ^ 0101 = 0010

  +----+----+----+-----+
  |  1 |  0 |  1 |  0  |   <-- eq[0] ^ eq[1]
  |  0 |  1 |  0 |  1  |   <-- pivot
  |  0 |  0 |  1 |  0  |   <-- eq[2] ^ eq[1]
  +----+----+----+-----+

  Step 3: Pivot on column 2 (x3).
          Pivot row = eq[2]. Eliminate x3 from eq[0]:
          eq[0] ^= eq[2]:  1010 ^ 0010 = 1000

  +----+----+----+-----+
  |  1 |  0 |  0 |  0  |   <-- eq[0] ^ eq[2]
  |  0 |  1 |  0 |  1  |
  |  0 |  0 |  1 |  0  |   <-- pivot
  +----+----+----+-----+

  Reduced row echelon form! Read off:
    x1 = 0,  x2 = 1,  x3 = 0

  Verify:  0 ^ 1 ^ 0 = 1  OK
           0     ^ 0 = 0  OK
              1 ^ 0 = 1  OK

  Rank = 3 = number of variables --> unique solution.
  All row operations were XOR -- no division, no fractions!

The Invariant

After processing the first k pivot columns, the leading k columns of the working matrix form an upper-triangular (echelon) structure: every row below a pivot has a zero in that pivot's column.

For the XOR basis, the invariant is:

basis[i] is either 0 (empty) or a value whose highest set bit is exactly bit i, and no other basis element has bit i set.

This guarantees each basis vector is "responsible" for exactly one bit position, making greedy queries and rank computation trivial.

When to Reach for This

Trigger phrase in problem statementTechnique
"solve system of linear equations"Gauss over reals / mod p
"XOR basis" / "linear basis"Gauss over GF(2)
"maximum XOR of a subset"XOR basis + greedy
"number of distinct XOR values"2rank
"can we represent x as XOR of a subset?"Insert into basis, check if x reduces to 0
"matrix rank"Count nonzero pivots after elimination
"count solutions mod p"Gauss mod p, free variables = nrank
"minimum cost to make all XOR zero"XOR basis on differences / constraints

Insight Gaps

XOR basis = Gaussian elimination over GF(2). The XOR basis (also called linear basis) is the result of Gaussian elimination over GF(2). Given a set of integers, you insert them one by one: if the current number's highest unprocessed bit matches an existing basis vector, XOR them and continue; if no match, add to basis. The resulting basis has the property that every element of the original set is a XOR combination of basis vectors, and no basis vector is a XOR of others.

The XOR basis enables maximum XOR queries. Given a set of numbers, what is the maximum XOR of any subset? Build the basis, then greedily: starting from the highest bit, if XOR-ing the current result with the basis vector increases the result, do it. The XOR basis makes this O(log(maxvalue)) per query after O(nlog(maxvalue)) build.

Gaussian elimination over other fields. Over R: solve systems of linear equations, compute matrix inverse, find determinant. Over Zp (prime modulus): same operations, but using modular inverse instead of real division. Over GF(2): use XOR. Each field has the same algorithm structure but different arithmetic rules.

What the Code Won't Teach You

The code mechanically zeros out entries below each pivot. What it hides is the geometric meaning: each row operation preserves the solution set. The matrix changes, but the set of solutions stays the same.

GEOMETRIC MEANING -- 2 equations in 2 unknowns:

  2x + y = 5       <- line in the (x,y)-plane
   x - y = 1       <- another line

  After row ops:    x = 2, y = 1  <- intersection point

  The row operations TILT the lines but keep them
  intersecting at the same point.

  Before:               After:
    y                     y
    |  /   /              |      |
    | /  /                |      |
    |/ /                  |------*--- x = 2
    *------- x            |      (2,1)
   (2,1)                  |

Over GF(2), the "geometry" is different -- you're working in a vector space where every coordinate is 0 or 1, and addition is XOR. The XOR basis spans exactly the set of reachable XOR combinations. The code builds this basis; understanding the linear algebra tells you why it works.

ROW REDUCTION STAGES -- 3x3 system:

  Original:           After step 1:        After step 2:
  ┌            ┐      ┌            ┐       ┌            ┐
  │ 2  1  -1│ 8│      │ 2  1  -1│ 8│       │ 2  1  -1│ 8│
  │-3 -1   2│-11│     │ 0  ½   ½│ 1│       │ 0  ½   ½│ 1│
  │-2  1   2│-3│      │ 0  2   1│ 5│       │ 0  0  -1│ 1│
  └            ┘      └            ┘       └            ┘
     augmented          eliminated            upper
     matrix             column 1              triangular

  Back-substitute: z = -1, y = 3, x = 2
XOR BASIS INSERTION -- building basis from {5, 3, 6, 7}:

  5 = 101    basis = [_, _, _]
             bit 2 empty -> basis[2] = 101
             basis = [101, _, _]

  3 = 011    bit 1 empty -> basis[1] = 011
             basis = [101, 011, _]

  6 = 110    bit 2 has 101 -> 110 ⊕ 101 = 011
             bit 1 has 011 -> 011 ⊕ 011 = 000
             -> 6 is redundant (expressible from basis)

  7 = 111    bit 2 has 101 -> 111 ⊕ 101 = 010
             bit 1 has 011 -> 010... bit 1 of 010 is 1
             -> 010 ⊕ 011 = 001
             bit 0 empty -> basis[0] = 001
             basis = [101, 011, 001]  rank = 3

  Reachable XOR values: 2^3 = 8 (all 3-bit numbers!)

The cycle space of a graph is a vector space over GF(2). The code won't teach you this: in any undirected graph, assign each edge a bit. A cycle is a subset of edges where every vertex has even degree -- equivalently, a vector in {0,1}|E| where each vertex's "incidence constraint" is satisfied mod 2. The set of all cycles forms a vector space over GF(2) of dimension |E||V|+c (where c is the number of connected components). You can find a basis for all cycles by building a spanning tree: each non-tree edge creates one fundamental cycle. XOR basis / Gaussian elimination over GF(2) lets you answer questions like "can these edges form a cycle?" or "how many distinct cycles exist?" (2|E||V|+c1, excluding the empty set).

The k-th smallest XOR value needs reduced row echelon form (RREF). The code gives you row echelon form, but for the k-th smallest XOR value achievable from a set, you need the basis in reduced form: each basis vector has a unique leading bit, and that bit is cleared in all other basis vectors. Then each achievable value corresponds to a unique subset of basis vectors, and you can binary-decompose k to select which basis vectors to XOR -- giving you the k-th smallest in O(B) time.

RREF vs REF -- why it matters for k-th smallest XOR:

  Row Echelon Form (REF):       Reduced Row Echelon (RREF):
  basis[2] = 101                basis[2] = 100
  basis[1] = 011                basis[1] = 010
  basis[0] = 001                basis[0] = 001

  To get RREF from REF:         (clear bit 2 in basis[1]: 011 has bit 2=0, ok)
  basis[2] ^= basis[0]          (clear bit 1 in basis[2]: 101 has bit 1=0, ok)
  -> already reduced here        (clear bit 0 in basis[2]: 101^001=100 OK)
                                 (clear bit 0 in basis[1]: 011^001=010 OK)

  With RREF, the k-th smallest XOR (1-indexed):
    k=1: 001 = basis[0]                        -> XOR = 1
    k=2: 010 = basis[1]                        -> XOR = 2
    k=3: 011 = basis[1] ^ basis[0]             -> XOR = 3
    k=4: 100 = basis[2]                        -> XOR = 4
    k=5: 101 = basis[2] ^ basis[0]             -> XOR = 5
    k=6: 110 = basis[2] ^ basis[1]             -> XOR = 6
    k=7: 111 = basis[2] ^ basis[1] ^ basis[0]  -> XOR = 7

  Pattern: binary representation of k selects which basis vectors to XOR.
  This ONLY works with RREF -- with REF, the mapping is not monotone.

The XOR basis has a greedy structure that makes it feel like a different algorithm. The standard XOR basis insertion -- try to insert a number, XOR away the highest set bit using existing basis elements, and either add the result to the basis or discard it -- looks nothing like Gaussian elimination at first glance. But it IS Gaussian elimination, row by row, over GF(2). Each "basis slot" is a pivot position. The greedy "maximum XOR subset" query walks the basis from high bit to low, which is back-substitution in reduced row echelon form. Recognizing this connection lets you apply linear algebra reasoning to XOR problems and vice versa.

XOR BASIS = GAUSSIAN ELIMINATION (side by side):

  Numbers to insert: [6, 3, 5]  (binary: 110, 011, 101)

  ┌──── XOR Basis View ────────┬──── GaussElim GF(2) View ──────┐
  │                            │                                 │
  │ Insert 6 = 110             │ Row 1: [1 1 0]                 │
  │   basis[2] = 110           │   Pivot at col 2                │
  │                            │                                 │
  │ Insert 3 = 011             │ Row 2: [0 1 1]                 │
  │   bit 2 not set, skip      │   Col 2: no match, skip        │
  │   bit 1 set, basis empty   │   Pivot at col 1                │
  │   basis[1] = 011           │                                 │
  │                            │                                 │
  │ Insert 5 = 101             │ Row 3: [1 0 1]                 │
  │   bit 2 set -> XOR basis[2] │   Eliminate col 2: XOR row 1   │
  │   101 ^ 110 = 011          │   [1 0 1] ^ [1 1 0] = [0 1 1] │
  │   bit 1 set -> XOR basis[1] │   Eliminate col 1: XOR row 2   │
  │   011 ^ 011 = 000          │   [0 1 1] ^ [0 1 1] = [0 0 0] │
  │   Zero -> discard           │   Zero row -> rank unchanged     │
  │                            │                                 │
  │ Basis: {110, 011}          │ Rank: 2                         │
  │ Max XOR = 110^011 = 101    │ 2^rank = 4 distinct XOR values │
  └────────────────────────────┴─────────────────────────────────┘

3. Core Concepts

Fields and Row Reduction

Gaussian elimination works over any field -- a set with +, , ×, ÷ (except by zero). The three fields you will use in CP:

FieldElementsAdditionMultiplicationDivision
R (doubles)floating-point+*/
Z/pZ{0,1,,p1}(a+b)%p(a*b)%pa * inv(b) % p
GF(2){0,1}XORAND(1 is self-inverse)

Pivoting

Partial pivoting (swap the row with the largest absolute value in the current column to the pivot position) is essential over reals to avoid dividing by near-zero values. Over Z/p or GF(2), any nonzero element works -- just pick the first one.

Rank, Free Variables, and Solution Count

After elimination the matrix has r nonzero rows (= rank). For an n-variable system:

  • r=n: unique solution.
  • r<n and system is consistent: nr free variables, pnr solutions over Z/p (or 2nr over GF(2)).
  • Inconsistent (a row of form [0 0  0c] with c0): no solution.

Rank-deficient cases -- geometric view in 3D:

  CASE 1: UNIQUE SOLUTION (rank = 3)       CASE 2: NO SOLUTION (rank < 3,
  Three planes intersect at one point.     inconsistent -- parallel planes)

       /\                                        +-----------+
      /  \                                       |  Plane 1  |
     / *  \   <-- intersection point             +-----------+
    /------\                                          |  (gap)
   /   /\   \                                    +-----------+
  /   /  \   \                                   |  Plane 2  |
     /    \                                      +-----------+
    /  P3  \                                          |  (gap)
                                                 +-----------+
   All 3 planes share exactly                   |  Plane 3  |
   one common point.                             +-----------+

                                                 Planes are parallel:
                                                 no common point exists.
                                                 After elimination: a row
                                                 [0 0 0 | c], c != 0.

  CASE 3: INFINITE SOLUTIONS (rank < n, consistent)

       +----------------------------------+
       |         Plane 1                   |
       |     +---------------------------+ |
       |     |       Plane 2             | |
       |     |   *=========*            | |
       |     |   | line of |            | |
       |     |   | solutions|           | |
       |     +---|---------|------------+ |
       +---------|---------|-------------+
                 *=========*
                      ^
                      |
          Two planes intersect in a line;
          third plane contains that line.
          Free variables = n - rank.

XOR Basis (Linear Basis over GF(2))

An XOR basis is a set of at most log2(max_val) numbers such that:

  • Every XOR-reachable value from the original set is also reachable from the basis.
  • The basis vectors are linearly independent over GF(2).
  • basis[i] has its highest set bit at position i.

Insertion is O(logmax_val). Querying max/min XOR, checking membership, and counting distinct values are all O(logmax_val).

GF(2) VECTOR SPACE FOR 3 BITS -- the XOR cube:

  Vertices = all 3-bit numbers. Edge = XOR by a basis vector.

         110 ─────────── 111          Given basis = {101, 011}:
         /|              /|           Reachable from 000:
        / |             / |             000 ⊕ 101 = 101
      010 ─────────── 011  |            000 ⊕ 011 = 011
       |  |            |  |             101 ⊕ 011 = 110
       | 100 ──────────| 101            011 ⊕ 101 = 110 (same)
       | /             | /              000, 101, 011, 110  <- 2^2 = 4 values
       |/              |/
      000 ─────────── 001           NOT reachable: 001, 010, 100, 111
                                    (these form a parallel coset: {001, 100, 010, 111})
  Bold edges show XOR by basis
  vectors from reachable set.       rank = 2 -> 2^2 = 4 reachable values
                                    The reachable set is a SUBSPACE (contains 000).
                                    Any coset x ⊕ subspace has the same size.
XOR BASIS vs RAW SET -- what the basis buys you:

  Raw set S = {13, 10, 6, 3, 5}
  (5 numbers -- 2^5 = 32 possible subsets to try for max XOR?)

  After XOR basis insertion:
  ┌─────────────────────────────────────────────────┐
  │  bit 3: basis[3] = 1101  (13)                   │
  │  bit 2: ────────  (empty, no vector has bit 2   │
  │                    as highest after reduction)   │
  │  bit 1: basis[1] = 0011  (3)                    │
  │  bit 0: basis[0] = 0001  (derived)              │
  └─────────────────────────────────────────────────┘
  rank = 3 -> 2^3 = 8 distinct reachable XOR values

  Max XOR (greedy from high bit):
    ans = 0000
    bit 3: ans ^ 1101 = 1101 > 0000 -> take it. ans = 1101
    bit 1: ans ^ 0011 = 1110 > 1101 -> take it. ans = 1110
    bit 0: ans ^ 0001 = 1111 > 1110 -> take it. ans = 1111

  Max XOR = 1111 = 15.  Found in O(B) = O(4), not O(2^5).

Implementation

Gaussian Elimination over Reals

cpp
#include <bits/stdc++.h>
using namespace std;

// Solves Ax = b where A is n x n augmented into n x (n+1) matrix a.
// Returns number of solutions: 0 (none), 1 (unique), 2 (infinitely many).
// If unique, solution is written into ans.
const double EPS = 1e-9;

int gauss_real(vector<vector<double>>& a, vector<double>& ans) {
    int n = (int)a.size();
    int m = (int)a[0].size() - 1;  // number of variables
    ans.assign(m, 0);

    vector<int> pivot_col(n, -1);
    int col = 0, row = 0;
    for (; col < m && row < n; col++) {
        // Partial pivoting: find row with largest |a[r][col]|
        int best = row;
        for (int r = row + 1; r < n; r++)
            if (abs(a[r][col]) > abs(a[best][col]))
                best = r;
        if (abs(a[best][col]) < EPS) continue;  // skip zero column
        swap(a[row], a[best]);
        pivot_col[row] = col;

        // Eliminate below and above (full Gauss-Jordan)
        for (int r = 0; r < n; r++) {
            if (r == row) continue;
            if (abs(a[r][col]) < EPS) continue;
            double factor = a[r][col] / a[row][col];
            for (int j = col; j <= m; j++)
                a[r][j] -= factor * a[row][j];
        }
        row++;
    }
    int rank = row;

    // Check consistency
    for (int r = rank; r < n; r++)
        if (abs(a[r][m]) > EPS) return 0;  // no solution

    if (rank < m) return 2;  // infinitely many

    // Read off unique solution
    for (int r = 0; r < rank; r++)
        ans[pivot_col[r]] = a[r][m] / a[r][pivot_col[r]];
    return 1;
}

Gaussian Elimination modulo a Prime

cpp
long long binpow(long long b, long long e, long long mod) {
    long long r = 1; b %= mod;
    while (e > 0) {
        if (e & 1) r = r * b % mod;
        b = b * b % mod;
        e >>= 1;
    }
    return r;
}
long long modinv(long long a, long long mod) {
    return binpow(a, mod - 2, mod);
}

// Solves n equations in m unknowns modulo prime p.
// Returns rank. If unique solution, writes it into ans.
int gauss_mod(vector<vector<long long>>& a, vector<long long>& ans,
              long long p) {
    int n = (int)a.size();
    int m = (int)a[0].size() - 1;
    ans.assign(m, 0);

    vector<int> pivot_col(n, -1);
    int col = 0, row = 0;
    for (; col < m && row < n; col++) {
        // Find first nonzero in this column from row downward
        int piv = -1;
        for (int r = row; r < n; r++) {
            if (a[r][col] != 0) { piv = r; break; }
        }
        if (piv == -1) continue;
        swap(a[row], a[piv]);
        pivot_col[row] = col;

        long long inv = modinv(a[row][col], p);
        for (int j = col; j <= m; j++)
            a[row][j] = a[row][j] * inv % p;

        for (int r = 0; r < n; r++) {
            if (r == row || a[r][col] == 0) continue;
            long long factor = a[r][col];
            for (int j = col; j <= m; j++)
                a[r][j] = (a[r][j] - factor * a[row][j] % p + p) % p;
        }
        row++;
    }
    int rank = row;

    // Check consistency
    for (int r = rank; r < n; r++)
        if (a[r][m] != 0) return -1;  // no solution, signal with -1

    // Unique if rank == m
    for (int r = 0; r < rank; r++)
        ans[pivot_col[r]] = a[r][m];
    return rank;
}

XOR Basis (Gaussian Elimination over GF(2))

cpp
struct XorBasis {
    static constexpr int BITS = 60;  // for values up to ~1e18
    long long basis[BITS] = {};
    int rank = 0;

    // Insert val into basis. Returns true if val was independent.
    bool insert(long long val) {
        for (int i = BITS - 1; i >= 0; i--) {
            if (!((val >> i) & 1)) continue;
            if (!basis[i]) {
                basis[i] = val;
                rank++;
                return true;
            }
            val ^= basis[i];
        }
        return false;  // val is linearly dependent
    }

    // Maximum XOR achievable from the basis.
    long long max_xor() const {
        long long ans = 0;
        for (int i = BITS - 1; i >= 0; i--)
            ans = max(ans, ans ^ basis[i]);
        return ans;
    }

    // Maximum XOR starting from an initial value.
    long long max_xor_with(long long init) const {
        long long ans = init;
        for (int i = BITS - 1; i >= 0; i--)
            ans = max(ans, ans ^ basis[i]);
        return ans;
    }

    // Minimum nonzero XOR achievable from the basis.
    long long min_xor() const {
        for (int i = 0; i < BITS; i++)
            if (basis[i]) return basis[i];
        return 0;  // only if basis is empty
    }

    // Can val be represented as XOR of a subset?
    bool contains(long long val) const {
        for (int i = BITS - 1; i >= 0; i--) {
            if (!((val >> i) & 1)) continue;
            if (!basis[i]) return false;
            val ^= basis[i];
        }
        return true;  // val reduced to 0
    }

    // Number of distinct XOR values representable (including 0).
    long long count_distinct() const {
        return 1LL << rank;
    }

    // Merge another basis into this one.
    void merge(const XorBasis& other) {
        for (int i = BITS - 1; i >= 0; i--)
            if (other.basis[i]) insert(other.basis[i]);
    }
};

Gaussian Elimination on Bitmask Equations

For systems of XOR equations (e.g., "given that xa1xa2=bi for each constraint"), represent each equation as a bitmask and run elimination:

cpp
// Solve system of m XOR equations in n unknowns (n <= 63).
// eq[i] has bits 0..n-1 as coefficients, bit n as RHS.
// Returns rank, or -1 if inconsistent.
int gauss_xor_system(vector<long long>& eq, int n) {
    int m = (int)eq.size();
    int row = 0;
    for (int col = 0; col < n && row < m; col++) {
        int piv = -1;
        for (int r = row; r < m; r++) {
            if ((eq[r] >> col) & 1) { piv = r; break; }
        }
        if (piv == -1) continue;
        swap(eq[row], eq[piv]);
        for (int r = 0; r < m; r++) {
            if (r != row && ((eq[r] >> col) & 1))
                eq[r] ^= eq[row];
        }
        row++;
    }
    // Check consistency: any row with all coeff bits zero but RHS = 1
    for (int r = row; r < m; r++)
        if ((eq[r] >> n) & 1) return -1;
    return row;  // rank
}

5. Complexity

OperationTimeSpace
Gauss over reals / mod p (n equations, m unknowns)O(nmmin(n,m))O(nm)
XOR basis insert (single value, B bits)O(B)O(B)
XOR basis build (n values, B bits)O(nB)O(B)
XOR basis query (max/min/contains)O(B)O(1)
Gauss on bitmask equations (m eqs, n vars)O(mn)O(m)

Typical contest bounds: n,m500 for dense systems (runs in <0.5s), B=30 or 60 for XOR basis.

PARTIAL PIVOTING -- why the largest element matters (reals only):

  Without pivoting:                With partial pivoting:
  ┌                  ┐             ┌                  ┐
  │ 0.0001  1 │  1   │             │ 1      0.5 │ 0.5 │  <- swapped row 2 up
  │ 1      0.5│ 0.5  │             │ 0.0001  1  │  1  │
  └                  ┘             └                  ┘

  Multiply row 1 by 10000:        Multiply row 2 by ~0.0001:
  ┌                      ┐        ┌                      ┐
  │ 0.0001  1    │  1    │        │ 1    0.5       │ 0.5     │
  │ 0      -9999 │-9999  │        │ 0    0.99995   │ 0.99995 │
  └                      ┘        └                         ┘

  -> y ~= 1, x = (1 - 10000)/0.0001  -> y ~= 1, x ~= 0
  -> catastrophic cancellation!       -> stable result OK

  RULE: Always swap the row with max |a[i][col]| to pivot position.
        Over Z/p or GF(2), any nonzero works -- no stability issue.

6. Patterns

Pattern 1 -- Solve a Linear System

Problem shape: Given n equations in n unknowns, find the solution (or report no/infinite solutions).

Approach: Build augmented matrix, run gauss_real or gauss_mod.

cpp
// Example: solve 2x + y - z = 8, -3x - y + 2z = -11, -2x + y + 2z = -3
vector<vector<double>> a = {
    { 2,  1, -1,   8},
    {-3, -1,  2, -11},
    {-2,  1,  2,  -3}
};
vector<double> ans;
int sols = gauss_real(a, ans);
// sols == 1, ans = {2, 3, -1}

Pattern 2 -- Maximum XOR Subset

Problem shape: Given an array of n integers, find the maximum XOR of any nonempty subset.

Approach: Insert all values into an XOR basis, then call max_xor().

cpp
vector<long long> vals = {6, 5, 3, 2};
XorBasis xb;
for (long long v : vals) xb.insert(v);
cout << xb.max_xor() << "\n";  // 7
MAXIMUM XOR SUBSET -- greedy on the basis:

  Basis (after elimination):
    basis[3] = 1010  (bit 3)
    basis[2] = 0110  (bit 2)
    basis[0] = 0011  (bit 0)

  Greedy walk (start with result = 0):

    Bit 3: result = 0000
           result ^ basis[3] = 1010 > 0000?  YES -> result = 1010
    
    Bit 2: result = 1010
           result ^ basis[2] = 1100 > 1010?  YES -> result = 1100
    
    Bit 1: no basis element for bit 1 -> skip
    
    Bit 0: result = 1100
           result ^ basis[0] = 1111 > 1100?  YES -> result = 1111

  Maximum XOR = 1111 = 15

  ┌────────────────────────────────────────────┐
  │ The greedy works because each basis vector │
  │ controls a unique highest bit. Choosing to │
  │ XOR it in can only set that bit (never     │
  │ unset higher bits already set).            │
  └────────────────────────────────────────────┘

Pattern 3 -- Maximum XOR Path in a Graph

Problem shape: Weighted undirected graph, find the path from s to t maximizing XOR of edge weights.

Key idea: Any path from s to t differs from any other by a set of cycles. Find any path's XOR, collect all cycle XORs into a basis, then maximize using max_xor_with(path_xor).

cpp
// After DFS from s, dist[v] = XOR distance from s to v.
// For each back edge (u, v, w): cycle XOR = dist[u] ^ dist[v] ^ w.
XorBasis cycles;
// ... collect cycle XORs during DFS ...
long long ans = cycles.max_xor_with(dist[t]);

Pattern 4 -- Count Distinct XOR Values

Problem shape: Given n numbers, how many distinct values can be formed by XOR-ing subsets (including the empty subset = 0)?

Approach: Build basis, answer is 2rank.

COUNTING XOR-REACHABLE VALUES:

  Given set S = {5, 3, 6}  (101, 011, 110)
  
  XOR basis rank = 2 (only 2 independent vectors)
  Number of distinct achievable XOR values = 2² = 4

  Enumeration (verify):
    XOR of {}      = 000 = 0
    XOR of {5}     = 101 = 5
    XOR of {3}     = 011 = 3
    XOR of {5,3}   = 110 = 6
    XOR of {6}     = 110 = 6  <- same as {5,3}!
    XOR of {5,6}   = 011 = 3  <- same as {3}!
    XOR of {3,6}   = 101 = 5  <- same as {5}!
    XOR of {5,3,6} = 000 = 0  <- same as {}!

  Unique values: {0, 3, 5, 6} -> exactly 4 = 2²  OK

  ┌──────────────────────────────────────────────┐
  │ Formula: # distinct XOR values = 2^(rank)    │
  │ This counts the EMPTY subset (XOR=0) too.    │
  │ If the problem excludes empty: subtract 1    │
  │ (but 0 might still be achievable by a non-   │
  │ empty subset -- check if rank < |S|).         │
  └──────────────────────────────────────────────┘

Pattern 5 -- Matrix Rank

Problem shape: Given an n×m binary matrix, find its rank over GF(2).

Approach: Run XOR Gaussian elimination on the rows (each row is a bitmask). The rank equals the number of nonzero rows after elimination.

Pattern 6 -- Counting Solutions mod p

Problem shape: Given a system of linear equations over Z/p, count the number of solutions.

Approach: Run gauss_mod. If inconsistent, answer is 0. Otherwise the number of free variables is mrank, and the number of solutions is pmrank.


7. Cheat Sheet

+--------------------------------------------------------------+
|             GAUSSIAN ELIMINATION CHEAT SHEET                  |
+--------------------------------------------------------------+
| WHEN TO USE:                                                  |
|  - "system of equations"      --> gauss over R / mod p        |
|  - "max XOR of subset"        --> XOR basis + greedy          |
|  - "count distinct XORs"      --> 2^(basis rank)              |
|  - "XOR path in graph"        --> any path + cycle basis      |
|  - "matrix rank over GF(2)"   --> Gauss on bitmask rows       |
|  - "can we XOR to get x?"     --> basis.contains(x)           |
+--------------------------------------------------------------+
| KEY CODE:                                                     |
|  XorBasis xb;                                                 |
|  xb.insert(val);             // O(BITS) per insert            |
|  xb.max_xor();              // greedy top-down                |
|  xb.contains(x);            // reduce x through basis         |
|  xb.count_distinct();       // 1LL << rank                    |
+--------------------------------------------------------------+
| COMPLEXITIES:                                                 |
|  Dense system (n eqs, m vars): O(n * m * min(n,m))            |
|  XOR basis build (n vals)    : O(n * BITS)                    |
|  XOR basis query             : O(BITS)                        |
+--------------------------------------------------------------+
| FIELD CHOICE:                                                 |
|  Reals  : use EPS = 1e-9, partial pivoting, doubles           |
|  Mod p  : p must be prime, use Fermat inverse                 |
|  GF(2)  : XOR = add, AND = multiply, basis[i] or bitmask     |
+--------------------------------------------------------------+
| PITFALL: over reals, EPS too large --> wrong pivots           |
|          over reals, EPS too small --> div-by-near-zero        |
|          mod p, p not prime --> modinv fails silently          |
+--------------------------------------------------------------+

Gotchas & Debugging

  • Floating-point precision. Over reals, use EPS = 1e-9 and partial pivoting. If the system is ill-conditioned (pivot close to zero even after row swaps), the answer may be garbage. Switch to modular arithmetic if the problem has integer coefficients and a prime modulus.

  • Forgetting to normalize the pivot row (mod p). After selecting a pivot, divide the entire row by the pivot value. Skipping this means elimination produces wrong coefficients.

  • Wrong loop bounds in XOR basis. If values can be up to 1018, you need BITS = 60. Using BITS = 30 silently ignores high bits.

  • Empty basis edge case. max_xor() on an empty basis returns 0. If the problem says "nonempty subset", check that the basis rank is 1 and handle the case where all values are 0 separately.

  • Off-by-one in GF(2) equation bit. When encoding the RHS of an XOR equation, store it in bit n (the (n+1)-th bit). Confusing bit n with bit n1 corrupts the entire system.

  • Not checking consistency. After elimination, always scan the zero rows for nonzero RHS. Returning "rank = r" without this check can produce wrong "number of solutions" answers.

  • Modular inverse of zero. If a pivot candidate is 0(modp), skip it. Calling modinv(0, p) returns 0 via Fermat -- silently wrong.

  • XOR basis merge is O(B2). Merging two bases of size B costs O(B2) because each of the B insertions is O(B). For segment tree with XOR basis nodes, total is O(nB2logn) -- make sure B and n are small enough.

  • Greedy min XOR. min_xor() returns the smallest basis element, which is the minimum nonzero XOR. The minimum XOR including 0 is always 0 (empty subset). Read the problem statement carefully.

Mental Traps

  • "Gaussian elimination is for solving linear equations. This problem isn't about equations." Many problems that aren't explicitly about linear equations reduce to linear algebra. XOR basis problems are Gaussian elimination over GF(2). Counting vectors that can be expressed as XOR combinations of a given set = dimension of the linear span. Even some graph problems (cycle space of a graph) are linear algebra over GF(2).

  • "The solution to Ax = b is unique if the system is consistent." Only if the matrix has full rank. If the rank is less than n, the solution space is a coset with 2nr elements (over GF(2)) or a translation of a subspace (over reals). Many problems ask for the number of solutions, not a specific one.

  • "I can just reduce to row echelon form without tracking pivots." Pivot tracking is essential for: determining rank (number of pivot rows), identifying free variables (non-pivot columns), and reconstructing the solution. If you don't track which column each pivot came from, you can't extract meaningful information from the reduced form.

  • "I'll just use real-valued Gaussian elimination for everything." Different fields demand different arithmetic. Over R (doubles) you do floating-point division with partial pivoting. Over Z/p you must use modular inverse (a * pow(b, p-2, p) % p), not floating division -- one float division mod a prime and your answer is garbage. Over GF(2) you must use XOR and bitsets -- there is no "division," and using int arithmetic with % 2 everywhere is 64x slower than bitset operations and invites sign bugs with negative numbers.

  • "The XOR basis is something completely different from Gaussian elimination." It is literally the same algorithm, just specialized to GF(2). Recognizing this unification is powerful: every technique for Gaussian elimination (rank, free variables, solution counting, basis extraction) transfers directly to XOR basis problems, and vice versa.

THE UNIFICATION -- same algorithm, three fields:

  REALS (doubles)              Z/p (mod prime)             GF(2) (XOR basis)
  ─────────────────            ─────────────────           ─────────────────
  for each column:             for each column:            for each bit position:
    find pivot (max |val|)       find pivot (any !=0)         find vector with this
    swap rows                    swap rows                     bit set
    scale pivot row to 1         scale by mod inverse        (no scaling needed)
    subtract from others         subtract mod p              XOR with others
                                 from others

  Division: a / b              Division: a * b^(p-2) % p   "Division": XOR
  Danger: floating-point       Danger: forgetting           Danger: using % 2
          precision                    modular inverse              instead of ^

  Result: solution to Ax=b     Result: solution mod p       Result: XOR basis
          rank, determinant            rank, det mod p              rank = dim(span)

  ALL THREE produce row echelon form. ALL THREE give you rank.
  ALL THREE let you count solutions: inf / p^(n-r) / 2^(n-r).

Trap: "The rank of the matrix equals the number of non-zero rows after elimination."

Only if your implementation correctly avoids counting rows that were supposed to be pivots but got zeroed out. A subtle bug: if you process rows in order without properly swapping the pivot row to the current position, you may end up with non-zero rows that aren't actually independent pivots. Always count the number of pivot columns found, not just non-zero rows.

RANK COUNTING -- correct vs incorrect:

  After elimination:
  
  OK CORRECT (tracking pivot columns):
    [1  0  1  0]  <- pivot in col 0
    [0  1  1  0]  <- pivot in col 1
    [0  0  0  1]  <- pivot in col 3
    [0  0  0  0]
    Rank = 3 (three pivot columns: 0, 1, 3)
    Free variables: col 2 (no pivot)

  X INCORRECT (counting non-zero rows naively):
    [1  0  1  0]  <- looks like a pivot row
    [0  0  1  0]  <- non-zero, but col 1 had no pivot!
    [0  0  0  1]  <- pivot
    [0  0  0  0]
    Counting non-zero = 3, but is col 1 a pivot? NO.
    This matrix wasn't properly reduced.

  ┌─────────────────────────────────────────────────┐
  │ FIX: Track pivot_col[] during elimination.      │
  │ Rank = len(pivot_col).                          │
  │ Free variables = columns NOT in pivot_col.      │
  └─────────────────────────────────────────────────┘

Trap: "I can mix GF(2) operations with regular integer arithmetic."

In GF(2), 1 + 1 = 0. In integers, 1 + 1 = 2. If you're building an XOR basis and accidentally use + instead of ^, or compare with == 0 after integer subtraction instead of XOR, results are wrong but may accidentally pass small tests (since XOR and subtraction agree on single bits). Always use bitwise operations for GF(2) -- never +, -, or * on row values.

GF(2) vs INTEGER ARITHMETIC -- the silent bug:

  Operation on rows [1,0,1] and [1,1,0]:

  Integer add:    [1,0,1] + [1,1,0] = [2,1,1]  <- WRONG for GF(2)
  GF(2) XOR:      [1,0,1] ^ [1,1,0] = [0,1,1]  <- CORRECT

  As bitsets:
    101 + 110 = ??? (meaningless as bitwise add)
    101 ^ 110 = 011  OK

  In code:
    X row[i] = row[i] + row[j];     // integer add
    OK row[i] = row[i] ^ row[j];     // GF(2) elimination
    OK bs[i] ^= bs[j];               // bitset version

  ┌─────────────────────────────────────────────┐
  │ GF(2) has ONLY two elements: {0, 1}         │
  │ Addition = XOR, Multiplication = AND         │
  │ There is no "2" in GF(2).                   │
  └─────────────────────────────────────────────┘

Self-Test

Before moving on, you should be able to answer these without looking at the notes:

  1. Write the XOR basis insertion function from memory -- the loop from high bit to low, the XOR with existing basis vector or assignment to empty slot.

  2. Compute the maximum XOR of any subset of {3, 5, 6, 8} -- build the XOR basis and run the greedy algorithm.

  3. State how many distinct XOR values are achievable from a set with XOR basis of rank r -- and explain why.

  4. Write Gaussian elimination over GF(2) using bitsets -- rows as bitsets, XOR for row operations, pivot tracking.

  5. Explain how to detect an inconsistent system -- what does the matrix look like after elimination when there's no solution?

  6. Given a 3x3 system over the reals, perform elimination by hand -- including partial pivoting and back-substitution. What changes if you work mod a prime p?

  7. Explain the relationship between rank, free variables, and solution count -- for systems over R, Z/p, and GF(2).

  • [ ] Convert an XOR basis from REF to RREF and use it to find the k-th smallest XOR -- given basis = {101, 011, 001}, reduce to RREF and find the 5th smallest achievable XOR value.

  • [ ] Model a graph's cycle space as a vector space over GF(2) -- given a graph with V vertices, E edges, and c components, state the dimension of the cycle space and explain how non-tree edges form the fundamental cycles.

  • [ ] Identify when to use partial pivoting vs. any-nonzero pivoting -- explain why partial pivoting is critical over R (numerical stability) but unnecessary over Z/p and GF(2).

  • [ ] Debug an off-by-one rank error -- if you compute rank = r but the correct rank is r+1 (or r1), how does this affect the solution count? Trace through a concrete 3-variable GF(2) system where this matters.

  • [ ] Determine whether a given value is XOR-reachable from a set -- insert the set into an XOR basis, then check membership. Explain what happens at each step and why a zero remainder means "reachable."

  • [ ] Build the XOR basis for the set {12, 10, 6, 5, 3} (binary: 1100, 1010, 0110, 0101, 0011) -- show the state of basis[] after each insertion.

  • [ ] Compute the maximum XOR of any subset of {7, 3, 5, 9} by building the XOR basis and running the greedy from the highest bit.

  • [ ] Solve the system {x1x3=1x1x2=0x2x3=1 over GF(2) -- set up the augmented matrix, eliminate, and determine all solutions.

  • [ ] Given a graph with 5 edges, find the cycle space dimension (number of independent cycles = mn+c where c = connected components) and explain how this relates to the null space of the incidence matrix over GF(2).

  • [ ] Explain why partial pivoting is needed for floating-point Gaussian elimination but not for GF(2) -- give a numerical example where a small pivot causes catastrophic cancellation over R.


Practice Problems

#ProblemSourceDifficultyKey Idea
1XOR on SegmentCF 1847C1500XOR prefix + basis intuition warmup
2Xor-MSTCF 888G2100XOR basis + divide-and-conquer on trie
3Linear Algebra TestCF 832E2100Gauss over GF(2), maximize subset XOR
4Maximum Xor SubarraySPOJ XMAX1800Classic XOR basis -- max XOR subset
5Xor SumCF Gym1900XOR basis construction + greedy query
6Yet Another MinimizationCF 1516D1800XOR + basis rank for counting
7XOR EquationCF 1710B1800System of XOR constraints
8Toy TrainCF 1556E2000XOR path in graph + cycle basis
9Basis ChangeCF 895C1900Count subsets with XOR divisible by 3 (rank + counting)
10GukiZ and MatricesCF 551E2100Gauss over Z/p, matrix rank

Further Reading


The Aha Moment

Gaussian elimination is just systematic substitution -- the same thing you did in middle school algebra, but organized as row operations so a computer can do it. The magic is that it works identically over any field: real numbers, integers modulo a prime, or GF(2) (where "addition" is XOR and "multiplication" is AND). This universality means one algorithm unlocks three families of problems.

The real "aha" for CP: over GF(2), Gaussian elimination becomes the XOR basis -- row-reduce a set of binary numbers to find the maximum XOR, count XOR-reachable values, or solve systems of XOR equations. This is the same algorithm, just with ^ instead of + and & instead of *.


Pattern Fingerprint

Constraint / SignalTechnique
System of n linear equations, n500Standard Gaussian elimination O(n3)
"How many solutions?" / "Does a solution exist?"Row reduce and check rank vs. augmented rank
"Maximum XOR of a subset"Gaussian elimination over GF(2) (XOR basis)
"Number of XOR-distinct values reachable"2rank where rank = number of basis vectors
Binary strings/masks with XOR constraintsGF(2) Gaussian elimination on bitmasks
"Solve modular equations modp" (p prime)Gaussian elimination with modular inverse
"Rank of a matrix"Count non-zero rows after elimination
"Determinant of a matrix"Product of diagonal after elimination (track sign flips)

Rating Progression

CF RatingWhat You Should Be Able To Do
1200Solve 2×2 systems by hand; understand what "no solution" and "infinite solutions" mean
1500Implement basic Gaussian elimination over reals; compute matrix rank
1800Gaussian elimination over GF(2) for XOR problems; find maximum XOR subset; count solutions
2100Solve underdetermined systems (free variables); combine Gaussian elimination with DP or probability; handle modular arithmetic fields

Before You Code Checklist

  • [ ] Choose the field -- reals (use double with epsilon), integers mod prime (use modular inverse), or GF(2) (use bitwise XOR)
  • [ ] Partial pivoting for reals -- always swap in the row with the largest absolute value in the pivot column to minimize floating-point error
  • [ ] Check for no-solution -- after elimination, a row like [000|c] with c0 means inconsistency
  • [ ] Count free variables -- free=variablesrank; number of solutions = pfree over Fp, or 2free over GF(2)
  • [ ] For GF(2): use bitmask representation -- each equation is a long long, XOR replaces addition, no need for division

What Would You Do If...?

Scenario 1: "The system has more variables than equations." The system is underdetermined. After Gaussian elimination, some variables are "free" -- you can set them to anything and back-substitute to get the rest. The number of solutions is p(varsrank) if working over Fp.

Scenario 2: "The matrix entries are floating-point and you're getting wrong answers." Floating-point Gaussian elimination is numerically unstable without partial pivoting. Always swap in the row with the largest |a[i][col]| as the pivot. Also use an epsilon (e.g., 109) for comparisons instead of == 0. If precision is critical, consider working with fractions or modular arithmetic.

Scenario 3: "You need to solve the same system with many different right-hand sides." Row-reduce the coefficient matrix once (LU decomposition), then for each right-hand side, just do forward and back substitution -- O(n2) per RHS instead of O(n3). Alternatively, compute the inverse matrix once (O(n3)) and multiply by each RHS.


The Mistake That Teaches

Story: You're solving a problem over GF(2): given n bitmasks, find the maximum XOR of any subset. You implement Gaussian elimination on the bitmasks. It works on small tests. You submit -- Wrong Answer.

After debugging, you realize: when you find a pivot row for bit b, you XOR it into rows below -- but you forget to also XOR it into rows above. Standard Gaussian elimination for solving systems only eliminates downward, but for the XOR basis, you want reduced row echelon form (eliminate in both directions) so each basis vector has a unique leading bit.

Wait, actually that's not the issue -- for maximum XOR, you greedily pick the highest bit, which works even without full reduction. The real bug: you're iterating bits from low to high instead of high to low. The greedy construction of the maximum XOR needs the highest bit first.

The fix: Iterate from bit 59 (or 29) downward, and for the "maximum XOR" query, greedily XOR in each basis vector if it increases the result.

Lesson: The direction of iteration matters -- for maximum XOR, process high bits first. And always double-check whether you need REF or RREF.


Debug This

Bug 1:

Find the bug in this Gaussian elimination (reals)
cpp
// Solve Ax = b, n equations, n unknowns
for (int col = 0; col < n; col++) {
    // Find pivot
    int pivot = -1;
    for (int row = col; row < n; row++) {
        if (a[row][col] != 0) { pivot = row; break; }
    }
    if (pivot == -1) continue;  // no pivot in this column
    swap(a[col], a[pivot]);
    // Eliminate below
    for (int row = col + 1; row < n; row++) {
        double factor = a[row][col] / a[col][col];
        for (int j = col; j <= n; j++)
            a[row][j] -= factor * a[col][j];
    }
}

Bug: No partial pivoting -- it picks the first nonzero element, not the largest. For floating-point arithmetic, this leads to numerical instability when the pivot is very small (e.g., 1015), causing division by a tiny number and amplified errors.

Fix: Replace the pivot search with:

cpp
int pivot = col;
for (int row = col + 1; row < n; row++)
    if (fabs(a[row][col]) > fabs(a[pivot][col])) pivot = row;
if (fabs(a[pivot][col]) < 1e-9) continue;  // effectively zero

Bug 2:

Find the bug in this GF(2) Gaussian elimination
cpp
// XOR basis: given array of bitmasks, build basis
vector<long long> basis;
for (long long val : nums) {
    for (long long b : basis) {
        val = min(val, val ^ b);  // reduce val by basis
    }
    if (val > 0) basis.push_back(val);
}

Bug: min(val, val ^ b) is wrong. The intention is to eliminate the highest set bit of val, but min doesn't do that -- it just picks the numerically smaller value, which doesn't guarantee bit elimination. The correct approach is to check if XORing with the basis element reduces the highest bit.

Fix: Use the standard approach -- iterate from the highest bit:

cpp
for (long long b : basis) {
    val = min(val, val ^ b);
}

Actually, min(val, val ^ b) does work as a valid insertion strategy (it's a known trick that always reduces the value), but the basis must be sorted or the iteration order must guarantee convergence. The safer standard approach:

cpp
long long basis[60] = {};
for (long long val : nums) {
    for (int i = 59; i >= 0; i--) {
        if (!(val >> i & 1)) continue;
        if (!basis[i]) { basis[i] = val; break; }
        val ^= basis[i];
    }
}

Bug 3:

Find the bug in this solution counter
cpp
// Count solutions to Ax = b over GF(2)
int rank = 0;
for (int col = 0; col < m; col++) {
    int pivot = -1;
    for (int row = rank; row < n; row++)
        if (a[row][col]) { pivot = row; break; }
    if (pivot == -1) continue;
    swap(a[rank], a[pivot]);
    for (int row = 0; row < n; row++) {
        if (row != rank && a[row][col])
            for (int j = 0; j <= m; j++) a[row][j] ^= a[rank][j];
    }
    rank++;
}
// Check consistency
for (int row = rank; row < n; row++)
    if (a[row][m]) { cout << 0; return; }  // no solution

cout << (1 << (m - rank));  // 2^(free variables)

Bug: 1 << (m - rank) overflows if m - rank > 30 (since 1 is an int, which is 32 bits). For m=60 variables and rank 0, this is undefined behavior.

Fix: Use 1LL << (m - rank) to get a long long result.


Common Mistakes

  • No partial pivoting over reals. Always swap in the row with the largest absolute value in the pivot column to avoid numerical instability.
  • EPS too large or too small. Too large drops valid pivots; too small keeps noise. Use 1e-9 as a starting point.
  • Overflow in 1 << (m - rank). Use 1LL << (m - rank) when counting solutions over GF(2).
  • Forgetting to check consistency. After elimination, check if any row has all-zero coefficients but a non-zero constant -- that means no solution.

Historical Origin

The method is named after Carl Friedrich Gauss, who used it in the early 1800s for astronomical calculations (specifically, fitting orbital parameters for the asteroid Ceres). However, the technique was known to Chinese mathematicians as early as 179 AD in the Jiuzhang Suanshu ("Nine Chapters on the Mathematical Art"), making it one of the oldest algorithms still in everyday use.


Mnemonic Anchor

"Row reduce over any field -- reals for systems, GF(2) for XOR, rank for counting."


ConceptLink
Math Fundamentals01-math-fundamentals
Bit Manipulation07-bit-manipulation
Linear Basis (XOR)11-linear-basis-xor
Matrix Exponentiation09-matrix-exponentiation
PracticeLadder - Mixed

Built for the climb from Codeforces 1100 to 2100.