using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Lab5._3
{
    class Program
    {
         struct Vec
        {
            public double x1;
            public double x2;
            public double x3;

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

         struct Matr
         {
             public double x11;
             public double x12;
             public double x21;
             public double x22;
             public double x31;
             public double x32;

             public 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;
             }
         }

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

         static private Vec F(double[] X)
         {
             Vec v = new Vec(0, 0, 0);
             v.x1 = X[0] + 3 * Math.Log(X[0]) - X[1] * X[1];
             v.x2 = 2 * X[0] * X[0] - 5 * X[0] - X[0] * X[1] + 1;
             v.x3 = Math.Cos(X[0] * X[1]) - X[1] + 0.5;
             return v;
         }

        //Возвращает линеаризованную систему 
         static private 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] * Math.Sin(X[0] * X[1]);
             m.x32 = -X[0] * Math.Sin(X[0] * X[1]) - 1;
             return m;
         }

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

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

        static private 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];
        }

        static private 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
        static private 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
        static private 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
        static private 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];
            }
        }

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

        //Решение переопределённой системы
        static private 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];
            }

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

        static void Main(string[] args)
        {
            double eps = 0.000001;
            Pereopred(eps);
    
        }
    }
}