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

#define N 8

void print_matrix(double matrix[N][N]) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            printf("%f ", matrix[i][j]);
        }
        printf("\n");
    }
}

double compute_determinant(double matrix[N][N], int n) {
    double sub_matrix[N][N];
    double determinant = 0;
    int i, j, k, sign;

    if (n == 2) {
        return matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];
    } else {
        for (i = 0; i < n; i++) {
            sign = (i % 2 == 0) ? 1 : -1;

            for (j = 1; j < n; j++) {
                for (k = 0; k < n; k++) {
                    if (k < i) {
                        sub_matrix[j - 1][k] = matrix[j][k];
                    } else if (k > i) {
                        sub_matrix[j - 1][k - 1] = matrix[j][k];
                    }
                }
            }

            determinant += sign * matrix[0][i] * compute_determinant(sub_matrix, n - 1);
        }
    }

    return determinant;
}

void compute_inverse(double matrix[N][N], double inverse[N][N]) {
    double determinant = compute_determinant(matrix, N);
    double adjugate[N][N];
    double sub_matrix[N][N];
    int i, j, k, l, sign;

    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            sign = ((i + j) % 2 == 0) ? 1 : -1;

            for (k = 0; k < N - 1; k++) {
                for (l = 0; l < N - 1; l++) {
                    if (k < i) {
                        if (l < j) {
                            sub_matrix[k][l] = matrix[k][l];
                        } else {
                            sub_matrix[k][l] = matrix[k][l + 1];
                        }
                    } else {
                        if (l < j) {
                            sub_matrix[k][l] = matrix[k + 1][l];
                        } else {
                            sub_matrix[k][l] = matrix[k + 1][l + 1];
                        }
                    }
                }
            }

            adjugate[j][i] = sign * compute_determinant(sub_matrix, N - 1);
        }
    }

    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            inverse[i][j] = adjugate[i][j] / determinant;
        }
    }
}

int is_identity_matrix(double matrix[N][N]) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            if ((i == j && (matrix[i][j] < 0.99 || matrix[i][j] > 1.01)) ||
                (i != j && (matrix[i][j] < -0.01 || matrix[i][j] > 0.01))) {
                return 0; // 非単位行列
            }
        }
    }
    return 1; // 単位行列
}

int main() {
    srand(time(NULL));

    double random_matrix[N][N];
    double inverse_matrix[N][N];
    double result[N][N];

    // ランダムな値で初期化された8x8行列の生成
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            random_matrix[i][j] = (double)rand() / RAND_MAX;
        }
    }

    // 逆行列の計算
    compute_inverse(random_matrix, inverse_matrix);

    // 元の行列と逆行列の掛け算
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            double sum = 0;

            for (int k = 0; k < N; k++) {
                sum += random_matrix[i][k] * inverse_matrix[k][j];
            }

            result[i][j] = sum;
        }
    }

    // 結果の表示
    printf("Random Matrix:\n");
    print_matrix(random_matrix);

    printf("\nInverse Matrix:\n");
    print_matrix(inverse_matrix);

    printf("\nResult (Random Matrix * Inverse Matrix):\n");
    print_matrix(result);

    // 単位行列であるかの判断
    int is_identity = is_identity_matrix(result);
    printf("\nIs Result an Identity Matrix? %s\n", is_identity ? "Yes" : "No");

    return 0;
}