#include <iostream>
#include <iomanip>
#include <cmath>

using namespace std;

const double eps = 1e-12; // otoczenie zera wynoszące 10^-12

// Funkcja dokonuje rozkładu LU macierzy A opisana uprzednio w algorytmie

bool rozklad(int n, double A[4][4])
{
    int i,j,k; //dane pomocnicze do przechodzenia macierzy

    for(k = 0; k < n - 1; k++)
    {
        if(fabs(A[k][k]) < eps) return false;

        for(i = k + 1; i < n; i++)
        {
            A[i][k] /= A[k][k]; //normalizacja kolumny
        }

        for(i = k + 1; i < n; i++)
        {
            for(j = k + 1; j < n; j++)
            {
                 A[i][j] -= A[i][k] * A[k][j];//modyfikacja podmacirzy
            }
        }
    }

    return true;
}

// Funkcja wyznacza wektor X na podstawie A i B

bool rozwiaz(int n, double A[4][4], double B[4], double * X)
{
    int    i,j;
    double s; // zmienna przechowywująca sumę iloczynów

    X[0] = B[0]; //  wyliczamy wyrazy wektora Y

    for(i = 1; i < n; i++)
    {
        s = 0; // zerujemy sumę iloczynow

        for(j = 0; j < i; j++)
        {
            s += A[i][j] * X[j]; // obliczamy sumę iloczynów l(i,j) x y(j)
        }
        X[i] = B[i] - s; // obliczamy y(i)
    }

    if(fabs(A[n-1][n-1]) < eps) return false; // błąd dzielenia przez 0, zakończ

    X[n-1] /= A[n-1][n-1]; // obliczamy wektor X

    for(i = n - 2; i >= 0; i--)
    {
        s = 0; // zerujemy sumę

        for(j = i + 1; j < n; j++)
        {
            s += A[i][j] * X[j]; // obliczamy sumę iloczynów u(i,j) x x(j)
        }
        if(fabs(A[i][i]) < eps) return false; // blad dzielenia przez 0

        X[i] = (X[i] - s) / A[i][i]; //  obliczamy x(i)
    }

    return true; // zmienna boolowska przyjmuje wartosc true - sukces
}

// Program główny

int  main()
{
    double *X;
    int n,i,j;

    cout << setprecision(4) << fixed; // podanie wyników z dokaldnością do 4 miejsc po przecinku

    n=4;  // liczba niewiadomych w rownaniu

    // wczytujemy macierze macierze A, B oraz inicjujemy macierz X

    X = new double [n];

    // macierz A wspolczynnikow stojacyh przy niewiadomych

    double A[4][4] =
    {
        {2.0, 3.0, 2.0, 4.0},
        {1.0, 2.0, -3.0, -2.0},
        {3.0, 2.0, -2.0, 1.0},
        {-1.0, 4.0, 1.0, -5.0}
    };

    // macierz B wyrazow wolnych
    double B[4] = {8.0,-4.0,2.0,-5.0};

    // rozwiązujemy układ i wyświetlamy wyniki
    cout<<"Układ równań rozwiązany metodą LU ma następujące wyniki: "<<endl<<endl;

    if(rozklad(n,A) && rozwiaz(n,A,B,X))
    {
        for(i = 0; i < n; i++)
            cout << "x" << i + 1 << " = " << X[i] << endl;
    }
    else cout << "DZIELNIK ZERO\n"; // blad dzielenia przez 0


    return 0;
}