Gaussian Elimination for Solving Linear Systems
Gaussian elimination solves systems of linear equations in $O(n^3)$ time.
Consider a system of $n$ linear equations in $n$ unknowns:
$$ a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n = b_1 $$ $$ \vdots $$ $$ a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n = b_n $$
Gauss-Jordan Elimination
Algorithm:
- Select an unpicked unknown as pivot, and an equation containing this pivot.
- Scale the pivot equation so the coefficient of the pivot becomes 1.
- Eliminate this unknown from all other equations using subtraction.
- Repeat until each equation contains a distinct pivot.
After completion, the system becomes: $$ c_1 x_1 = d_1 $$ $$ \vdots $$ $$ c_n x_n = d_n $$
The solution is $x_i = d_i / c_i$.
- If $c_i = 0$ but $d_i \neq 0$, there is no solution.
- If $c_i = 0$ and $d_i = 0$, there are infinitely many solutions.
Code Examples
Basic Gaussian Elimination (P3389)
This implementation transforms the matrix into reduced row echelon form. If no unique solution exists, output "No Solution".
#include <bits/stdc++.h>
#define N 105
#define db double
using namespace std;
int n;
db a[N][N];
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++)
for(int j=1;j<=n+1;j++)
scanf("%lf",&a[i][j]);
for(int i=1;i<=n;i++){
int mx=i;
for(int j=i+1;j<=n;j++)
if(fabs(a[j][i])>fabs(a[mx][i])) mx=j;
for(int j=1;j<=n+1;j++)
swap(a[i][j],a[mx][j]);
if(!a[i][i]){
printf("No Solution\n");
return 0;
}
for(int j=1;j<=n;j++){
if(j==i) continue;
db factor = a[j][i] / a[i][i];
for(int k=i+1;k<=n+1;k++)
a[j][k] -= a[i][k] * factor;
}
}
for(int i=1;i<=n;i++)
printf("%.2lf\n",a[i][n+1]/a[i][i]);
return 0;
}
System with Possible Infinite/No Solutions (P2455)
This version handles non-unique cases: output -1 for no solution, 0 for infinite solutions.
#include <bits/stdc++.h>
#define N 55
#define db double
using namespace std;
int n, cur;
db a[N][N];
int main(){
scanf("%d",&n);
for(int i=0;i<n;i++)
for(int j=0;j<n+1;j++)
scanf("%lf",&a[i][j]);
for(int i=0;i<n;i++){
int mx = cur;
for(int j=cur+1;j<n;j++)
if(fabs(a[j][i]) > fabs(a[mx][i])) mx=j;
if(!a[mx][i]) continue;
for(int j=0;j<n+1;j++)
swap(a[cur][j], a[mx][j]);
for(int j=0;j<n;j++){
if(j==cur) continue;
db factor = a[j][i] / a[cur][i];
for(int k=i+1;k<n+1;k++)
a[j][k] -= a[cur][k] * factor;
}
cur++;
}
if(cur < n){
for(int i=cur;i<n;i++)
if(a[i][n]){ printf("-1\n"); return 0; }
printf("0\n");
return 0;
}
for(int i=0;i<n;i++)
printf("x%d=%.2lf\n",i+1,a[i][n]/a[i][i]);
return 0;
}
XOR System – Lights Out (P2962)
Given an undirected graph, each node is initially 0. Flipping a node toggles itself and its neighbors. Find minimal flips to make all nodes 1. This reduces to solving a linear system over GF(2). Free variables are enumerated with DFS and pruning.
#include <bits/stdc++.h>
#define N 40
using namespace std;
int n, m;
int a[N][N], ans = N;
int f[N];
void gauss(){
for(int i=1;i<=n;i++){
int mx=i;
for(int j=i+1;j<=n;j++)
if(a[j][i]){ mx=j; break; }
for(int j=1;j<=n+1;j++)
swap(a[mx][j], a[i][j]);
for(int j=i+1;j<=n;j++)
if(a[j][i]){
for(int k=1;k<=n+1;k++)
a[j][k] ^= a[i][k];
}
}
}
void dfs(int x, int tot){
if(tot > ans) return;
if(!x){
ans = min(ans, tot);
return;
}
if(a[x][x]){
f[x] = a[x][n+1];
for(int i=n;i>x;i--)
f[x] ^= f[i] & a[x][i];
dfs(x-1, tot+f[x]);
} else {
f[x]=0; dfs(x-1, tot);
f[x]=1; dfs(x-1, tot+1);
}
}
int main(){
scanf("%d%d",&n,&m);
for(int u,v;m;m--){
scanf("%d%d",&u,&v);
a[u][v]=a[v][u]=1;
}
for(int i=1;i<=n;i++)
a[i][i]=a[i][n+1]=1;
gauss();
dfs(n,0);
printf("%d\n",ans);
return 0;
}
Determining Unique Solution (P2447)
Equations are added sequentially; bitset helps to accelerate. If a unique solution exists, output the step when its determined.
#include <bits/stdc++.h>
#define N 1010
#define M 2010
#define bs bitset<N>
using namespace std;
int n, m;
bs a[M];
int ans;
int getBit(){
char ch;
while(!isdigit(ch=getchar()));
return ch-'0';
}
int main(){
scanf("%d%d",&n,&m);
for(int i=1;i<=m;i++){
for(int j=1;j<=n+1;j++)
a[i][j]=getBit();
}
for(int i=1;i<=n;i++){
int cur=i;
while(cur<=m && !a[cur][i]) cur++;
if(cur==m+1){ printf("Cannot Determine\n"); return 0; }
ans = max(ans, cur);
swap(a[cur], a[i]);
for(int j=1;j<=m;j++){
if(j==i) continue;
if(a[j][i]) a[j] ^= a[i];
}
}
printf("%d\n",ans);
for(int i=1;i<=n;i++)
printf(a[i][n+1]?"?y7M#\n":"Earth\n");
return 0;
}
Parsing Problem – JSOI 2012 (P6126)
Partition vertices into two sets such that for each vertex, the number of its outgoing neighbors in the same set is even. Model each vertex's set membership as variable $x_i \in {0,1}$. For vertex $i$, let $S_i$ be its adjacency list. The condition is:
$$ \bigoplus_{j \in S_i} (x_j \oplus x_i \oplus 1) = 0 $$
which simplifies to:
$$ \left(\bigoplus_{j \in S_i} x_j\right) \oplus (|S_i| \bmod 2 \cdot x_i) = |S_i| \bmod 2 $$
Use bitset to solve the XOR system.
#include <bits/stdc++.h>
#define N 2010
using namespace std;
int n;
bitset<N> a[N];
int ans;
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++){
int m; scanf("%d",&m);
a[i][i]=a[i][n+1]=m&1;
while(m--){
int v; scanf("%d",&v);
a[i][v]=1;
}
}
int cnt=0;
for(int i=1;i<=n;i++){
int mx=i;
for(int j=mx;j<=n;j++)
if(a[j][i]){ mx=j; break; }
swap(a[mx], a[i]);
if(!a[i][i]) continue;
cnt++;
for(int j=1;j<=n;j++)
if(j!=i && a[j][i]) a[j] ^= a[i];
}
if(cnt < n){
for(int i=cnt+1;i<=n;i++)
if(!a[i][i] && a[i][n+1]){ printf("Impossible\n"); return 0; }
}
for(int i=1;i<=n;i++) ans += a[i][n+1];
printf("%d\n",ans);
for(int i=1;i<=n;i++) if(a[i][n+1]) printf("%d ",i);
printf("\n");
return 0;
}
Probability DP with Gaussian Elimination – BJOI 2018 (P4457)
This is a hard problem combining DP and Gaussian elimination. The hero's health is $p$, max health $n$. There are $m$ infinite-health minions. Each round: randomly heal one non-full friendly unit by 1, then randomly deal 1 damage $k$ times. Find expected number of rounds until hero dies.
Define $P(i)$ = probability of taking $i$ damage in one round:
$$ P(i) = \frac{\binom{k}{i} m^{k-i}}{(m+1)^k} $$
Let $f(x)$ be expected rounds starting from health $x$. Recurrence:
$$ f(x) = \frac{1}{m+1}\sum_{i=0}^{x+1} f(i) P(x+1-i) + \frac{m}{m+1}\sum_{i=0}^{x} f(i) P(x-i) + 1 $$
$$ f(n) = \sum_{i=0}^{n} P(i) f(n-i) + 1 $$
Boundary: $f(0)=0$.
Rearranged in to linear equations, each involving only $f(i)$ and $f(i+1)$, enabling $O(n^2)$ elimination.
#include <bits/stdc++.h>
#define N 1505
#define Mod 1000000007
using namespace std;
int qpow(int a, int b, int p){
int res=1; a%=p;
while(b){
if(b&1) res = 1ll*res*a%p;
a = 1ll*a*a%p; b>>=1;
}
return res;
}
int inv(int x){ return qpow(x, Mod-2, Mod); }
int T, n, p, m, k;
int P[N], a[N][N];
int main(){
scanf("%d",&T);
while(T--){
scanf("%d%d%d%d",&n,&p,&m,&k);
if(k==0){ printf("-1\n"); continue; }
if(m==0){
if(k==1) printf("-1\n");
else{
if(p!=n) printf("%d\n", (p-1)/(k-1)+1);
else{
p-=k;
if(p>0) printf("%d\n", (p-1)/(k-1)+2);
else printf("1\n");
}
}
continue;
}
memset(P,0,sizeof(P));
memset(a,0,sizeof(a));
int im = inv(m), im1 = inv(m+1);
P[0] = qpow(1ll*m*im1%Mod, k, Mod);
for(int i=1;i<=min(n,k);i++)
P[i] = 1ll*P[i-1]*(k-i+1)%Mod*inv(i)%Mod*im%Mod;
for(int i=1;i<n;i++){
for(int j=1;j<=i;j++)
a[i][j] = (1ll*P[i+1-j] + 1ll*m*P[i-j]) % Mod * im1 % Mod;
a[i][i+1] = 1ll*P[0]*im1%Mod;
a[i][i] = (a[i][i]-1+Mod)%Mod;
a[i][n+1] = Mod-1;
}
for(int i=1;i<=n;i++) a[n][i] = P[n-i];
a[n][n] = (a[n][n]-1+Mod)%Mod;
a[n][n+1] = Mod-1;
for(int i=1;i<=n;i++){
if(!a[i][i]) continue;
int inv_pivot = inv(a[i][i]);
for(int j=1;j<=n;j++){
if(i==j || !a[j][i]) continue;
int factor = 1ll*a[j][i]*inv_pivot%Mod;
if(i<n) a[j][i+1] = (a[j][i+1] - 1ll*a[i][i+1]*factor % Mod + Mod) % Mod;
a[j][n+1] = (a[j][n+1] - 1ll*a[i][n+1]*factor % Mod + Mod) % Mod;
}
}
printf("%d\n", 1ll*a[p][n+1] * inv(a[p][p]) % Mod);
}
return 0;
}
