/* 
 * File:   main.cpp
 * Author: Bangonkali
 *
 * Created on March 23, 2012, 8:14 PM
 */
#include <iostream>
#include <cmath>
#include <cstdlib>
 
 
// example taken from wikipedia 
// http://en.wikipedia.org/wiki/Gaussian_elimination#Example
 
using namespace std;
 
//void printMatrix(
//        int N,
//        float A[20][21]
//)
//{
//    for (int i = 0; i < N; i++)
//    {
//        for (int j = 0; j <= N; j++)
//        {
//            cout << A[i][j] << '\t';
//        }
//        
//        cout << endl;
//    }
//    
//    cout << endl;
//}
 
int gauss(
        int N, // number of unknowns
        float A [20] [21], // coefficients and constants
        float r[20]
)
{
//    printMatrix(N, A);
    float multiplier, divider;
    bool fin = false;
 
    // forward substitution
    for (int m=0; m<=N; m++)
    {
        for (int i=m+1; i<=N; i++)
        {
            multiplier = A[i][m]; // current row - root
            divider = A[m][m]; // first row - root
 
            if (divider == 0)
            {
                return 1;
            }
 
            for (int j=m; j<=N; j++)
            {
                if (i == N) 
                {
                    break;
                    fin = true;
                }
 
                A[i][j] = (A[m][j] * multiplier / divider) - A[i][j];
//                cout << "A[" << i << "][" << j << "]  = (A[" << m << "][" << j << "] * " << multiplier << " / " << divider << ") - A[" << i << "][" << j << "]" << endl;
//                cout << "Current cell " << i << j << endl;
//                cout << "Current set " << m << m << endl;
//                printMatrix(N, A);
            }
 
            for (int j=m; j<=N; j++)
            {
                A[i-1][j] = A[i-1][j] / divider;
            }
 
            if (fin)
                break;
        }
    }
 
//    printMatrix(N, A);
 
    // back substitution
    float s = 0;
    r[N-1] = A[N-1][N];
    int y = 0;
    for (int i = N-2; i >= 0; i--)
    {
        s = 0;
        y++;
        for (int x = 0; x < N; x++)
        {
            s = s + (A[i][N-1-x] * r[N-(x+1)]);
//            cout << "A[" << i << "][" << N-1-x << "] = " << A[i][N-1-x] << " ";
//            cout << " r[" << N-(x+1) << "] = " << r[N-(x+1)] << " " << endl;
//            cout << "s = " << s << endl;
 
            if (y == x+1) break;
        }
        r[i] = A[i][N] - s;
//        cout << "r[" << i << "] = A[" << i << "][" << N << "] - " << s << " = " << r[i] << endl;
    }
 
//    for (int i = 0; i < N; i++)
//    {
//        cout << r[i] << endl;
//    }
 
 
}
/*
 * 
 */
int main(int argc, char** argv) {
    float A[20][21];
    float X[20];
    bool err;
    A[0][0] = 2; 
    A[0][1] = 1; 
    A[0][2] = -1; 
    A[0][3] = 8;
 
    A[1][0] = -3; 
    A[1][1] = -1; 
    A[1][2] = 2; 
    A[1][3] = -11; 
 
    A[2][0] = -2; 
    A[2][1] = 1; 
    A[2][2] = 2; 
    A[2][3] = -3; 
 
    gauss(3, A, X);
    for (int i=0; i<3; i++) cout << X[i] << "  ";
 
    return 0;
}