import numpy as np


def ternary_search(a, b, x, func, grad, H):
    while b - a > 1e-3:
        c = a + (b - a) / 3
        d = a + (b - a) / 3 * 2
        if func(x - c * H @ grad(x)) >= func(x - d * H @ grad(x)):
            a = c
        else:
            b = d
    return (b + a) / 2


def bfgs(x0, grad, func, step):
    xk = x0.copy()
    n = xk.shape[0]
    H = np.eye(n)
    res = [func(xk)]
    for i in range(1, step):
        hk = ternary_search(0.0000001, 1, xk, func, grad, H)
        x_last = xk.copy()
        xk = xk - hk * H @ grad(xk)
        yk = grad(xk) - grad(x_last)
        delta = xk - x_last
        Bk = 1 + (yk @ delta) / ((H @ yk) @ yk)
        H = H + (H * yk * delta.T + delta * yk.T * H) / ((H @ yk) @ yk)
        H -= Bk * ((H * yk * yk.T * H) / ((H @ yk) @ yk))
        res.append(func(xk))
    return xk, res