#include <iostream>
#include <stdio.h>
#include <math.h>

using namespace std;

  class Program
    {
         struct Vec
        {
		public:
			double x1;
            double x2;
            double x3;

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

        struct Matr
         {
		 public:
			 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;
             }
         };

	public:

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

          Vec* F(double X[])
         {
             Vec* v = new Vec(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;
         }

        //Возвращает линеаризованную систему 
          Matr* LinF(double X[])
         {
             Matr* m = new Matr(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;
         }

        //Бесконечная норма матрицы
         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, 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];
            }
        }

        //Транспонированная матрица
        static private 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, n]; //Якобиан в точке Х0
            double* AT = new double[n, m];
            double* PA = new double[n, n]; //Произведение АT*А
            double* Pb = new double[n]; //Произведение АT*b
            double* b = new double[m]; // Правая часть

            Matr* M = new Matr(0, 0, 0, 0, 0, 0);
            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);
            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<<"Result:" <<"\n";
            cout<<"x = " << X[0] <<"\n";
            cout<<"y = " << X[1] <<"\n";
            v = F(X);
            cout<<"F1 = " << v->x1 <<"\n";
            cout<<"F2 = " << v->x2 <<"\n";
            cout<<"F3 = " << v->x3 <<"\n";
        }
};	
        void Main()
        {
            double eps = 0.000001;
			Program Test;
            Test.Pereopred(eps);
   
        }