import math
import sys
from sympy import diff, symbols, sin
import numpy

ideal = (1 + math.sqrt(5)) / 2


def bisection(f, a, b, eps, i=0):
    c = (a + b) / 2
    dist = math.fabs(a - b)
    i += 1
    if dist < 2 * eps:
        return c, i
    elif f(c) < f(c + eps):
        return bisection(f, a, c, eps, i)
    else:
        return bisection(f, c, b, eps, i)


def helper(b, a):
    global ideal
    return (a - b) / ideal


def golden(f, x0, x1, eps, i=0):
    alpha = x1 - helper(x0, x1)
    betha = x0 + helper(x0, x1)
    i += 1
    if math.fabs(x1 - x0) < eps:
        return (x0 + x1) / 2, i
    elif f(alpha) >= f(betha):
        return golden(f, alpha, x1, eps, i)
    else:
        return golden(f, x0, betha, eps, i)


def newton():
    x, y = symbols('x y')
    func = x ** 2 + 10 * (y - sin(x)) ** 2
    f_x = diff(func, x)
    f_y = diff(func, y)
    f_xx = diff(func, x, x)
    f_xy = diff(func, x, y)
    f_yx = diff(func, y, x)
    f_yy = diff(func, y, y)
    val = [1, 1]

    def grad(a, b):
        s = {x: a, y: b}
        return numpy.array([f_x.evalf(subs=s), f_y.evalf(subs=s)], dtype=numpy.float64)

    def gesse(a, b):
        s = {x: a, y: b}
        return numpy.array([[f_xx.evalf(subs=s), f_xy.evalf(subs=s)],
                            [f_yx.evalf(subs=s), f_yy.evalf(subs=s)]],
                           dtype=numpy.float64)

    p = numpy.linalg.solve(gesse(1, 1), grad(1, 1))
    next = val-(numpy.linalg.solve(gesse(1, 1), grad(1, 1)))


print(newton())