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

#define N 4

// Function to print a 4x4 matrix
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. Initialize data
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. Gaussian Elimination Results ---\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];
}
}
}
}

// Output Gauss-Jordan Results
printf("--- 2. Gauss-Jordan Method Results ---\n");
for (int i = 0; i < N; i++) {
printf("x%d = %10.4f\n", i + 1, A_GJ[i][N]);
}
printf("\n");

// ==========================================
// 4. Inverting Matrix A (using Gauss-Jordan)
// ==========================================
float A_Inv[N][N] = {
{1, 0, 0, 0},
{0, 1, 0, 0},
{0, 0, 1, 0},
{0, 0, 0, 1}
};
// Initialize identity operations in a copy of matrix 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. Calculated Inverse Matrix (A^-1) ---\n");
printMatrix(A_Inv);

// ==========================================
// 5. Verification (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. Verification (A * A^-1) ---\n");
printMatrix(Identity_Check);

return 0;
}