#include <stdio.h>
#define N 4

int main() {
double A[N][N] = {
{3, 1.5, -6, 4.8},
{1, 1.5, -2, -2.4},
{0, -1.5, -2, -1},
{2, 4, -1.8, -0.6}
};

double b[N] = {1.2, 0.6, -2.4, 0};
double L[N][N] = {0};
double U[N][N];
double y[N] = {0};
double x[N] = {0};
double LU[N][N] = {0};
double Ax[N] = {0};

int i, j, k;

for(i = 0; i < N; i++) {
for(j = 0; j < N; j++) {
U[i][j] = A[i][j];
}
}

for(i = 0; i < N; i++) {
L[i][i] = 1.0;
}

for(i = 0; i < N - 1; i++) {
for(j = i + 1; j < N; j++) {

double factor = U[j][i] / U[i][i];
L[j][i] = factor;

for(k = 0; k < N; k++) {
U[j][k] = U[j][k] - factor * U[i][k];
}
}
}
for(i = 0; i < N; i++) {

y[i] = b[i];

for(j = 0; j < i; j++) {
y[i] -= L[i][j] * y[j];
}
}

for(i = N - 1; i >= 0; i--) {

x[i] = y[i];

for(j = i + 1; j < N; j++) {
x[i] -= U[i][j] * x[j];
}

x[i] /= U[i][i];
}
printf("L行列:\n");

for(i = 0; i < N; i++) {
for(j = 0; j < N; j++) {
printf("%8.3f ", L[i][j]);
}
printf("\n");
}
printf("\nU行列:\n");

for(i = 0; i < N; i++) {
for(j = 0; j < N; j++) {
printf("%8.3f ", U[i][j]);
}
printf("\n");
}

printf("\nyベクトル:\n");

for(i = 0; i < N; i++) {
printf("%8.3f\n", y[i]);
}


printf("\n連立方程式の解 x:\n");

for(i = 0; i < N; i++) {
printf("x%d = %8.3f\n", i + 1, x[i]);
}

for(i = 0; i < N; i++) {
for(j = 0; j < N; j++) {

LU[i][j] = 0;

for(k = 0; k < N; k++) {
LU[i][j] += L[i][k] * U[k][j];
}
}
}

printf("\nLU行列（LU = A の確認）:\n");

for(i = 0; i < N; i++) {
for(j = 0; j < N; j++) {
printf("%8.3f ", LU[i][j]);
}
printf("\n");
}

for(i = 0; i < N; i++) {

Ax[i] = 0;

for(j = 0; j < N; j++) {
Ax[i] += A[i][j] * x[j];
}
}

printf("\nAxベクトル（Ax = b の確認）:\n");

for(i = 0; i < N; i++) {
printf("%8.3f\n", Ax[i]);
}

printf("\n元のbベクトル:\n");

for(i = 0; i < N; i++) {
printf("%8.3f\n", b[i]);
}
return 0;
}