#include <iostream>
using namespace std;

void print_r(double **mat, int n, int m) {
	for (int i=0;i<n;i++){
		for (int j=0;j<m;j++) {
			cout << mat[i][j]<< "  ";
		}
		cout << endl;
	}
}
//swapping two rows
void swap_row(double **mat, int n, int m, int i, int j)
{
	cout << "Swapped rows " << i << "<->" << j << endl;
	for (int k=0;k<m;k++) {
		double temp=mat[i][k];
		mat[i][k]=mat[j][k];
		mat[j][k]=temp;
	}
}
//forward elimination
//reduce matrix --> ref
int felim(double **mat, int n, int m) {
	for (int k=0;k<n;k++) {
	//initialize max value and index of pivot
		int i_max=k;
		int v_max=mat[i_max][k];
		//find the greatest amplitude for pivot if 
		//any
		for (int i=k+1;i<n;i++) {
			if (abs(mat[i][k])>v_max) {
				v_max=mat[i][k],i_max=i;
			}
		}
		//if a principal diagonal element is zero,
		// it denotes that matrix is singular
		//and will lead to a division by zero later
		if (!mat[k][i_max]) {
			return k;
		}
		/*swap the  greatest value row with current
		row*/
		if (i_max!=k) {
			swap_row(mat,n,m,k,i_max);
		}
		/*factor f to set current row kth element
		to 0 and remaining column to 0*/
		for (int i=k+1;i<n;i++) {
			double f=mat[i][k]/mat[k][k];
			/*subtract fth multiple of corr. kth
			row element*/
			for (int j=k+1;j<m;j++) {
				mat[i][j]-=mat[k][j]*f;
			}
			//filling the row triangular with 0
			mat[i][k]=0;
		}
	}
	return -1;
}
//function to calculate values of uknowns
double *backusb(double **mat, int n, double *y) {
	/*start calculation of last eq back to first*/
	for (int i=n-1;i>=0;i--) {
	//start with th rhs of eq.
		y[i]=mat[i][n];
	//Initialize j to i+1 since matrix is upper triang.
		for (int j=i+1;j<n;j++) {
			cout <<"(" << i <<"," <<j <<")"<<",";
			/*subtract all the lhs values
			except the coefficient of the variabl
			whose value is being calculated*/
			y[i]-=mat[i][j]*y[j];
		}
		cout << endl;
		/*divide the rhs coeff */
		y[i]=y[i]/mat[i][i];
	}	
	return y;
}
	
//Determinant calculation of matrix n,n
//#include <iostream>
//#include <boost/chrono.hpp>
//using namespace std;
//function to get cofactor of mat[p][q] in temp[][].n
void getCofactor(int **mat, int **temp, int p, int q, int n){
	int i=0,j=0;
	for (int row=0;row<n;row++) {
		for (int col=0;col<n;col++) {
			//Copy into temp  the elements
			//which are not in given row,column
			if (row!=p && col!=q) {
				temp[i][j++]=mat[row][col];
			
			//row is filled so increase row index
			// reset col index
				if (j==n-1) {
					j=0;
					i++;
				}
			}
		}
	}
}
/*Recursive function to get det of matrix*/
int determinantOfMatrix(int **mat, int n) {
	int D=0;
	if (n==1) 
		return mat[0][0];
	int **temp=new int*[n]; //to store the cofactors
	for (int i=0;i<n;i++) 
		temp[i]=new int[n];
	int sign=1; // sign multiplier
	//iterate for each element of first row
	for (int f=0;f<n;f++) {
		getCofactor(mat,temp,0,f,n);
		D+=sign*mat[0][f]*
			determinantOfMatrix(temp,n-1);
		sign=-sign;
	}
	return D;
}
void adjoint(int **,int **, int);
//Function to calculate inverse
//returns false if singular matrix
bool inverse(int **a, float **inv,int n) {
	//find determinant
	int det=determinantOfMatrix(a,n);
	if (det==0) {
		cout << "Mi antistrepsimos matrix !" <<endl;
		return false;
	}
	//find adjoint
	int **adj;
	for (int i=0;i<n;i++) {
		adj=new int *[n];
	}
	for (int i=0;i<n;i++) {
		adj[i]=new int[n];
	}
	adjoint(a,adj,n);
	//find inverse
	for (int i=0;i<n;i++) {
		for (int j=0;j<n;j++) {
			inv[i][j]=adj[i][j]/float(det);
		}
	}
	return true;
}

void print_ri(int **mat, int row, int col) {
	for (int i=0;i<row;i++) {
		for (int j=0;j<col;j++) {
			cout << mat[i][j] << " " ;
		}
		cout << endl;
	}
} 
void print_rf(float **mat, int row, int col) {
	for (int i=0;i<row;i++) { 
		for (int j=0;j<col;j++) {
			cout <<mat[i][j] << "  ";
		}
		cout << endl;
	}
}
//Adjoint
void adjoint(int **a, int **adj, int n) {
	if (n==1) {
		adj[0][0]=1;
		return;
	}
	int sign=1; int **temp;
	for (int i=0;i<n;i++) {
		temp=new int*[n];
	}
	for (int i=0;i<n;i++) {
		temp[i]=new int[n];
	}
	//Cofactor
	for (int i=0;i<n;i++) {
		for (int j=0;j<n;j++) {
			getCofactor(a,temp,i,j,n);
			sign=((i+j)%2==0)?1:-1;
			adj[j][i]=(sign)*
				determinantOfMatrix(temp,n-1);
		}
	}
	for (int i=0;i<n;i++) 
		delete []temp[i];

}