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

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

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

#include 
#include 

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');
    }
}

参考文献