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

double myPersonalNumber = 7 * 0.2;

int matrixSize = 4;
//Мое k, номер в журнале
//8.30   2.62+a   4.10  1.90
//3.92   8.45   7.78-a    2.46
//3.77  7.21+a  8.04  2.28
//2.21  3.65-a   1.69  6.99


//1.Методом Гаусса найти решение слау.вычислить определитель.
//найти обратную матрицу .сделать проверку
//2. Методом квадратного корня найти решение системы
//3.методом простой итерации найти решение. Точность E=10^-5. Число итераций, вектор ???невязки???(я хз, я ее не слышу)
//4. Гаусса-Зейбеля Е=10^-5, число итераций

template <typename T> int sgn(T val) {
    return (T(0) < val) - (val < T(0));
}

using namespace std;


void printMatrix(vector<vector<double>> matrix){
    for(int i = 0; i < matrix.size(); i++){
        for (int j = 0; j < matrix[i].size(); j++) {

            if (j == matrix.size()) {
                cout << left << dec << fixed << "| " << setw(8) << setprecision(3) << matrix[i].back();
            } else {
                cout << left << dec << fixed << setw(8) << setprecision(3) << matrix[i][j];
            }

        }
        cout << endl;

    }
    cout << endl;
}

void printMatrix(vector<vector<double>> matrix, vector<double> vector){
    for(int i = 0; i < matrix.size(); i++){
        for (int j = 0; j <= matrix[i].size(); j++) {
            if (j == matrix.size()) {
                cout << left << dec << fixed << "| " << setw(8) << setprecision(3) << vector[i];
            } else {
                cout << left << dec << fixed << setw(8) << setprecision(3) << matrix[i][j];
            }

        }
        cout << endl;

    }
    cout << endl;
}



