import sympy as s


def sub_x(f_, val):
    return f_.subs(x, val)


def sub_y(f_, val):
    return f_.subs(y, val)


def sub_h(f_, val):
    return f_.subs(h, val)


def sub(f_, val_x, val_y):
    return sub_y(sub_x(f_, val_x), val_y)


x, y, h = s.symbols('x y h')

f = (x - 3)**2 + (y - 2)**2
# f = 2*x*y + 2*x - x**2 - 2*y**2

f_x = s.diff(f, x)

f_y = s.diff(f, y)


def steepest(x_initial, y_initial, max_error, iterations):
    x0, y0 = x_initial, y_initial
    ite = 0
    error = 1

    while error > max_error and ite < iterations:
        ite += 1

        X = x0 + sub(f_x, x0, y0) * h
        Y = y0 + sub(f_y, x0, y0) * h

        if not(h in X.free_symbols) and not(h in Y.free_symbols):
            break

        g = sub(f, X, Y)
        g_h = s.diff(g, h)

        h_root = s.solve(g_h)[0]

        x_root = sub_h(X, h_root)
        y_root = sub_h(Y, h_root)

        error_x = abs(x_root - x0) / x_root
        error_y = abs(y_root - y0) / y_root

        error = max(error_x, error_y)

        x0, y0 = x_root, y_root

    return [float(x0), float(y0), float(sub(f, x0, y0)), error, ite]


res = steepest(1, 1, 1e-2, 100)

print('The optimum of f(x) is', res[2], 'at x =', res[0], 'and y =', res[1])