// Arup Guha
// 3/26/2012

// First Attempt at a linear equation solver.
// This is set to work for up to 10 equations in 10 variables, where
// the equations are all "nice" and have a unique solution.

#include <stdio.h>

#define SIZE 10
#define DEBUG 0

void makeLowerDiag(double aug_matrix[][SIZE+1], int n);
void print(double aug_matrix[][SIZE+1], int n);
void finishSolvingmakeLowerDiag(double aug_matrix[][SIZE+1], int n);

int main() {

    double aug_matrix[SIZE][SIZE+1];

    // Open the input file.
    FILE* ifp = fopen("equations.txt", "r");

    // Read in the number of sets of equations to solve.
    int numprobs;
    fscanf(ifp, "%d", &numprobs);

    // Go through each set.
    int loop;
    for (loop=1; loop<=numprobs; loop++) {

        int n;
        fscanf(ifp, "%d", &n);

        // Read in augmented matrix.
        int i,j;
        for (i=0; i<n; i++) {
            for (j=0; j<n+1; j++) {
                fscanf(ifp,"%lf", &aug_matrix[i][j]);
            }
        }

        // First put in lower diagonal form.
        makeLowerDiag(aug_matrix, n);

        if (DEBUG) {
            print(aug_matrix, n);
            printf("\n\n");
        }

        // Finish solving
        finishSolvingmakeLowerDiag(aug_matrix,n);

        if (DEBUG) {
            print(aug_matrix, n);
            printf("\n\n");
        }

        // Output results.
        printf("Set of Equations #%d:\n",loop);
        for (i=0; i<n; i++)
            printf("Var %d = %lf.\n", i, aug_matrix[i][n]);
        printf("\n\n");

    }

    return 0;
}

// Prints out aug_matrix.
void print(double aug_matrix[][SIZE+1], int n) {

    int i,j;
    for (i=0; i<n; i++) {

        for (j=0; j<n+1; j++)
            printf("%6.2lf", aug_matrix[i][j]);
        printf("\n");
    }
}

// Changes aug_matrix to lower diagonal form, changing
// the equations to an equivalent set of equations.
void makeLowerDiag(double aug_matrix[][SIZE+1], int n) {

    int i,j,k;

    // Using row (i-1) as a pivot.
    for (i=1; i<n; i++) {

        // This is the row we are simplifying.
        for (j=i; j<n; j++) {

            // This is the factor we multiply each term by in equation i-1.
            double mult = aug_matrix[j][i-1]/aug_matrix[i-1][i-1];

            // Update all entries on row j using the multiplicative factor.
            for (k=0; k<n+1; k++)
                aug_matrix[j][k] -= (mult*aug_matrix[i-1][k]);
        }
    }
}

// Pre-condition: aug_matrix is in lower diagonal form.
// Post-condition: Solves the system, setting aug_matrix to be the identity
void finishSolvingmakeLowerDiag(double aug_matrix[][SIZE+1], int n) {

    int i,j,k;

    // i = row that is currently being solved for.
    for (i=n-1; i>=0; i--) {

        // Sum up the known terms on the LHS of this equation to backtrack.
        double knownterms = 0;
        for (j=i+1; j<n; j++)
            knownterms += (aug_matrix[i][j]*aug_matrix[j][n]);

        // Solve for variable i, using this equation.
        double ans = (aug_matrix[i][n] - knownterms)/aug_matrix[i][i];

        // Reset this entire row to be what you know it should be.
        aug_matrix[i][n] = ans;
        for (j=0; j<n; j++) {
            if (i!=j)
                aug_matrix[i][j] = 0;
            else
                aug_matrix[i][j] = 1;
        }
    }
}