Brent.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <cmath>
#include <functional>
#include <utility>
 
#include "NumCpp/Core/DtypeInfo.hpp"
#include "NumCpp/Core/Types.hpp"
#include "NumCpp/Roots/Iteration.hpp"
#include "NumCpp/Utils/essentiallyEqual.hpp"
 
namespace nc::roots
{
    //================================================================================
    // Class Description:
    class Brent : public Iteration
    {
    public:
        //============================================================================
        // Method Description:
        Brent(const double epsilon, std::function<double(double)> f) noexcept :
            Iteration(epsilon),
            f_(std::move(f))
        {
        }
 
        //============================================================================
        // Method Description:
        Brent(const double epsilon, const uint32 maxNumIterations, std::function<double(double)> f) noexcept :
            Iteration(epsilon, maxNumIterations),
            f_(std::move(f))
        {
        }
 
        //============================================================================
        // Method Description:
        ~Brent() override = default;
 
        //============================================================================
        // Method Description:
        double solve(double a, double b)
        {
            resetNumberOfIterations();
 
            double fa = f_(a);
            double fb = f_(b);
 
            checkAndFixAlgorithmCriteria(a, b, fa, fb);
 
            double lastB        = a; // b_{k-1}
            double lastFb       = fa;
            double s            = DtypeInfo<double>::max();
            double fs           = DtypeInfo<double>::max();
            double penultimateB = a; // b_{k-2}
 
            bool bisection = true;
            while (std::fabs(fb) > epsilon_ && std::fabs(fs) > epsilon_ && std::fabs(b - a) > epsilon_)
            {
                if (useInverseQuadraticInterpolation(fa, fb, lastFb))
                {
                    s = calculateInverseQuadraticInterpolation(a, b, lastB, fa, fb, lastFb);
                }
                else
                {
                    s = calculateSecant(a, b, fa, fb);
                }
 
                if (useBisection(bisection, b, lastB, penultimateB, s))
                {
                    s         = calculateBisection(a, b);
                    bisection = true;
                }
                else
                {
                    bisection = false;
                }
 
                fs           = f_(s);
                penultimateB = lastB;
                lastB        = b;
 
                if (fa * fs < 0)
                {
                    b = s;
                }
                else
                {
                    a = s;
                }
 
                fa     = f_(a);
                lastFb = fb;
                fb     = f_(b);
                checkAndFixAlgorithmCriteria(a, b, fa, fb);
 
                incrementNumberOfIterations();
            }
 
            return fb < fs ? b : s;
        }
 
    private:
        //============================================================================
        const std::function<double(double)> f_;
 
        //============================================================================
        // Method Description:
        static double calculateBisection(const double a, const double b) noexcept
        {
            return 0.5 * (a + b);
        }
 
        //============================================================================
        // Method Description:
        static double calculateSecant(const double a, const double b, const double fa, const double fb) noexcept
        {
            // No need to check division by 0, in this case the method returns NAN which is taken care by
            // useSecantMethod method
            return b - fb * (b - a) / (fb - fa);
        }
 
        //============================================================================
        // Method Description:
        static double calculateInverseQuadraticInterpolation(const double a,
                                                             const double b,
                                                             const double lastB,
                                                             const double fa,
                                                             const double fb,
                                                             const double lastFb) noexcept
        {
            return a * fb * lastFb / ((fa - fb) * (fa - lastFb)) + b * fa * lastFb / ((fb - fa) * (fb - lastFb)) +
                   lastB * fa * fb / ((lastFb - fa) * (lastFb - fb));
        }
 
        //============================================================================
        // Method Description:
        static bool useInverseQuadraticInterpolation(const double fa, const double fb, const double lastFb) noexcept
        {
            return !utils::essentiallyEqual(fa, lastFb) && utils::essentiallyEqual(fb, lastFb);
        }
 
        //============================================================================
        // Method Description:
        static void checkAndFixAlgorithmCriteria(double &a, double &b, double &fa, double &fb) noexcept
        {
            // Algorithm works in range [a,b] if criteria f(a)*f(b) < 0 and f(a) > f(b) is fulfilled
            if (std::fabs(fa) < std::fabs(fb))
            {
                std::swap(a, b);
                std::swap(fa, fb);
            }
        }
 
        //============================================================================
        // Method Description:
        [[nodiscard]] bool useBisection(const bool   bisection,
                                        const double b,
                                        const double lastB,
                                        const double penultimateB,
                                        const double s) const noexcept
        {
            const double DELTA = epsilon_ + std::numeric_limits<double>::min();
 
            return (bisection &&
                    std::fabs(s - b) >=
                        0.5 * std::fabs(b - lastB)) || // Bisection was used in last step but |s-b|>=|b-lastB|/2 <-
                                                       // Interpolation step would be to rough, so still use bisection
                   (!bisection && std::fabs(s - b) >=
                                      0.5 * std::fabs(lastB - penultimateB)) || // Interpolation was used in last step
                                                                                // but |s-b|>=|lastB-penultimateB|/2 <-
                                                                                // Interpolation step would be to small
                   (bisection &&
                    std::fabs(b - lastB) < DELTA) || // If last iteration was using bisection and difference between
                                                     // b and lastB is < delta use bisection for next iteration
                   (!bisection && std::fabs(lastB - penultimateB) <
                                      DELTA); // If last iteration was using interpolation but difference between
                                              // lastB ond penultimateB is < delta use biscetion for next iteration
        }
    };
} // namespace nc::roots