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

#define N 4

// 4x4行列を表示する関数
void printMatrix(float matrix[N][N]) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            printf("%10.4f ", matrix[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}

int main() {
    // 1. データの初期化
    float A_GE[N][N + 1] = {
        {3.0,  1.5, -6.0,  4.8,  1.2},
        {1.0,  1.5, -2.0, -2.4,  0.6},
        {0.0, -1.5, -2.0, -1.0, -2.4},
        {2.0,  4.0, -1.8, -0.6,  0.0}
    };

    float A_GJ[N][N + 1] = {
        {3.0,  1.5, -6.0,  4.8,  1.2},
        {1.0,  1.5, -2.0, -2.4,  0.6},
        {0.0, -1.5, -2.0, -1.0, -2.4},
        {2.0,  4.0, -1.8, -0.6,  0.0}
    };
    
    float A_Original[N][N] = {
        {3.0,  1.5, -6.0,  4.8},
        {1.0,  1.5, -2.0, -2.4},
        {0.0, -1.5, -2.0, -1.0},
        {2.0,  4.0, -1.8, -0.6}
    };

    float x[N];
    float ratio, sum;

    // ==========================================
    // 2. ガウスの消去法 (Gaussian Elimination)
    // ==========================================
    for (int i = 0; i < N - 1; i++) {
        for (int j = i + 1; j < N; j++) {
            ratio = A_GE[j][i] / A_GE[i][i];
            for (int k = 0; k < N + 1; k++) {
                A_GE[j][k] -= ratio * A_GE[i][k];
            }
        }
    }

    // 後退代入 (Back-Substitution)
    x[N - 1] = A_GE[N - 1][N] / A_GE[N - 1][N - 1];
    for (int i = N - 2; i >= 0; i--) {
        sum = 0;
        for (int j = i + 1; j < N; j++) {
            sum += A_GE[i][j] * x[j];
        }
        x[i] = (A_GE[i][N] - sum) / A_GE[i][i];
    }

    printf("--- 1. ガウスの消去法による計算結果 ---\n");
    for (int i = 0; i < N; i++) {
        printf("x%d = %10.4f\n", i + 1, x[i]);
    }
    printf("\n");

    // ==========================================
    // 3. ガウス・ジョルダン法 (Gauss-Jordan Method)
    // ==========================================
    for (int i = 0; i < N; i++) {
        float pivot = A_GJ[i][i];
        for (int j = 0; j < N + 1; j++) {
            A_GJ[i][j] /= pivot;
        }

        for (int k = 0; k < N; k++) {
            if (k != i) {
                float factor = A_GJ[k][i];
                for (int j = 0; j < N + 1; j++) {
                    A_GJ[k][j] -= factor * A_GJ[i][j];
                }
            }
        }
    }

    // ガウス・ジョルダン法の結果を出力
    printf("--- 2. ガウス・ジョルダン法による計算結果 ---\n");
    for (int i = 0; i < N; i++) {
        printf("x%d = %10.4f\n", i + 1, A_GJ[i][N]);
    }
    printf("\n");

    // ==========================================
    // 4. 逆行列の計算 (ガウス・ジョルダン法を使用)
    // ==========================================
    float A_Inv[N][N] = {
        {1, 0, 0, 0},
        {0, 1, 0, 0},
        {0, 0, 1, 0},
        {0, 0, 0, 1}
    };
    
    // Aの値を操作用の一時行列にコピー
    float A_Temp[N][N];
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            A_Temp[i][j] = A_Original[i][j];
        }
    }

    for (int i = 0; i < N; i++) {
        float pivot = A_Temp[i][i];
        for (int j = 0; j < N; j++) {
            A_Temp[i][j] /= pivot;
            A_Inv[i][j] /= pivot;
        }
        for (int k = 0; k < N; k++) {
            if (k != i) {
                float factor = A_Temp[k][i];
                for (int j = 0; j < N; j++) {
                    A_Temp[k][j] -= factor * A_Temp[i][j];
                    A_Inv[k][j] -= factor * A_Inv[i][j];
                }
            }
        }
    }

    printf("--- 3. 算出された逆行列 (A^-1) ---\n");
    printMatrix(A_Inv);

    // ==========================================
    // 5. 検算 (A * A^-1)
    // ==========================================
    float Identity_Check[N][N] = {0};
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            for (int k = 0; k < N; k++) {
                Identity_Check[i][j] += A_Original[i][k] * A_Inv[k][j];
            }
        }
    }

    printf("--- 4. 検算結果 (A * A^-1) ---\n");
    printMatrix(Identity_Check);

    return 0;
}