Gaussian Elimination

Gaussian elimination transforms a linear system into an equivalent simpler one by row operations. Pivot variables can be solved from the transformed equations; missing pivots identify free variables, and a zero coefficient row with nonzero right side reveals inconsistency. We implement scaled partial-pivot Gauss–Jordan elimination over floating-point real numbers, explaining both the exact algebra and the limits of tolerance-based numerical classification.

Prerequisites: Matrix Multiplication concepts; linear equations; floating-point precision.

Build the idea

1. 1. Build the augmented matrix

   For equations Ax=b, append b as a final column. Swapping rows, multiplying a row by a nonzero scalar, or adding a multiple of one row to another preserves the real solution set. A pivot is a nonzero coefficient selected to eliminate a variable from other rows. Rectangular systems are valid: the number of equations need not equal the number of variables.

   $$Ax=b$$

2. 2. Select and normalize pivots

   The code first scales each nonzero coefficient row by its largest absolute coefficient, then chooses the largest absolute available entry in the current column. This partial pivoting avoids unnecessarily dividing by a tiny entry when a better one exists. Swap into the next pivot row, divide by the pivot, and subtract suitable multiples from every other row. where[col] records which row solves that variable.

3. 3. Classify by rank and consistency

   In exact arithmetic, a row 0x₁+…+0xₖ=c with c≠0 makes the system impossible. Otherwise rank equal to the variable count gives a unique solution. If rank is smaller, at least one variable is free and there are infinitely many real solutions. Setting free variables to zero produces one particular solution; it does not enumerate the whole affine family.

4. 4. Work all three outcomes

   Equations 2x+y=5 and x−y=1 give unique x=2,y=1. Equations x+y=2 and 2x+2y=4 give infinitely many solutions, such as x=2,y=0 when y is chosen free. Replacing the second right side by 5 makes the rows inconsistent. The program’s sample prints UNIQUE and approximate coordinates 2.000000,1.000000.

5. 5. Separate exact theory from finite precision

   The implementation treats scaled magnitudes≤10^−12 as zero and checks a computed solution’s relative residual against original equations. A failure prints NUMERICAL FAILURE. Neither a small residual nor partial pivoting guarantees small solution error in an ill-conditioned system; rank/NO decisions near the tolerance are numerical judgments, not exact certificates. Over a prime modular field, use nonzero pivots and modular inverses instead of floating division, with exact zero tests. A composite modulus has noninvertible nonzero elements, so this field algorithm cannot be copied unchanged.

This is a complete program. Compile with C++17 and enter the input described below.

#include<iostream>
#include<vector>
#include<cmath>
#include<algorithm>
#include<iomanip>
int main(){
    int rows,variables;
    std::cin>>rows>>variables;
    using Real=long double;
    const Real EPS=1e-12L;
    std::vector<std::vector<Real>>a(rows,std::vector<Real>(variables+1));
    for(auto&row:a){
        for(auto&x:row){
            std::cin>>x;
        }
    }
    auto original=a;
    // Scale each nonzero coefficient row before partial pivoting.
    for(auto&row:a){
        Real scale=0;
        for(int j=0;j<variables;++j){
            scale=std::max(scale,std::abs(row[j]));
        }
        if(scale>0){
            for(auto&x:row){
                x/=scale;
            }
        }
    }
    std::vector<int>where(variables,-1);
    int rank=0;
    for(int col=0;col<variables&&rank<rows;++col){
        int pivot=rank;
        for(int i=rank+1;i<rows;++i){
            if(std::abs(a[i][col])>std::abs(a[pivot][col])){
                pivot=i;
            }
        }
        if(std::abs(a[pivot][col])<=EPS){
            continue;
        }
        std::swap(a[pivot],a[rank]);
        where[col]=rank;
        Real value=a[rank][col];
        for(int j=col;j<=variables;++j){
            a[rank][j]/=value;
        }
        for(int i=0;i<rows;++i){
            if(i!=rank){
                Real factor=a[i][col];
                for(int j=col;j<=variables;++j){
                    a[i][j]-=factor*a[rank][j];
                }
            }
        }
        ++rank;
    }
    for(int i=0;i<rows;++i){
        bool zero=true;
        for(int j=0;j<variables;++j){
            if(std::abs(a[i][j])>EPS){
                zero=false;
            }
        }
        if(zero&&std::abs(a[i][variables])>EPS){
            std::cout<<"NO\n";
            return 0;
        }
    }
    std::vector<Real>answer(variables,0); // Free variables chosen as zero.
    for(int j=0;j<variables;++j){
        if(where[j]!=-1){
            answer[j]=a[where[j]][variables];
        }
    }
    for(int i=0;i<rows;++i){
        Real value=0,scale=1+std::abs(original[i][variables]);
        for(int j=0;j<variables;++j){
            value+=original[i][j]*answer[j];
            scale+=std::abs(original[i][j]*answer[j]);
        }
        if(std::abs(value-original[i][variables])>1e-8L*scale){
            std::cout<<"NUMERICAL FAILURE\n";
            return 0;
        }
    }
    std::cout<<(rank==variables?"UNIQUE":"INFINITE")<<'\n';
    std::cout<<std::fixed<<std::setprecision(6);
    for(int j=0;j<variables;++j){
        std::cout<<(std::abs(answer[j])<5e-7L?0.0L:answer[j])<<(j+1==variables?'\n':' ');
    }
}

Complexity analysis

O(r·c·min(r,c)) time for r equations and c variables, O(r(c+1)) memory including an original copy for residual checks. Square n×n elimination is O(n³). Rows and variables are 1..30 with finite moderate-scale coefficients/right sides (|value|≤10^6); floating classification remains tolerance-dependent, especially for nearly singular inputs.

Example input and output

Read rows variables (1..30), then rows augmented lines containing variables real coefficients and one right-hand side, all finite with magnitude≤10^6. Print numerical NO, UNIQUE or INFINITE; successful cases also print one approximate solution with free variables zero. NUMERICAL FAILURE reports a failed residual check. Near-singular cases need caution as explained.

Input

2 2
2 1 5
1 -1 1

Output

UNIQUE
2.000000 1.000000

Adding the equations gives 3x=6, so x=2 and y=1; the residuals are zero.

Practical applications

Use for small real linear systems, recurrence coefficient inference and linear-algebra subproblems with a stated tolerance. For exact integer/rational classification use exact rational techniques with overflow/size management. For finite fields use a separately implemented modular elimination under a prime; the scalar arithmetic changes even though row logic is similar.

Check your understanding

In exact arithmetic, what causes no solution?

Put it into practice

Test a unique, dependent-consistent and contradictory system, then a nearly singular one.