Gauss–Jordan Elimination

Gauss–Jordan Elimination in C++

This is my first attempt to solve linear systems with unique solutions in O(N³).

Algorithm:
1. Read augmented matrix for the system.
2. Perform Gaussian Elimination to get the matrix in row echelon form.
3. Transform the matrix to reduced row echelon form.
4. Write results.

Preview:

[click for a clear view]

My C++ solution (change MAX to desired size, it is the maximum limit of how many equations it can handle, change the macro DEBUG to "if(1)" if you want to see the intermediate steps):

#include <iostream>
#include <cstdlib>
#include <cassert>
#include <algorithm>
using namespace std;

#define DEBUG if(0)
#define MAX 5

struct MAT {
    int R[MAX+1];
} M[MAX];

int N;

int gcd(int a, int b) {
    while(b) b ^= a ^= b ^= a %= b;
    return a;
}

bool comp(MAT A, MAT B) {
    for(int k=0; k<N; k++) {
        if(A.R[k] > B.R[k]) return true;
        if(A.R[k] < B.R[k]) return false;
    }
    return true;
}

void print(void) {
    int i, j;
    cout << "..." << endl;
    for(i=0; i<N; i++) {
        for(j=0; j<N; j++) cout << M[i].R[j] << '\t';
        cout << M[i].R[j] << endl;
    }
    cout << "..." << endl;
}

void modify(MAT &A, MAT B, int a, int b) {
    for(int r=0; r<=N; r++) A.R[r] = a*B.R[r] - b*A.R[r];
}

void eliminate(void) {
    int i, g, k;
    for(i=0; i<N-1; i++) {
        if(!M[i].R[i]) {
            sort(&M[i], M+N, comp);
            assert(M[i].R[i]);
        }
        for(k=i+1; k<N; k++) {
            g = gcd(abs(M[i].R[i]), abs(M[k].R[i]));
            modify(M[k], M[i], M[k].R[i]/g, M[i].R[i]/g);
        }
        DEBUG print();
    }
}

void reduce(void) {
    int i, g, k;
    for(i=N-1; i; i--) {
        for(k=i-1; k>=0; k--) {
            g = gcd(abs(M[i].R[i]), abs(M[k].R[i]));
            modify(M[k], M[i], M[k].R[i]/g, M[i].R[i]/g);
        }
        DEBUG print();
    }
}

void printsol(void) {
    double h, l;
    cout << "Solve for " << N << " variables:" << endl;
    for(int i=0; i<N; i++) {
        assert(M[i].R[i]);
        h = M[i].R[i];
        l = M[i].R[N];
        cout << 'x' << i+1 << " = " << l/h << endl;
    }
    cout << "..." << endl;
}

int main(void) {
    int i, j, g;
    while(cin >> N) {
        if(!N) break;
        memset(M, 0, sizeof(M));
        for(i=0; i<N; i++) {
            for(j=0; j<N; j++) cin >> M[i].R[j];
            cin >> M[i].R[j];
            for(j=g=0; j<=N; j++) g = gcd(abs(M[i].R[j]), g);
            for(j=0; j<=N; j++) M[i].R[j] /= g;
        }
        eliminate();
        reduce();
        printsol();
    }
    return 0;
}

However, I'm trying to do it in a more efficient way. I'll post that when I'm done.