int main() {

    vector<vector<double>> givenMatrix{
        {8.30, 2.62 + myPersonalNumber, 4.10,                    1.90, 10.65 + myPersonalNumber},
        {3.92, 8.45,                    7.78 - myPersonalNumber, 2.46, 12.21},
        {3.77, 7.21 + myPersonalNumber, 8.04,                    2.28, 15.45 - myPersonalNumber},
        {2.21, 3.65 - myPersonalNumber, 1.69,                    6.99, -8.35},
    };

    vector<vector<double>> givenMatrixCopy = givenMatrix;


    cout << "GIVEN MATRIX:" << endl;

    printMatrix(givenMatrix);

    cout << "START OF GAUSS" << endl;
    cout << endl;

    vector<vector<double>> transponedGiven = {
           {0, 0, 0, 0},
           {0, 0, 0, 0},
           {0, 0, 0, 0},
           {0, 0, 0, 0},
       };

    for(int i = 0; i < matrixSize; ++i){
        for(int j = 0; j < matrixSize; ++j){
            transponedGiven[j][i] = givenMatrix[i][j];
        }
    }


    double determinant = 1;

    for(int i = 0; i < matrixSize - 1; i++) {
        for(int j = i + 1; j < matrixSize; j++) {
            double temp = givenMatrix[j][i] / givenMatrix[i][i];
            for(int k = 0; k < matrixSize + 1; k++){
                givenMatrix[j][k] -= givenMatrix[i][k] * temp;
            }
        }
    }
    cout << "DIAGONAL:" << endl;
    printMatrix(givenMatrix);


    for(int i = 0; i < matrixSize; i++) {
        determinant *= givenMatrix[i][i];
    }

    cout << "DETERMINANT : " << determinant << endl;
    cout << endl;



    for(int i = 0; i < matrixSize; i++) {
        double temp = givenMatrix[i][i];
        for(int j = i; j < matrixSize + 1; j++) {
            givenMatrix[i][j]/=temp;
        }
    }

    cout << "INIT:" << endl;
    printMatrix(givenMatrix);

    vector<double> answers = {1,1,1,1};

    for ( int i = matrixSize - 1; i >= 0; i-- ) {
        double temp = 0;
        for (int j = i + 1; j < matrixSize; j++ ) {
            temp += givenMatrix[i][j] * answers[j];
        }
        answers[i] = (givenMatrix[i][4] - temp) / givenMatrix[i][i];
    }

    cout << "Solution of the system is : " << endl;

    for(int i = 0; i < matrixSize; i++) {
        cout << "X" << i+1 << ": " << answers[i] << endl;
    }
    cout << endl;
    cout << "END OF GAUSS" << endl;
    cout << endl;


    vector<vector<double>> tempMatrix = givenMatrix;
    vector<vector<double>> inversedMatrix = {
        {0, 0, 0, 0},
        {0, 0, 0, 0},
        {0, 0, 0, 0},
        {0, 0, 0, 0},
    };

    for (int i = 0; i < matrixSize; i++) {
        for (int j = 0; j < matrixSize; j++) {
            if (i != j) {
                inversedMatrix[i][j] = 0;
            }
            else {
                inversedMatrix[i][j] = 1;
            }
        }
    }

    double tempElement = 0;
    for (int i = 0; i < matrixSize; i++) {
        tempElement = tempMatrix[i][i];
        for (int j = 0; j < matrixSize; j++) {
            tempMatrix[i][j] = tempMatrix[i][j] / tempElement;
            inversedMatrix[i][j] = inversedMatrix[i][j] / tempElement;
        }
        for (int j = i + 1; j < matrixSize; j++) {
            tempElement = tempMatrix[j][i];
            for (int k = 0; k < matrixSize; k++) {
                tempMatrix[j][k] = tempMatrix[j][k] - tempMatrix[i][k] * tempElement;
                inversedMatrix[j][k] = inversedMatrix[j][k] - inversedMatrix[i][k] * tempElement;
            }
        }
    }

    for (int i = matrixSize - 1; i >= 0; i--) {
        for (int j = i - 1; j >= 0; j--) {
            tempElement = tempMatrix[j][i];
            for (int k = 0; k < matrixSize; k++) {
                tempMatrix[j][k] = tempMatrix[j][k] - tempMatrix[i][k] * tempElement;
                inversedMatrix[j][k] = inversedMatrix[j][k] - inversedMatrix[i][k] * tempElement;
            }
        }
    }



    cout << "INVERSED:" << endl;
    printMatrix(inversedMatrix);

    vector<vector<double>> proofOfInversibility = {
            {0, 0, 0, 0},
            {0, 0, 0, 0},
            {0, 0, 0, 0},
            {0, 0, 0, 0},
        };

    for (int i = 0; i < matrixSize; i++) {
        for (int j = 0; j < matrixSize; j++) {
            int c = 0;
            for (int k = 0; k < matrixSize; k++)
            {
                c += givenMatrix[i][k] * inversedMatrix[k][j];
            }
            proofOfInversibility[i][j] = c;
        }
    }

    cout << "A * A^(-1):" << endl;
    printMatrix(proofOfInversibility);

    //end of gauss method


    cout << endl;
    cout << "START OF SQUARE ROOT METHOD" << endl;
    cout << endl;


    vector<vector<double>> forSquareMethodMatrix = {
        {0, 0, 0, 0},
        {0, 0, 0, 0},
        {0, 0, 0, 0},
        {0, 0, 0, 0},
    };

    cout << "TRANSPONED" << endl;
    printMatrix(transponedGiven);



    vector<double> toSquareMethodVector = {10.65 + myPersonalNumber, 12.21, 15.45 - myPersonalNumber, -8.35};
    vector<double> forSquareMethodVector = {0, 0, 0, 0};



    for (int i = 0; i < matrixSize; i++) {
        for (int j = 0; j < matrixSize; j++)
        {
            forSquareMethodVector[i] += toSquareMethodVector[j] * transponedGiven[i][j];
        }
    }

    for (int i = 0; i < matrixSize; i++) {
        for (int j = 0; j < matrixSize; j++) {
            double c = 0;
//            cout << "C = ";
            for (int k = 0; k < matrixSize; k++)
            {
//                cout << transponedGiven[k][i] << " * " << givenMatrixCopy[j][k];
                c += transponedGiven[k][i] * givenMatrixCopy[j][k];
//                cout << " + ";
            }
//            cout << " = " << c << endl;
            forSquareMethodMatrix[i][j] = c;
        }
    }

    cout << endl;
    cout << "SYMMETRICAL: (A^T * A)"<< endl;
    printMatrix(forSquareMethodMatrix, forSquareMethodVector);



    double temp;

    vector<double> x = {0, 0, 0, 0};
    vector<double> y = {0, 0, 0, 0};
    vector<vector<double>> UpperTriangle = {
        {0, 0, 0, 0},
        {0, 0, 0, 0},
        {0, 0, 0, 0},
        {0, 0, 0, 0},
    };



    for (int i = 0; i < matrixSize; i++)
    {
        temp = 0;
        for (int k = 0; k < i; k++)
            temp = temp + UpperTriangle[k][i] * UpperTriangle[k][i];
        UpperTriangle[i][i] = sqrt(forSquareMethodMatrix[i][i] - temp);
        for (int j = i; j < matrixSize; j++)
        {
            temp = 0;
            for (int k = 0; k < i; k++)
                temp = temp + UpperTriangle[k][i] * UpperTriangle[k][j];
            UpperTriangle[i][j] = (forSquareMethodMatrix[i][j] - temp) / UpperTriangle[i][i];
        }
    }

    cout << "UPPER TRIANGLE:" << endl;
    printMatrix(UpperTriangle);

    for (int i = 0; i < matrixSize; i++)
    {
        temp = 0;
        for (int k = 0; k < i; k++)
            temp = temp + UpperTriangle[k][i] * y[k];
        y[i] = (forSquareMethodVector[i] - temp) / UpperTriangle[i][i];
    }
    for (int i = matrixSize - 1; i >= 0; i--)
    {
        temp = 0;
        for (int k = i + 1; k < matrixSize; k++)
            temp = temp + UpperTriangle[i][k] * x[k];
        x[i] = (y[i] - temp) / UpperTriangle[i][i];
    }

    cout << "Solution of the system is : " << endl;
    for (int i = 0; i < matrixSize; i++){
        cout << "X" << i + 1 << " = " << x[i] << endl;
    }

    cout << endl;
    cout << "END OF SQUARE ROOT METHOD" << endl;
    cout << endl;










    return 0;
}