行列式を計算する方法はいろいろありますが、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');
}
}