// Arup Guha
// 9/15/09
// Honors Mathematical Modeling
// An equation solver (ONLY FOR a system of 3 equations with a unique solution.

double det(double mat[][3]);
void replacecol(double mat[][3], double col[], int colnum);

int main() {
    
    // These are hard-coded values from the dem-rep-ind problem in the text.
    // I added the equation r + d + i = 1, so these stand for the fraction of
    // each in the steady state solution.
    double colx[] = {-.25, .05, 1};
    double coly[] = {.2, -.4, 1};
    double colz[] = {.4, .2, 1};
    double colsol[] = {0, 0, 1};
    
    double mat[3][3];
    
    // Fill in initial matrix.
    replacecol(mat, colx, 0);
    replacecol(mat, coly, 1);
    replacecol(mat, colz, 2);
    
    // Calculate its determinant.
    double detMat = det(mat);
    
    // Now, replace the first column (x) with the solution column to get detX.
    replacecol(mat, colsol, 0);
    double detX = det(mat);
    
    // Put back original column 1, and replace column 2 with the solution col.
    replacecol(mat, colx, 0);
    replacecol(mat, colsol, 1);
    double detY = det(mat);
    
    // Put back column 2 and put the replace column 3 with the solution col.
    replacecol(mat, coly, 1);
    replacecol(mat, colsol, 2);
    double detZ = det(mat);
    
    // Finish up Kramer's Rule.
    double x = detX/detMat;
    double y = detY/detMat;
    double z = detZ/detMat;
    
    // Print out the final answer as percentages.
    printf("Rep= %.2lf, Dem=%.2lf, Ind=%.2lf\n", 100*x, 100*y, 100*z);
    
    system("PAUSE");
    return 0;
    
    
}

// Returns the determinant of mat.
// Note: Only works for 3x3 matrix.
double det(double mat[][3]) {
       
   double answer = 0;
   int i;
   
   // Add up forward diagonals in basket-weaving algorithm.
   for (i=0; i<3; i++)     
       answer += mat[0][i]*mat[1][(i+1)%3]*mat[2][(i+2)%3];
       
   // Subtract out backward diagonals.
   for (i=0; i<3; i++)
       answer -= mat[2][i]*mat[1][(i+1)%3]*mat[0][(i+2)%3];
       
   return answer;
       
}

// Replaces the colnum column of mat with col.
void replacecol(double mat[][3], double col[], int colnum) {

    int i;
    
    // Just write over this column (colnum) with col.
    for (i=0; i<3; i++)
        mat[i][colnum] = col[i]; 
}
