// OM_1.cpp : Этот файл содержит функцию "main". Здесь начинается и заканчивается выполнение программы.
//

#include <iostream>
using namespace std;

//double k = 1, N = 0;
const int M = 50;
int k = 0;
using gradient = pair <double, double>;
using dot = pair <double, double>;
using vec = pair <double, double>;
using matrix = pair <dot, dot>;
double F(dot x) {
    return x.first * x.first + 8 * x.second * x.second + x.first * x.second + x.first;
}
gradient Gradient(dot x) {
    return make_pair(2 * x.first + x.second + 1, 16 * x.second + x.first);
}
matrix Hessian() {
    return  make_pair(make_pair(2, 1), make_pair(1, 16));
}
matrix InverseHessian() {
    return  make_pair(make_pair(16.0 / 31, -1.0 / 31), make_pair(-1.0 / 31, 2.0 / 31));
}
double Norm(dot x) {
    return sqrt(x.first * x.first + x.second * x.second);
}
vec Multiply(matrix m, vec v) {
    return make_pair(m.first.first * v.first + m.first.second * v.second, m.second.first* v.first + m.second.second * v.second);
}
double Checkpositive(matrix m) {

    double d1 = m.first.first;
    double d2 = m.first.first * m.second.second - m.first.second * m.second.second;
    return d1 > 0 && d2 > 0;
}
//}
double φ(double t, const dot x, gradient grad) {
    double first = t * grad.first;
    double second = t * grad.second;
    first = x.first - first;
    second = x.second - second;
    return F(make_pair(first, second));
}
//a = 0, b = 5, ε = 0.01, l = 0.1
//a=0 потому что шаги назад делать не надо
double Dichotomy(double a, double b, double eps, double l, const dot x, gradient grad) {
    double yk, zk, φyk, φzk; int counter=0;
    while (true) {
        yk = (a + b - eps) / 2;
        φyk = φ(yk,x,grad);
        zk = (a + b + eps) / 2;
        φzk = φ(zk,x, grad);
        if (φyk <= φzk)b = zk;
        else a = yk;
        double dif = abs(a - b);
        if (dif <= l) return (a + b) / 2;
        else  counter++;
    }
}
//
//
//dot GradientDescent(dot x, double eps1, double eps2) {
//    int k = 0; bool flag = false;
//
//    while (true) {
//        gradient grad = Gradient(x);
//        if (Norm(grad) < eps1)
//            return x;
//        if (k >= M)
//            return x;
//        double t = Dichotomy(0, 5, 0.01, 0.1, x, grad);
//
//        dot oldx = x;
//
//        double first = t * grad.first;
//        double second = t * grad.second;
//        first = x.first - first;
//        second = x.second - second;
//        x.first = first; x.second = second;
//
//        if (Norm(make_pair(x.first - oldx.first, x.second - oldx.second)) < eps2 && abs(f(x) - f(oldx)) < eps2) {
//
//            if (flag == true)
//                return x;
//            flag = true;
//
//        }
//        else {
//            k++;
//            flag = false;
//        }
//    }
//}

dot Newton(dot x0, double eps1, double eps2) {

    dot x = x0;
    double t; dot newx;
    double oldF = F(x);
    bool flag = false;
    while (true) {
        gradient grad = Gradient(x);
        if (Norm(grad) < eps1 || k >= M)
            return x;
        matrix ihessian = InverseHessian();
        vec d = make_pair(-grad.first, -grad.second);
        if (Checkpositive(ihessian)) {
            d = Multiply(ihessian, d);
            //newx = make_pair(x.first + d.first, x.second + d.second);
        }
        t = Dichotomy(0, 5, 0.01, 0.1, x, grad);
        newx.first = t * d.first; newx.second = t * d.second;
        newx.first += x.first; newx.second += x.second;
        double newF = F(x);
        /*double newF = F(x);
        if (t * Norm(d) < eps2 && abs(oldF - newF) < eps2) {
            if (flag) return x;
            else flag = true;
        }
        else flag = false;
        oldF = newF;
        k += 1;*/

        if (Norm(make_pair(newx.first - x.first, newx.second - x.second)) < eps2 && abs(F(newx) - F(x)) < eps2) {


            if (flag)
                return x;
            flag = true;

            }
            else {
                k++;
                flag = false;
            }
        x = newx;
        oldF = newF;
    }
}

int main()
{
    //dot x = GradientDescent(make_pair(1.5, 0.5), 0.1, 0.15);
    //cout << x.first << "\t" << x.second;
    dot x = Newton(make_pair(1.5, 0.5), 0.1, 0.15);
    cout << x.first << "\t" << x.second << endl;
    cout << k;
}

// Запуск программы: CTRL+F5 или меню "Отладка" > "Запуск без отладки"
// Отладка программы: F5 или меню "Отладка" > "Запустить отладку"

// Советы по началу работы
//   1. В окне обозревателя решений можно добавлять файлы и управлять ими.
//   2. В окне Team Explorer можно подключиться к системе управления версиями.
//   3. В окне "Выходные данные" можно просматривать выходные данные сборки и другие сообщения.
//   4. В окне "Список ошибок" можно просматривать ошибки.
//   5. Последовательно выберите пункты меню "Проект" > "Добавить новый элемент", чтобы создать файлы кода, или "Проект" > "Добавить существующий элемент", чтобы добавить в проект существующие файлы кода.
//   6. Чтобы снова открыть этот проект позже, выберите пункты меню "Файл" > "Открыть" > "Проект" и выберите SLN-файл.