def f(x):
    return x ** 3 - 2 * x + 1


def solve(a, b, epsilon):
    p = a + 0.382 * (b - a)
    q = a + 0.618 * (b - a)
    phip = f(p)
    phiq = f(q)
    while True:
        if phip <= phiq:
            if abs(b - p) <= epsilon:
                return a, q
                break
            else:
                b = q
                phiq = phip
                q = p
                p = a + 0.382 * (b - a)
                phip = f(p)
        if phip > phiq:
            if abs(b - p) <= epsilon:
                return p, b
                break
            else:
                a = p
                phip = phiq
                p = q
                q = a + 0.618 * (b - a)
                phiq = f(q)


if __name__ == '__main__':
    a, b = solve(0, 3, 5e-10)
    print('a: {}\nb: {}\nf((a+b)/2): {}'.format(a, b, f((a + b) / 2)))