#include <iostream>
#include <cstdlib>
#include <string>
using namespace std;


class Matrix {
private:
    long double** MATRIX;
    int M, N;
    long double** Create()
    {
        MATRIX = new long double* [M];
        for (int z = 0; z < M; z++)
            MATRIX[z] = new long double[N];
        return MATRIX;
    }
public:

    Matrix(int i, int j) : M(i), N(j)   // конструктор матрицы с нулевыми элементами
    {
        Create();
        for (int i = 0; i < M; i++)
        {
            for (int j = 0; j < N; j++)
            {
                MATRIX[i][j] = 0;
            }
        }
    }
    ~Matrix()    // деструктор матрицы
    {
        for (int i = 0; i < M; i++)
        {
            delete[] MATRIX[i];
        }
        delete[] MATRIX;
    }
    void setJK(int j, int k, long double r) { MATRIX[j][k] = r; } // ввод kj элемента матрицы
    int giveM() { return M; }  //вывод количества строк
    int giveN() { return N; }  //вывод количества столбцов
    long double giveJK(int j, int k) { return MATRIX[j][k]; }

    void CreateByHand() {   // заполнение вручную
        int i, j;
        for (i = 0; i < M; i++) {
            for (j = 0; j < N; j++) {
                cout << "[" << i << "][" << j << "]= ";
                cin >> MATRIX[i][j];
            }

        }
        cout << " " << endl;
    }
    void CreateRandom() {  // рандомное заполнение
        int i, j;
        for (i = 0; i < M; i++) {
            for (j = 0; j < N; j++) {
                MATRIX[i][j] = 100 * (double)(rand()) / RAND_MAX;
            }
        }
    }
    long double Element(int j, int m)  // элемент с номером j, m
    {
        return MATRIX[j][m];
    }

    void Print() // функция вывода матрицы
    {
        int i, j; // счетчики
        for (i = 0; i < M; i++) {
            for (j = 0; j < N; j++)
                cout << MATRIX[i][j] << " ";
            cout << endl;
        }
    }
};



// Функция копирования матрицы
Matrix EqualsMatrix(Matrix A)
{
    Matrix F(A.giveM(), A.giveN());
    int i, j;
    for (i = 0; i < A.giveM(); i++) {
        for (j = 0; j < A.giveN(); j++) {
            F.setJK(i, j, A.giveJK(i, j));
        }
    }
    return  F;
}



// Получение матрицы без i-й строки и j-го столбца
Matrix GetMatr(Matrix A, int i, int j) {
    Matrix G(A.giveM() - 1, A.giveN() - 1);
    int ki, kj, di, dj;
    di = 0;
    for (ki = 0; ki < A.giveM() - 1; ki++) { // проверка индекса строки
        if (ki == i) di = 1;
        dj = 0;
        for (kj = 0; kj < A.giveN() - 1; kj++) { // проверка индекса столбца
            if (kj == j) dj = 1;
            G.setJK(ki, kj, A.giveJK(ki + di, kj + dj));
        }
    }
    return G;
}



// Рекурсивное вычисление определителя
long double Determinant(Matrix A) {
    int i, k;
    long double d;
    d = 0;
    k = 1; //(-1) в степени i
    if (A.giveM() == 1) {
        d = A.giveJK(0, 0);
        return(d);
    }
    if (A.giveM() == 2) {
        d = A.giveJK(0, 0) * A.giveJK(1, 1) - (A.giveJK(1, 0) * A.giveJK(0, 1));
        return(d);
    }
    if (A.giveM() > 2) {
        for (i = 0; i < A.giveM(); i++) {
            Matrix J = GetMatr(A, i, 0);
            d = d + k * A.giveJK(i, 0) * Determinant(J);
            k = -k;
        }
    }
    return(d);
}



// промежуточные матрицы с замененными стобцами
Matrix PartialMatrix(Matrix V, Matrix K, int t)
{
    int z;
    Matrix F = EqualsMatrix(V);
    for (z = 0; z < V.giveM(); z++) {
        F.setJK(z, t, K.giveJK(z, 0));
    }
    return F;
}



// Функция получения решения методом краммера
void GetsolutionKrammer(Matrix A, Matrix B)
{
    long double u = Determinant(A);
    cout << " " << endl << "Determinant of the system matrix equals " << u << endl;
    if (u == 0) {
        cout << " " << endl << "System cannot be solved";
    }
    else {

        int e;
        Matrix D(A.giveM(), 1); // вектор решения
        for (e = 0; e < A.giveM(); e++) {
            Matrix C = PartialMatrix(A, B, e);
            cout << " " << endl << "Partial matrix " << e << " is" << endl << " " << endl;
            C.Print();
            long double Y = Determinant(C);
            cout << " " << endl << " Determinant of partial matrix " << e << " is  " << Y << endl;
            D.setJK(0, e, Y / u);

        }
        for (e = 0; e < A.giveM(); e++) {
            cout << " " << endl << "Solution x" << e << " :" << D.Element(0, e) << endl;
        }
    }
}



