#include #include #include #include #include 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; }