#include <iostream>
#include <cmath>

using namespace std;

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

Vec F(double *X)
{
    Vec v(0, 0, 0, 0);
    v.x1 = tan(X[0]*X[1]-0.2) -  pow(X[0] ,2)+3;
    v.x2 = 2 *pow( X[0],4) +pow( X[1],2) -  1;
    v.x3 = X[0] + X[1] -  3;
	v.x4 = 0.3*X[0] - X[1]+exp(0.2*(pow(X[0],2)+(pow(X[1],2))));
    return v;
}


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

//Возвращает линеаризованную систему
Matr LinF(double *X)
{
    Matr m(0, 0, 0, 0, 0, 0, 0, 0);
    m.x11 = 1 /pow(cos(X[1]),2)-2*X[0];
    m.x12 = 1 /pow(cos(X[0]),2);
    m.x21 = 8 *pow(X[0],3);
    m.x22 = 2*X[1];
    m.x31 = 1;
    m.x32 = 1;
	m.x41 = 0.3+0.4*X[0]*exp(0.2*(pow(X[0],2)+(pow(X[1],2))));
    m.x42 = -1+0.4*X[1]*exp(0.2*(pow(X[0],2)+(pow(X[1],2))));
    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 MainElem( int k, double **mas, int n,int m, int otv[] )
{
  int i, j, i_max = k, j_max = k;
  double temp;
  //Ищем максимальный по модулю элемент
  for ( i = k; i < m; i++ )
    for ( j = k; j < n; j++ )
      if ( fabs( mas[i_max] [j_max] ) < fabs( mas[i] [j] ) )
      {
        i_max = i;
        j_max = j;
      }
  //Переставляем строки
  for ( j = k; j < n ; j++ )
  {
    temp = mas[k] [j];
    mas[k] [j] = mas[i_max] [j];
    mas[i_max] [j] = temp;
  }
  //Переставляем столбцы
  for ( i = 0; i < m; i++ )
  {
    temp = mas[i] [k];
    mas[i] [k] = mas[i] [j_max];
    mas[i] [j_max] = temp;
  }
  //Учитываем изменение порядка корней
  i = otv[k];
  otv[k] = otv[j_max];
  otv[j_max] = i;
}

double* OneVector(double **A,double *B, int m,int n)
{
  double* x = new double [n];
  int i, j, k;
  int* otv = new int [n]; //Отвечает за порядок корней
  double **mas = new double *[n];
  for(int i = 0; i < m; i++)
    {
        mas[i] = new double[m+1];
    }
  for ( j = 0; j < m; j++ )
  {
	  for ( i = 0; i < n; i++ )
	  {
		  if (j!=m)
		  {
			  mas[i][j]=A[i][j];
		  }
		  else
		  {
			  mas[i][j]=B[i];
		  }
	  }
  }
  //Сначала все корни по порядку
  for ( i = 0; i < n; i++ )
    otv[i] = i;
  //Прямой ход метода Гаусса
  for ( k = 0; k < n; k++ )
  { //На какой позиции должен стоять главный элемент
    MainElem( k, mas, n,m, otv ); //Установка главного элемента
    for ( j = n; j >= k; j-- )
      mas[k] [j] /= mas[k] [k];
    for ( i = k + 1; i < n; i++ )
      for ( j = n; j >= k; j-- )
        mas[i] [j] -= mas[k] [j] * mas[i] [k];
  }
  //Обратный ход
  for ( i = 0; i < n; i++ )
    x[i] = mas[i] [n];
  for ( i = n - 2; i >= 0; i-- )
    for ( j = i + 1; j < n; j++ )
      x[i] -= x[j] * mas[i] [j];
  return x;
}

//Транспонированная матрица
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 += pow(A[i] ,2);
    sum = sqrt(sum);
    return sum;
}

//Решение переопределённой системы
void Pereopred(double eps)
{
    int m = 4;
    int n = 2;
    double X[] = { 0.5, 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;
	A[3][0] = M.x41;
    A[3][1] = M.x42;

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

    Transp(A, AT, m, n);
    MProizv(AT, A, PA, n, m, n);
    VProizv(AT, b, Pb, n, m);
    delt=OneVector(PA,Pb,m,n);

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

	cout << "x = " << delt[0] << endl;
    cout << "y = " << delt[1] << endl;
    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;
		A[3][0] = M.x41;
        A[3][1] = M.x42;

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

        Transp(A, AT, m, n);
        MProizv(AT, A, PA, n, m, n);
        VProizv(AT, b, Pb, n, m);
        delt=OneVector(PA,Pb,m,n);

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

	cout << "Result:" << 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;
	cout << "F4 = " << v.x4 << endl;
}

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

    return 0;
}