/****************************************************************************************
	filename: matrix.h	name:KIHA
	行列に関する関数集
****************************************************************************************/
#include<stdio.h>

#define MAT_RANK 10

/***************************************************************************************/
/*	プロトタイプ宣言一覧							       */
/***************************************************************************************/
/* 行列の出力 */
void print_int_matrix(int, int, int[][]);
void print_double_matrix(int, int, double[][]);

/* 逆行列 */
int inverse_matrix(int, double[][], double[][]);
int inverse_int_matrix_rank2(int[][], int[][]);
int inverse_double_matrix_rank2(double[][], double[][]);
int inverse_int_matrix_rank2_modulo(int[][], int[][], int);

/* 行列の乗算 */
int multiple_int_matrix(int, int, int[][], int, int, int[][], int[][]);
int multiple_double_matrix(int, int, double[][], int, int, double[][], double[][]);
int multiple_int_matrix_rank2(int[][], int[][], int[][]);

/* 行列式 */
double determinant(int, int, double[][]);		
int determinant_int_rank2(int[][]);
/***************************************************************************************/


/* 行列を表示する関数 */
void print_int_matrix(int row, int column, int matrix[MAT_RANK][MAT_RANK])
{
	int i, j;

	printf("--------------------------------------\n");
	for(i=0;i<row;i++){
		for(j=0;j<column;j++){
			printf("%10d ", matrix[i][j]);
		}
		printf("\n");
	}
	printf("--------------------------------------\n");
}


/* 行列を表示する関数(小数点以上 6 文字，以下 3 文字) */
void print_double_matrix(int row, int column, double matrix[MAT_RANK][MAT_RANK])
{
	int i, j;

	printf("--------------------------------------\n");
	for(i=0;i<row;i++){
		for(j=0;j<column;j++){
			printf("%10.3lf ", matrix[i][j]);
		}
		printf("\n");
	}
	printf("--------------------------------------\n");
}


/* 逆行列を求める関数．						*/
/* 第２引数で入力された行列の逆行列を第３引数に与える関数．	*/
/* 逆行列が存在すれば１を返し，存在しなければ０を返す．		*/
/* determinant関数を使用					*/
int inverse_matrix(int rank, double matrix1[MAT_RANK][MAT_RANK], double matrix2[MAT_RANK][MAT_RANK])
{
	double hakidashi[MAT_RANK][MAT_RANK*2];
	double waru, kakeru;
	double det;
	int i, j, k;
	int success_flag;

	det = determinant(rank, rank, matrix1);
	if(det < 0.0001 && det > -0.0001)	/* 丸め誤差の範囲 */
	{
		success_flag = 0;
	}else{
		/* 初期化 */
		for(i=0;i<rank;i++){
			for(j=0;j<rank;j++){
				hakidashi[i][j] = matrix1[i][j];
			}
			for(j=rank;j<rank*2;j++){
				if(j-rank == i){
					hakidashi[i][j] = 1;
				}
				else{
					hakidashi[i][j] = 0;
				}
			}
		}
	
		/* 掃き出し法により逆行列を求める */
		for(i=0;i<rank;i++){
			waru = hakidashi[i][i];
			for(j=i;j<rank*2;j++){
				hakidashi[i][j] = hakidashi[i][j] / waru;
			}
			for(k=0;k<rank;k++){
				if(i==k){
					continue;
				}else{
					kakeru = hakidashi[k][i];
					for(j=i;j<rank*2;j++){
						hakidashi[k][j] = hakidashi[k][j] - hakidashi[i][j] * kakeru;
					}
				}
			}
		}

		/* 結果を第３引数に代入 */
		for(i=0;i<rank;i++){
			for(j=0;j<rank;j++){
				matrix2[i][j] = hakidashi[i][j+rank];
			}
		}
		success_flag = 1;
	}
	return(success_flag);
}


/* 二次正方行列(整数)の逆行列を求める関数．			*/
/* 第１引数で入力された行列の逆行列を第２引数に与える関数．	*/
/* 逆行列が存在すれば１を返し，存在しなければ０を返す．		*/
/* determinant関数を使用					*/
int inverse_int_matrix_rank2(int matrix1[MAT_RANK][MAT_RANK], int matrix2[MAT_RANK][MAT_RANK])
{
	int det;
	int success_flag;

	det = determinant_int_rank2(matrix1);
	if(det < 0.0001 && det > -0.0001)	/* 丸め誤差の範囲 */
	{
		success_flag = 0;
	}else{
		matrix2[0][0] = matrix1[1][1] / det;
		matrix2[0][1] = (-1) * matrix1[0][1] / det;
		matrix2[1][0] = (-1) * matrix1[1][0] / det;
		matrix2[1][1] = matrix1[0][0] / det;
	}
	return(success_flag);
}


/* 二次正方行列(実数)の逆行列を求める関数．			*/
/* 第１引数で入力された行列の逆行列を第２引数に与える関数．	*/
/* 逆行列が存在すれば１を返し，存在しなければ０を返す．		*/
/* determinant関数を使用					*/
int inverse_double_matrix_rank2(double matrix1[MAT_RANK][MAT_RANK], double matrix2[MAT_RANK][MAT_RANK])
{
	double det;
	int success_flag;

	det = determinant(2, 2, matrix1);
	if(det == 0)
	{
		success_flag = 0;
	}else{
		matrix2[0][0] = matrix1[1][1] / det;
		matrix2[0][1] = (-1) * matrix1[0][1] / det;
		matrix2[1][0] = (-1) * matrix1[1][0] / det;
		matrix2[1][1] = matrix1[0][0] / det;
	}
	return(success_flag);
}


