#include <iostream>
#include <cmath>

using namespace std;

struct Vec
{
    double x1;
    double x2;
    double x3;
    Vec(double k1, double k2, double k3)
    {
        x1 = k1; x2 = k2; x3 = k3;
    }
};

Vec F(double *X)
{
    Vec v(0, 0, 0);
    v.x1 = X[0] + 3 * log(X[0]) - X[1] * X[1];
    v.x2 = 2 * X[0] * X[0] - 5 * X[0] - X[0] * X[1] + 1;
    v.x3 = cos(X[0] * X[1]) - X[1] + 0.5;
    return v;
}


struct Matr
{
    double x11;
    double x12;
    double x21;
    double x22;
    double x31;
    double x32;
    Matr(double k1, double k2, double k3, double k4, double k5, double k6)
    {
        x11 = k1; x21 = k2; x31 = k3;
        x12 = k4; x22 = k5; x32 = k6;
    }
};

//Возвращает линеаризованную систему
Matr LinF(double *X)
{
    Matr m(0, 0, 0, 0, 0, 0);
    m.x11 = 1 + 3 / X[0];
    m.x12 = -2 * X[1];
    m.x21 = 4 * X[0] - 5 - X[1];
    m.x22 = -X[0];
    m.x31 = -X[1] * sin(X[0] * X[1]);
    m.x32 = -X[0] * sin(X[0] * X[1]) - 1;
    return m;
}

void RevCopy(double **A, double *B, int n) //Функция из 3.3
{
    for (int i = 0; i < n; ++i)
        B[i] = A[i][i];
}

// Бесконечная норма матрицы
double InfNorm(double **A, int n)
{
    double max = 0;
    double sum = 0;
    for (int i = 0; i < n; ++i)
        max += abs(A[0][i]);

    for (int i = 1; i < n; ++i)
    {
        sum = 0;
        for (int j = 0; j < n; ++j)
            sum += abs(A[i][j]);
        if (max < sum)
            max = sum;
    }
    return max;
}

void Copy(double **A, double **B, int n)
{
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            B[i][j] = A[i][j];
}

void LU(double **A, double *b, double **L, int m)
{
    for (int i = 0; i < m; ++i)
        L[i][i] = 1;

    for (int i = 0; i < m - 1; ++i) //Вычисление LU-разложения.
    {
        for (int j = i + 1; j < m; ++j) //Вычисление масштабирующих множителей.
            L[j][i] = A[j][i] / A[i][i];

        for (int j = i + 1; j < m; ++j) //Вычитание i-ой строки из последующих.
            for (int k = i; k < m; ++k)
                A[j][k] = A[j][k] - L[j][i] * A[i][k];
    }
    for (int i = 0; i < m - 1; ++i)
    {
        for (int j = i + 1; j < m; ++j)
            b[j] = b[j] - L[j][i] * b[i];
    }
}

//Вычисляет решение системы Ax=b
void Solution(double **A, double *b, int m)
{
    double **L = new double *[m];
    for(int i = 0; i < m; i++)
    {
        L[i] = new double[m];
    }

    LU(A, b, L, m); //Разложение матрицы.
    double sum;

    for (int i = m - 1; i >= 0; --i)
    {
        sum = 0;
        for (int j = i + 1; j < m; ++j)
            sum += A[i][j] * A[j][j];
        A[i][i] = (b[i] - sum) / A[i][i];
    }
}

//Транспонированная матрица
void Transp(double **A, double **AT, int m, int n)
{
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < m; ++j)
            AT[i][j] = A[j][i];
}

//Произведение матриц: A nxm, B mxk, P nxk
void MProizv(double **A, double **B, double **P, int n, int m, int k)
{
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < k; ++j)
        {
            P[i][j] = 0;
            for (int r = 0; r < m; ++r)
                P[i][j] += A[i][r] * B[r][j];
        }
}

//Произведение матрицы на вектор: A nxm, b m, P m
void VProizv(double **A, double *B, double *P, int n, int m)
{
    for (int i = 0; i < n; ++i)
    {
        P[i] = 0;
        for (int r = 0; r < m; ++r)
            P[i] += A[i][r] * B[r];
    }
}

//Евклидова норма вектора
double ENorm(double *A, int n)
{
    double sum = 0;
    for (int i = 0; i < n; ++i)
        sum += A[i] * A[i];
    sum = sqrt(sum);
    return sum;
}

//Решение переопределённой системы
void Pereopred(double eps)
{
    int m = 3;
    int n = 2;
    double X[] = { 1, 1 }; //Начальное приближение
    double *delt = new double[n]; //Вектор невязки

    double **A = new double * [m]; //Якобиан в точке Х0
    for(int i = 0; i < m; i++)
    {
        A[i] = new double[n];
    }

    double **AT = new double *[n];
    for(int i = 0; i < m; i++)
    {
        AT[i] = new double[n];
    }
    double **PA = new double *[n]; //Произведение АT*А
    for(int i = 0; i < n; i++)
    {
        PA[i] = new double[n];
    }

    double *Pb = new double[n]; //Произведение АT*b
    double *b = new double[m]; // Правая часть

    // Matr M(0, 0, 0, 0, 0, 0);
    Matr M = LinF(X);
    A[0][0] = M.x11;
    A[0][1] = M.x12;
    A[1][0] = M.x21;
    A[1][1] = M.x22;
    A[2][0] = M.x31;
    A[2][1] = M.x32;

    // Vec v = new Vec(0, 0, 0);
    Vec v = F(X);
    b[0] = v.x1;
    b[1] = v.x2;
    b[2] = v.x3;

    Transp(A, AT, m, n);
    MProizv(AT, A, PA, n, m, n);
    VProizv(AT, b, Pb, n, m);
    Solution(PA, Pb, n);
    RevCopy(PA, delt, n);

    X[0] = X[0] - delt[0];
    X[1] = X[1] - delt[1];

    while (ENorm(delt,n) > eps)
    {
        M = LinF(X);
        A[0][0] = M.x11;
        A[0][1] = M.x12;
        A[1][0] = M.x21;
        A[1][1] = M.x22;
        A[2][0] = M.x31;
        A[2][1] = M.x32;

        v = F(X);
        b[0] = v.x1;
        b[1] = v.x2;
        b[2] = v.x3;

        Transp(A, AT, m, n);
        MProizv(AT, A, PA, n, m, n);
        VProizv(AT, b, Pb, n, m);
        Solution(PA, Pb, n);
        RevCopy(PA, delt, n);

        X[0] = X[0] - delt[0];
        X[1] = X[1] - delt[1];
    }

    cout << "Решение переопределённой системы:" << endl;
    cout << "x = " << X[0] << endl;
    cout << "y = " << X[1] << endl;
    v = F(X);
    cout << "F1 = " << v.x1 << endl;
    cout << "F2 = " << v.x2 << endl;
    cout << "F3 = " << v.x3 << endl;
}

int main()
{
    double eps = 0.000001;
    Pereopred(eps);

    return 0;
}