n次行列式を計算するプログラム

行列式を計算する方法はいろいろありますが、4次以上の行列については、上三角行列に変形し、対角成分の積をとる方法が有効です。

このプログラムでは、行列式に関する以下の性質を利用しています。

#include <stdio.h>
#include <stdlib.h>

double determinant(double *m, int N);
void inputMatrix(double *mat, int cols, int rows, char *mNameStr);
double *allocateMatrix(int *N);
void printMatrix(double *a, int N);

int main(int argc, char *argv[])
{
	int N;
	double *m;
	
	m = allocateMatrix(&N);
	inputMatrix(m, N, N, "m");
	printf("det = %6.4lf\n", determinant(m, N));

	return 0;
}

#define SWAP(a, b)	(a != b && (a += b, b = a - b, a -= b))
double determinant(double *m, int N)
{
	int x, y, i;
	double det = 1, r;

	// 上三角行列に変換しつつ、対角成分の積を計算する。
	for(y = 0; y < N - 1; y++){
		printMatrix(m, N);
		if(m[y * N + y] == 0){
			// 対角成分が0だった場合は、その列の値が0でない行と交換する
			for(i = y + 1; i < N; i++){
				if(m[i * N + y] != 0){
					break;
				}
			}
			if(i < N){
				for(x = 0; x < N; x++){
					SWAP(m[i * N + x], m[y * N + x]);
				}
				// 列を交換したので行列式の値の符号は反転する。
				det = -det;
			}
		}
		for(i = y + 1; i < N; i++){
			r = m[i * N + y] / m[y * N + y];
			for(x = y; x < N; x++){
				m[i * N + x] -= r * m[y * N + x];
			}
		}
		det *= m[y * N + y];
	}
	det *= m[y * N + y];

	return det;
}

void inputMatrix(double *mat, int cols, int rows, char *mNameStr)
{
	int x, y;
	
	for(y = 0; y < rows; y++){
		for(x = 0; x < cols; x++){
			if(rows != 1){
				printf("%s[%d][%d] = ", mNameStr, y, x);
			} else{	
				printf("%s[%d] = ", mNameStr, x);
			}
			scanf("%lf", &mat[y * cols + x]);
		}
	}
}

double *allocateMatrix(int *N)
{
	printf("N = ");
	scanf("%d", N);
	return malloc(*N * *N * sizeof(double));
}

void printMatrix(double *a, int N)
{
	int i, k;
	for(i = 0; i < N; i++){
		for(k = 0; k < N; k++){
			printf("%+6.2f ", a[i * N + k]);
		}
		putchar('\n');
	}
}

参考文献