#include <iostream>
#include <math.h>
#include <vector>
#include <cmath>
#include <iomanip>


void inversion(double **A, int N)
{
    double temp;

    double **E = new double *[N];

    for (int i = 0; i < N; i++)
        E[i] = new double [N];

    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
        {
            E[i][j] = 0.0;

            if (i == j)
                E[i][j] = 1.0;
        }

    for (int k = 0; k < N; k++)
    {
        temp = A[k][k];

        for (int j = 0; j < N; j++)
        {
            A[k][j] /= temp;
            E[k][j] /= temp;
        }

        for (int i = k + 1; i < N; i++)
        {
            temp = A[i][k];

            for (int j = 0; j < N; j++)
            {
                A[i][j] -= A[k][j] * temp;
                E[i][j] -= E[k][j] * temp;
            }
        }
    }

    for (int k = N - 1; k > 0; k--)
    {
        for (int i = k - 1; i >= 0; i--)
        {
            temp = A[i][k];

            for (int j = 0; j < N; j++)
            {
                A[i][j] -= A[k][j] * temp;
                E[i][j] -= E[k][j] * temp;
            }
        }
    }

    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
            A[i][j] = E[i][j];

    for (int i = 0; i < N; i++)
        delete [] E[i];

    delete [] E;
}

int main() {
    int n = 0;

    std::cout << "Enter matrix dimension (n):" << std::endl;
    std::cin >> n;

    float *diag_a = new float[n];
    float *diag_b = new float[n];
    float *diag_c = new float[n];
    float *vect_d = new float[n];
    float *vect_x = new float[n];

    float *vect_d1 = new float[n];
    float *vect_r = new float[n];
    float *vect_e = new float[n];

    float *alpha = new float[n];
    float *beta = new float[n];

//    int N;
//    std::cout << "Enter N: ";
//    std::cin >> N;
    double  **A = new double *[n];

    std::cout << "Enter main diagonal (a):" << std::endl;
    for (int i = 0; i < n; i++) {
        std::cin >> diag_a[i];
    }

    std::cout << "Enter side diagonal 1 (b):" << std::endl;
    for (int i = 0; i < n - 1; i++) {
        std::cin >> diag_b[i];
    }

    std::cout << "Enter side diagonal 2 (c):" << std::endl;
    for (int i = 0; i < n - 1; i++) {
        std::cin >> diag_c[i];
    }

    std::cout << "Enter side free vector (d):" << std::endl;
    for (int i = 0; i < n; i++) {
        std::cin >> vect_d[i];
    }


    //Check diagonal prevalence conditions

    for (int i = 1; i < n - 1; i++) {
        if (abs(diag_a[i]) < abs(diag_c[i - 1] + diag_b[i])) {
            std::cout << "The first condition of diagonal dominance is not satisfied!" << std::endl;
            return 0;
        } else std::cout << "Ok. The first condition of diagonal dominance is satisfied!" << std::endl;


        if ((abs(diag_b[i]) / abs(diag_a[i])) > 1) {
            std::cout << "The second condition of diagonal dominance is not satisfied!" << std::endl;
            return 0;
        } else std::cout << "Ok. The second condition of diagonal dominance is satisfied!" << std::endl;

        if ((abs(diag_c[i - 1]) / abs(diag_b[i])) > abs(diag_c[i - 1] + diag_b[i])) {
            std::cout << "The third condition of diagonal dominance is not satisfied!" << std::endl;
            return 0;
        } else std::cout << "Ok. The third condition of diagonal dominance is satisfied!" << std::endl;
    }


    //Find vector x
    if (diag_a[0] != 0) {
        alpha[0] = -(diag_b[0] / diag_a[0]);
        beta[0] = vect_d[0] / diag_a[0];
    } else return 0;


    for (int i = 1; i < n; i++) {
        alpha[i] = -diag_b[i] / (alpha[i - 1] * diag_c[i - 1] + diag_a[i]);
        beta[i] = (vect_d[i] - diag_c[i - 1] * beta[i - 1]) / (alpha[i - 1] * diag_c[i - 1] + diag_a[i]);
    }

    vect_x[n - 1] = beta[n - 1];

    for (int i = n - 2; i >= 0; i--) {
        vect_x[i] = (alpha[i] * vect_x[i + 1] + beta[i]);
        //vect_x[i] = (vect_d[i] - diag_c[i - 1] * beta[i - 1]) / (diag_a[i] + (diag_c[i - 1] * alpha[i - 1]));
    }

    std::cout << std::endl;
    std::cout << "Vector x (x*)" << std::endl;

    for (int i = 0; i < n; i++) {
        std::cout << std::fixed << std::setprecision(15) << vect_x[i] << std::endl;
    }

//    std::cout << std::endl;
//    std::cout << "Vector d*" << std::endl;
//
//    for (int i = 0; i < n; i++) {
//        vect_d1[i] = diag_a[i] * vect_x[i];
//        std::cout << std::fixed << std::setprecision(15) << vect_d1[i] << std::endl;
//
//    }



    std::cout << "Enter matrix A:" << std::endl;

    for (int i = 0; i < n; i++) {
        A[i] = new double[n];
    }

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            std::cin >> A[i][j];
        }
    }

    std::cout << std::endl;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            vect_d1[i] = vect_d1[i] + (A[i][j] * vect_x[i]);
        }
        //std::cout << vect_d1[i] << std::endl;
        std::cout << std::fixed << std::setprecision(15) << "d*" << vect_d1[i] << std::endl;

    }

    std::cout << std::endl;
    std::cout << "Vector r" << std::endl;

    for (int i = 0; i < n; i++) {
        vect_r[i] = vect_d[i] - vect_d1[i];
        std::cout << std::fixed << std::setprecision(15) << vect_r[i] << std::endl;

    }

    std::cout << "Matrix A^-1:" << std::endl;
    inversion(A, n);
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
            std::cout << A[i][j] << "  ";

        std::cout << std::endl;
    }

    std::cout << std::endl;
    std::cout << "Vector e" << std::endl;

    std::cout << std::endl;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            vect_e[i] = vect_e[i] + (A[i][j] * vect_r[i]);
        }
        //std::cout << vect_d1[i] << std::endl;
        std::cout << std::fixed << std::setprecision(15) << "e" << vect_e[i] << std::endl;

    }
}