// Функция задания системы
void SetsystemKrammer(Matrix A, Matrix B)
{
    string ans; // спрашиваем, заполнять все вручную или рандомно
    cout << "Would you like to use random system matrix and constant terms, or would you like to insert them manually? : " << endl;
    cout << "1.Manually" << endl << "2.Random" << endl;
    cin >> ans;
    if (ans == "1") {
        cout << "Type in system matrix: " << endl;
        A.CreateByHand();
        A.Print();
        cout << " " << endl << "Type in constant terms: " << endl;
        B.CreateByHand();
        B.Print();
    }
    else {
        A.CreateRandom();
        cout << " " << endl << "System matrix is: " << endl << " " << endl;
        A.Print();
        B.CreateRandom();
        cout << " " << endl << "Constant terms are: " << endl << " " << endl;
        B.Print();
    }
}


// Функция склеивания матриц
Matrix GlueMatrixes(Matrix A, Matrix B)
{
    Matrix L(A.giveM(), A.giveN() + B.giveN());
    int i, j;
    for (i = 0; i < A.giveM(); i++) {
        for (j = 0; j < A.giveN(); j++) {
            L.setJK(i, j, A.giveJK(i, j));
        }
    }
    for (i = 0; i < A.giveM(); i++) {
        for (j = A.giveN(); j < L.giveN(); j++) {
            L.setJK(i, j, B.giveJK(i, j - A.giveN()));
        }
    }
    return L;
}



// Решение Методом Гаусса
Matrix SolutionGauss(Matrix L)
{
    int i, j, k;
    Matrix U = EqualsMatrix(L);;
    for (i = 0; i < L.giveM(); i++)  // счетчик переменной, по которой производится вычитание
    {
        for (j = 0; j < i; j++) // счетчик строки
        {
            long double f = U.giveJK(j, i) / U.giveJK(i, i);
            for (k = 0; k < L.giveN(); k++) // счетчик столбца
            {
                U.setJK(j, k, (U.giveJK(j, k) - U.giveJK(i, k) * f));
            }
        }

        for (j = i + 1; j < L.giveM(); j++) // счетчик строки
        {
            long double f = U.giveJK(j, i) / U.giveJK(i, i);
            for (k = 0; k < L.giveN(); k++) // счетчик столбца
            {
                U.setJK(j, k, (U.giveJK(j, k) - U.giveJK(i, k) * f));
            }

        }

        if (i < U.giveM() - 1) {

            long double h = 0;
            for (k = 0; k < L.giveN(); k++)
            {
                h = U.giveJK(i + 1, k) * U.giveJK(i + 1, k); // проверяем, не получилась ли следующая строка (содержащая следующую переменную) нулевой
            }
            if (h == 0)
            {
                i = i + 1;  // если получилась, переходим к следующей
            }
        }
    }
    return U;
}

// Функция, выводящая результат решения системы методом гаусса
void GiveSolutionGauss(Matrix U)
{
    int i;
    long double solution;
    for (i = 0; i < U.giveM(); i++)
    {
        solution = U.giveJK(i, U.giveN() - 1) / U.giveJK(i, i);
        cout << " " << endl << "Solution x" << i << " :" << solution << endl;
    }
}

void SolveByGauss(Matrix A, Matrix B)
{
    Matrix L = GlueMatrixes(A, B);
    Matrix U = SolutionGauss(L);
    GiveSolutionGauss(U);
}


// Основная функция
int main() {

    int cycle = 1;
    while (cycle == 1) // выполняем программу несколько раз
    {
        int m;
        cout << "Enter the size of the system matrix: ";
        cin >> m;
        Matrix A(m, m); // матрица системы
        Matrix B(m, 1); // вектор свободных членов
        SetsystemKrammer(A, B);
        GetsolutionKrammer(A, B);
        int t;
        if (Determinant(A) != 0) {
            cout << " " << endl << "Would you like to use Gauss method on this system?" << endl;
            cout << "1.Yes" << endl << "2.No" << endl;
            cin >> t;
            if (t == 1) { SolveByGauss(A, B); }
        }
        cout << " " << endl << "Would you like to continue?" << endl;
        cout << "1.Yes" << endl << "2.No" << endl;
        cin >> cycle;
    }
    return 0;
}