/* 二次正方行列(整数)の逆行列を求める関数(法バージョン)．		*/
/* 第１引数で入力された行列の逆行列を第２引数に与える関数．	*/
/* 逆行列が存在すれば１を返し，存在しなければ０を返す．		*/
/* determinant関数を使用					*/
int inverse_matrix_rank2_modulo(int matrix1[MAT_RANK][MAT_RANK], int matrix2[MAT_RANK][MAT_RANK], int modulo)
{
	int det;
	int success_flag;

	det = determinant_int_rank2(matrix1);
	if(det == 0)
	{
		success_flag = 0;
	}else{
		det = det % modulo;
		matrix2[0][0] = (matrix1[1][1] / det) % modulo;
		matrix2[0][1] = (((-1) * matrix1[0][1] / det) % modulo) + modulo;
		matrix2[1][0] = (((-1) * matrix1[1][0] / det) % modulo) + modulo;
		matrix2[1][1] = (matrix1[0][0] / det) % modulo;
	}
	return(success_flag);

}


/* 行列(整数)の乗算を行う関数．						*/
/* 第３引数と第６引数の行列乗算を行い、演算結果を第７引数に与える関数．	*/
/* 乗算が可能であれば１を返し，不可能であれば０を返す			*/
int multiple_int_matrix(int row1, int column1, int matrix1[MAT_RANK][MAT_RANK], 
			int row2, int column2, int matrix2[MAT_RANK][MAT_RANK], 
			int matrix3[MAT_RANK][MAT_RANK])
{
	int i, j, k;
	int success_flag;

	if(column1 != row2){
		printf("演算不能(multiple_matrix関数内)：行数と列数が合っていません\n");
		success_flag = 0;
	}else{
		/* 初期化 */
		for(i=0;i<row1;i++){
			for(j=0;j<column2;j++){
				matrix3[i][j] = 0;
			}
		}

		/* 乗算の計算 */
		for(i=0;i<row1;i++){
			for(j=0;j<column2;j++){
				for(k=0;k<column1;k++){
					matrix3[i][j] = matrix3[i][j] + matrix1[i][k] * matrix2[k][j];
				}
			}
		}
		success_flag = 1;
	}
	
	return(success_flag);
}


/* 行列(実数)の乗算を行う関数．						*/
/* 第３引数と第６引数の行列乗算を行い、演算結果を第７引数に与える関数．	*/
/* 乗算が可能であれば１を返し，不可能であれば０を返す			*/
int multiple_double_matrix(int row1, int column1, double matrix1[MAT_RANK][MAT_RANK], 
		    int row2, int column2, double matrix2[MAT_RANK][MAT_RANK], 
		    double matrix3[MAT_RANK][MAT_RANK])
{
	int i, j, k;
	int success_flag;

	if(column1 != row2){
		printf("演算不能(multiple_matrix関数内)：行数と列数が合っていません\n");
		success_flag = 0;
	}else{
		/* 初期化 */
		for(i=0;i<row1;i++){
			for(j=0;j<column2;j++){
				matrix3[i][j] = 0;
			}
		}

		/* 乗算の計算 */
		for(i=0;i<row1;i++){
			for(j=0;j<column2;j++){
				for(k=0;k<column1;k++){
					matrix3[i][j] = matrix3[i][j] + matrix1[i][k] * matrix2[k][j];
				}
			}
		}
		success_flag = 1;
	}
	
	return(success_flag);
}


/* ２次正方行列(整数)の乗算を行う関数．					*/
/* 第３引数と第６引数の行列乗算を行い、演算結果を第７引数に与える関数．	*/
/* 乗算が可能であれば１を返し，不可能であれば０を返す			*/
int multiple_int_matrix_rank2(int matrix1[MAT_RANK][MAT_RANK], int matrix2[MAT_RANK][MAT_RANK], int matrix3[MAT_RANK][MAT_RANK])
{
	int i, j, k;
	int success_flag;

	/* 初期化 */
	for(i=0;i<2;i++){
		for(j=0;j<2;j++){
			matrix3[i][j] = 0;
		}
	}

	/* 乗算の計算 */
	for(i=0;i<2;i++){
		for(j=0;j<2;j++){
			for(k=0;k<2;k++){
				matrix3[i][j] = matrix3[i][j] + matrix1[i][k] * matrix2[k][j];
			}
		}
	}
	success_flag = 1;
	
	return(success_flag);
}


/* 行列式の計算結果を返す関数 */
double determinant(int row, int column, double matrix[MAT_RANK][MAT_RANK])
{
	double soinshi[MAT_RANK][MAT_RANK];
	double waru, kakeru;
	int i, j, k;
	double ans;
	
	if(row != column){
		printf("error(determinant()):仕様に合っていません\n");
		ans = 0;
	}else{
		/* 初期化 */
		for(i=0;i<row;i++){
			for(j=0;j<column;j++){
				soinshi[i][j] = matrix[i][j];
			}
		}
		ans = 1;
	
		/* 素因子展開により行列式の解を求める */
		for(i=0;i<row;i++){
			waru = soinshi[i][i];
			ans = ans * waru;
			for(j=i;j<column;j++){
				soinshi[i][j] = soinshi[i][j] / waru;
			}
			for(k=i+1;k<row;k++){
				kakeru = soinshi[k][i];
				for(j=i;j<column;j++){
					soinshi[k][j] = soinshi[k][j] - soinshi[i][j] * kakeru;
				}
			}
		}
	}
	return(ans);
}


/* 二次正方行列(整数)の行列式を計算する関数． */
int determinant_int_rank2(int matrix[MAT_RANK][MAT_RANK])
{
	return(matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]);
}