/*--------------------------------------------------------------------*/ 
/*  problem chapter7_22                                               */ 
/*                                                                    */ 
/*  This program computes a cofactor of a square matrix               */ 
/*  with 4 rows and 4 columns.                                        */ 

#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;    


/* Declare function prototypes */ 
double cofactor(double a[][4], int i, int j);


int main()
{ 
   /*  Declare variables */ 
   int i, j; 
   double a[4][4], value; 

   cout << "Compute a cofactor of 4 x 4 matrx. " << endl;
   /* input a 4 x 4 matrix */
   cout << "Enter a 4 x 4 matrix." <<endl;
   for (i = 0; i<4; i++)
     for (j = 0; j<4; j++)
       {  cin >> a[i][j];
       }

   /* output the matrix for check  */ 
   cout<<"Original matrix is"<<endl;
   for (i = 0; i<4; i++) 
      {
      for (j = 0; j<4; j++) 
        {
        cout<<setw(6)<<a[i][j];
        }
      cout<<endl;
      }

   /* input the (i,j) location where the cofactor will be computed. */
   cout<<"Enter the values of i and j: "<<endl;
   cin >> i >> j;

   value=cofactor(a, i, j);
   
   cout << i << "," << j << " cofactor is: "<< value << endl;
 
   return 0; 
} 


/*-------------------------------------------------------------*/
/*  This function returns the (i,j) cofactor of a matrix       */
/*-------------------------------------------------------------*/


double cofactor(double a[][4], int i, int j)
{
  double b[3][3], value;
 //Create new matrix.
 for (int k=0; k<4; k++) 
      for (int l=0; l<4; l++) 
            { 
               if (  k<i && l<j )   b[k][l]=a[k][l];
               if (  k>i && l<j )   b[k-1][l]=a[k][l];
               if (  k<i && l>j )   b[k][l-1]=a[k][l];
               if (  k>i && l>j )   b[k-1][l-1]=a[k][l];
            }
 
 /* compute the determinant of matrix b */
  value = b[0][0]*b[1][1]*b[2][2] + b[0][1]*b[1][2]*b[2][0]
        + b[0][2]*b[1][0]*b[2][1] - b[2][0]*b[1][1]*b[0][2]
        - b[2][1]*b[1][2]*b[0][0] - b[2][2]*b[1][0]*b[0][1];

  return (value*pow(-1.0, (double)(i+j)));

}
 
/*--------------------------------------------------------------------*/