BrentMethod

#include <algorithm>
#include <iostream>
#include <vector>
#include <cmath>
#include <exception>


double parab(double x1, double f1, double x2, double f2, double x3, double f3)
{
    return x2 - ((x2 - x1) * (x2 - x1) * (f2 - f3) - (x2 - x3) * (x2 - x3) * (f2 - f1)) / (2 * ((x2 - x1) * (f2 - f3) - (x2 - x3) * (f2 - f1)));
}

double  parabola(double a, double b, double (*f)(double x), int64_t n, double eps)
{
    if (a > b)
        std::swap(a, b);
    double e = (a + b) * 0.5;
    double x = (a + b) * 0.5;
    double res = x;
    double incor = b + 1;
    double c, d;
    for (int64_t i = 0; i < n; ++i)
    {
        res = parab(a, f(a), e, f(e), b, f(b));
        if (!(res >= a && res <= b))
            throw std::invalid_argument("Can't find min");
        if (abs(res - x) < eps)
            return res;
        x = res;
        c = std::min(e, x);
        d = std::max(e, x);
        if (f(c) < f(d))
        {
            b = d;
            e = c;
        }
        else
        {
            a = c;
            e = d;
        }
    }

    return res;
}

double brent(double a, double b, double (*f)(double x), int64_t n, double eps)
{
    if (a > b)
        std::swap(a, b);
    double x, v, w, u, dcur, dprev;
    static const double r = (3 - sqrt(5)) * 0.5;
    x = w = v = a + r * (b - a);
    dcur = dprev = b - a;
    for (int64_t i = 0; i < n; ++i) {
        if (std::max(x - a, b - x) < eps) {
            return x;
        }

        double g = dprev * 0.5;
        dprev = dcur;
        if (x == w || w == v || x == v) {
            if (x < (a + b) * 0.5) {
                u = x + r * (b - x);
                dprev = b - x;
            } else {
                u = x - r * (x - a);
                dprev = x - a;
            }
        } else {
            u = parab(x, f(x), w, f(w), v, f(v));
            if (!(u > a && u < b) || std::abs(u - x) > g) {
                if (x < (a + b) * 0.5) {
                    u = x + r * (b - x);
                    dprev = b - x;
                } else {
                    u = x - r * (x - a);
                    dprev = x - a;
                }
            }
        }

        dcur = std::abs(u - x);
        if (f(u) > f(x)) {
            if (u < x) {
                a = u;
            } else {
                b = u;
            }

            if (f(u) <= f(w) || w == x) {
                v = w;
                w = u;
            } else {
                if (f(u) <= f(v) || v == x || v == w) {
                    v = u;
                }
            }
        } else {
            if (u < x) {
                b = x;
            } else {
                a = x;
            }
            v = w;
            w = x;
            x = u;
        }

    }

    return x;
}