//
//  main.cpp
//  lab1
//
//  Created by Anastasia on 14.02.2020.
//  Copyright © 2020 Anastasia. All rights reserved.
//

#include <cmath>
#include <iostream>

void inversion(float** A, int N);

int main(int argc, const char * argv[]) {
    int n = 0;

    std::cout << "Enter size of matrix \n";
    std::cin >> n;

    float *diag_a = new float[n];
    float *diag_b = new float[n-1];
    float *diag_c = new float[n-1];
    float *free_d = new float[n];

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

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

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

    std::cout << "Enter free numbers d \n";
    for (int i=0; i<n; i++) {
        std::cin >> free_d[i];
    }

    float *alpha = new float[n-1];
    float *beta = new float[n-1];
    float *x = new float[n];

    //-------------------------------
    for (int i=1; i<n-1; i++) {
        if (std::abs(diag_a[i]) < std::abs(diag_c[i-1] + diag_b[i])) {
            std::cout << "Input is incorrect 1 stat ";
        }
        if (std::abs(diag_b[i])/(std::abs(diag_a[i])) > 1) {
            std::cout << "Input is incorrect 2 stat";
        }
        if (std::abs(diag_c[i-1]) / (std::abs(diag_b[i])) > 1) {
            std::cout << "Input is incorrect 3 stat";
        }
    }

    if (diag_a[0] != 0) {
        alpha[0] = - diag_b[0]/diag_a[0];
        std::cout << "alpha[" << "0" <<"]= " << alpha[0];
    } else {
        std::cout << "a[0] = 0 error";
        return 1;
    }
    beta[0] = free_d[0] / diag_a[0];
    std::cout << "beta[" << "0" <<"]= " << beta[0];

    for (int i=1; i<n-1; i++) {
        alpha[i] = -(diag_b[i])/(alpha[i-1]*diag_c[i-1] + diag_a[i]);
        beta[i] = (free_d[i]-diag_c[i-1]*beta[i-1])/(alpha[i-1]*diag_c[i-1] + diag_a[i]);
        std::cout << "alpha[" << i <<"]= " << alpha[i] << " ";
        std::cout << "beta[" << i <<"]= " << beta[i] << " ";
    }
    beta[n-1] = (free_d[n-1]-diag_c[n-2]*beta[n-2])/(alpha[n-2]*diag_c[n-2] + diag_a[n-1]);
    std::cout << "beta[" << n-1 <<"]= " << beta[n-1] << " ";

    alpha[n-1] = -(diag_b[n-1-1])/(alpha[n-1-1]*diag_c[n-1-1] + diag_a[n-1-1]);
    //std::cout << " --- "<< diag_b[n-1] << " " << alpha[n-1-1] << " " << diag_c[n-1-1] << " " << diag_a[n-1] << " --- ";
    std::cout << "alpha[" << n-1 <<"]= " << alpha[n-1] << std::endl;
    x[n-1] = beta[n-1];

    for (int i=n-2; i>=0; i--) {
        x[i] = alpha[i]*x[i+1] + beta[i];
    }

    for (int i=0; i<n; i++) {
        std::cout << "x[" << i << "]= " << x[i] << std::endl;
    }

    //---------2------------
    float *r = new float[n];
    float *e = new float[n];
    float *d = new float[n];

    float **matrix = new float*[n];
    for (int i = 0; i < n; i++) {
        matrix[i] = new float[n];
    }


    matrix[0][0] = 3; matrix[0][1] = 1; matrix[0][2] = 0; matrix[0][3] = 0;
    matrix[1][0] = 1; matrix[1][1] = 3; matrix[1][2] = 1; matrix[1][3] = 0;
    matrix[2][0] = 0; matrix[2][1] = 1; matrix[2][2] = 3; matrix[2][3] = 1;
    matrix[3][0] = 0; matrix[3][1] = 0; matrix[3][2] = 1; matrix[3][3] = 3;


    //matrix[0][0] = 0.2679; matrix[0][1] = -0.0717; matrix[0][2] = 0.0191; matrix[0][3] = -0.0047;
    //matrix[1][0] = -0.0717; matrix[1][1] = 0.2870; matrix[1][2] = -0.0765; matrix[1][3] = 0.0191;
    //matrix[2][0] = 0.0191; matrix[2][1] = -0.0765; matrix[2][2] = 0.2870; matrix[2][3] = -0.0717;
    //matrix[3][0] = -0.0047; matrix[3][1] = 0.0191; matrix[3][2] = -0.0717; matrix[3][3] = 0.2679;





    for (int i=0; i<n; i++) {
        for (int j=0; j<n; j++) {
            d[i] += matrix[i][j] * x[j];
        }
        std::cout << "d[" <<i<<"]= " << d[i] << std::endl;
    }

    for(int i=0; i<n; i++) {
        r[i] = free_d[i]-d[i];
        std::cout << "r[" <<i<<"]= " << r[i] << std::endl;
    }

    inversion(matrix, n);

    for (int i=0; i<4; i++) {
        for (int j=0; j<4; j++) {
            std::cout << matrix[i][j] << " ";
        }
        std::cout << std::endl;
    }

    for (int i=0; i<4; i++) {
        //sum = 0;
        for (int j=0; j<n; j++) {
            e[i]+= matrix[i][j] * r[j];
        }
    }

    std::cout << "error: " << std::endl;
    for (int i=0; i<4; i++) {
        printf("%.15f", e[i]);
        std::cout << " " << std::endl;
    }

    for (int i=0; i<4; i++) {
        float c =x[i]-e[i];
        std::cout << "x[" << i+1 << "]= ";
        printf("%.15f", c);
        std::cout << std::endl;
    }

    return 0;
}

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

        float **E = new float*[N];
        for (int i = 0; i < N; i++) {
            E[i] = new float[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];

    }