brent.hpp

Go to the documentation of this file.

/*******************************************************************************
 * Copyright (C) 2017-2026 Theodore Chang
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ******************************************************************************/
 
#ifndef BRENT_HPP
#define BRENT_HPP
 
#include <suanPan.h>
 
template<std::floating_point T, std::invocable<T> F> T brent(F&& func, const T x1, const T x2, const T tol) {
    static constexpr auto eps = std::numeric_limits<T>::epsilon();
    auto a = x1, b = x2, c = x2;
    T fa = func(a), fb = func(b), fc = fb;
 
    T d{}, e{}, p{}, q{};
 
    auto counter = 0u;
    while(true) {
        ++counter;
 
        if(fb * fc > T(0)) {
            c = a;
            fc = fa;
            e = d = b - a;
        }
        if(std::fabs(fc) < std::fabs(fb)) {
            a = b;
            b = c;
            c = a;
            fa = fb;
            fb = fc;
            fc = fa;
        }
        const auto comp_tol = T(2) * eps * std::fabs(b) + T(.5) * tol;
        const auto xm = T(.5) * (c - b);
        if(std::fabs(xm) <= comp_tol || std::fabs(fb) < tol) {
            suanpan_debug("Brent's method initial guess {:.5E} with {} iterations.\n", b, counter);
            return b;
        }
        if(std::fabs(e) >= comp_tol && std::fabs(fa) > std::fabs(fb)) {
            const auto s = fb / fa;
            if(a == c) {
                p = T(2) * xm * s;
                q = T(1) - s;
            }
            else {
                q = fa / fc;
                const auto r = fb / fc;
                p = s * (T(2) * xm * q * (q - r) - (b - a) * (r - T(1)));
                q = (q - T(1)) * (r - T(1)) * (s - T(1));
            }
            if(p > T(0)) q = -q;
            else p = std::fabs(p);
            if(T(2) * p < std::min(T(3) * xm * q - std::fabs(comp_tol * q), std::fabs(e * q))) {
                e = d;
                d = p / q;
            }
            else e = d = xm;
        }
        else e = d = xm;
        a = b;
        fa = fb;
        const auto dd = xm >= T(0) ? comp_tol : -comp_tol;
        fb = func(b += std::fabs(d) > comp_tol ? d : dd);
    }
}
 
#endif