bisection.h

class Bisection : public Iteration {
public:
    Bisection(double epsilon, const std::function<double (double)> &f) : Iteration(epsilon), mf(f) {}

    double solve(double a, double b) override {
        resetNumberOfIterations();

        checkAndFixAlgorithmCriteria(a, b);

        fmt::print("Bisection -> [{:}, {:}]\n", a, b);
        fmt::print("{:<5}|{:<20}|{:<20}|{:<20}|{:<20}\n", "K", "a", "b", "x", "f(x)");
        fmt::print("------------------------------------------------------------------------------------ \n");

        double x = 0.5*(a+b);
        double fx = mf(x);
        fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), a, b, x, fx);

        while(fabs(fx) >= epsilon()) {
            x = calculateX(x, a, b, fx);
            fx = mf(x);

            fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), a, b, x, fx);
        }

        fmt::print("\n");

        return x;
    }

private:
    void checkAndFixAlgorithmCriteria(double &a, double &b) const {
        //Algorithm works in range [a,b] if criteria f(a)*f(b) < 0 and f(a) > f(b) is fulfilled
        assert(mf(a)*mf(b) < 0);
        if(mf(a) < mf(b)) {
            std::swap(a,b);
        }
    }

    static double calculateX(double x, double &a, double &b, double fx) {
        if(fx >= 0) a = x;
        if(fx < 0) b = x;

        return 0.5*(a+b);
    }

    const std::function<double (double)> &mf;
};

brent.h

class Brent : public Iteration {
public:
    Brent(double epsilon, const std::function<double (double)> &f) : Iteration(epsilon), mf(f) {}

    double solve(double a, double b) override {
        resetNumberOfIterations();

        double fa = mf(a);
        double fb = mf(b);

        checkAndFixAlgorithmCriteria(a, b, fa, fb);

        fmt::print("Brent -> [{:}, {:}]\n", a, b);
        fmt::print("{:<5}|{:<20}|{:<20}|{:<20}|{:<20}|{:<20}\n", "K", "a", "b", "f(a)", "f(b)", "f(s)");
        fmt::print("---------------------------------------------------------------------------------------------------------- \n");
        fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), a, b, fa, fb, "");

        double lastB = a; // b_{k-1}
        double lastFb = fa;
        double s = std::numeric_limits<double>::max();
        double fs = std::numeric_limits<double>::max();
        double penultimateB = a; // b_{k-2}

        bool bisection = true;
        while(fabs(fb) > epsilon() && fabs(fs) > epsilon() && 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, a, b, lastB, penultimateB, s)) {
                s = calculateBisection(a, b);
                bisection = true;
            }
            else {
                bisection = false;
            }

            fs = mf(s);
            penultimateB = lastB;
            lastB = b;

            if(fa*fs < 0) {
                b = s;
            }
            else {
                a = s;
            }

            fa = mf(a);
            lastFb = fb;
            fb = mf(b);
            checkAndFixAlgorithmCriteria(a, b, fa, fb);

            fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), a, b, fa, fb, fs);
        }

        fmt::print("\n");

        return fb < fs ? b : s;
    }

private:
    static double calculateBisection(double a, double b) {
        return 0.5*(a+b);
    }

    static double calculateSecant(double a, double b, double fa, double fb) {
        //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);
    }

    static double calculateInverseQuadraticInterpolation(double a, double b, double lastB, double fa, double fb, double lastFb) {
        return a*fb*lastFb/((fa-fb)*(fa-lastFb)) +
               b*fa*lastFb/((fb-fa)*(fb-lastFb)) +
               lastB*fa*fb/((lastFb-fa)*(lastFb-fb));
    }

    static bool useInverseQuadraticInterpolation(double fa, double fb, double lastFb) {
        return fa != lastFb && fb != lastFb;
    }

    static void checkAndFixAlgorithmCriteria(double &a, double &b, double &fa, double &fb) {
        //Algorithm works in range [a,b] if criteria f(a)*f(b) < 0 and f(a) > f(b) is fulfilled
        assert(fa*fb < 0);
        if (fabs(fa) < fabs(fb)) {
            std::swap(a, b);
            std::swap(fa, fb);
        }
    }

    bool useBisection(bool bisection, double a, double b, double lastB, double penultimateB, double s) const {
        static const double DELTA = epsilon() + std::numeric_limits<double>::min();

        return (bisection && fabs(s-b) >= 0.5*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 && fabs(s-b) >= 0.5*fabs(lastB-penultimateB)) ||     //Interpolation was used in last step but |s-b|>=|lastB-penultimateB|/2 <- Interpolation step would be to small
               (bisection  && fabs(b-lastB) < DELTA) ||                         //If last iteration was using bisection and difference between b and lastB is < delta use bisection for next iteration
               (!bisection && fabs(lastB-penultimateB) < DELTA);                //If last iteration was using interpolation but difference between lastB ond penultimateB is < delta use biscetion for next iteration
    }

    const std::function<double (double)> &mf;
};

dekker.h

class Dekker : public Iteration {
public:
    Dekker(double epsilon, const std::function<double (double)> &f) : Iteration(epsilon), mf(f) {}

    double solve(double a, double b) override {
        resetNumberOfIterations();

        double fa = mf(a);
        double fb = mf(b);

        checkAndFixAlgorithmCriteria(a, b, fa, fb);

        fmt::print("Dekker -> [{:}, {:}]\n", a, b);
        fmt::print("{:<5}|{:<20}|{:<20}|{:<20}|{:<20}\n", "K", "a", "b", "f(a)", "f(b)");
        fmt::print("------------------------------------------------------------------------------------ \n");
        fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), a, b, fa, fb);

        double lastB = a;
        double lastFb = fa;

        while(fabs(fb) > epsilon() && fabs(b-a) > epsilon()) {
            const double s = calculateSecant(b, fb, lastB, lastFb);
            const double m = calculateBisection(a, b);

            lastB = b;

            b = useSecantMethod(b, s, m) ? s : m;

            lastFb = fb;
            fb = mf(b);

            if (fa * fb > 0 && fb * lastFb < 0) {
                a = lastB;
            }

            fa = mf(a);
            checkAndFixAlgorithmCriteria(a, b, fa, fb);

            fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), a, b, fa, fb);
        }

        fmt::print("\n");

        return b;
    }

private:
    static void checkAndFixAlgorithmCriteria(double &a, double &b, double &fa, double &fb) {
        //Algorithm works in range [a,b] if criteria f(a)*f(b) < 0 and f(a) > f(b) is fulfilled
        assert(fa*fb < 0);
        if (fabs(fa) < fabs(fb)) {
            std::swap(a, b);
            std::swap(fa, fb);
        }
    }

    static double calculateSecant(double b, double fb, double lastB, double lastFb) {
        //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-lastB)/(fb-lastFb);
    }

    static double calculateBisection(double a, double b) {
        return 0.5*(a+b);
    }

    static bool useSecantMethod(double b, double s, double m) {
        //Value s calculated by secant method has to be between m and b
        return (b > m && s > m && s < b) ||
               (b < m && s > b && s < m);
    }

    const std::function<double (double)> &mf;
};

newton.h

class Newton : public Iteration {
public:
    Newton(double epsilon, const std::function<double (double)> &f, const std::function<double (double)> &fPrime) : Iteration(epsilon), mf(f), mfPrime(fPrime) {}

    double solve(double x) override {
        resetNumberOfIterations();

        fmt::print("Newton -> [x0 = {:}]\n", x);
        fmt::print("{:<5}|{:<20}|{:<20}|{:<20}\n", "K", "x", "f(x)", "f'(x)");
        fmt::print("------------------------------------------------------------------ \n");

        double fx = mf(x);
        double fxPrime = mfPrime(x);
        fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), x, fx, fxPrime);

        while(fabs(fx) >= epsilon()) {
            x = calculateX(x, fx, fxPrime);

            fx = mf(x);
            fxPrime = mfPrime(x);

            fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), x, fx, fxPrime);
        }

        fmt::print("\n");

        return x;
    }

private:
    static double calculateX(double x, double fx, double fxPrime) {
        assert(fabs(fxPrime) >= std::numeric_limits<double>::min());

        return x - fx/fxPrime;
    }

    const std::function<double (double)> &mf;
    const std::function<double (double)> &mfPrime;
};

secant.h

class Secant : public Iteration {
public:
    Secant(double epsilon, const std::function<double (double)> &f) : Iteration(epsilon), mf(f) {}

    double solve(double a, double b) override {
        resetNumberOfIterations();

        fmt::print("Secant -> [{:}, {:}]\n", a, b);
        fmt::print("{:<5}|{:<20}|{:<20}|{:<20}|{:<20}\n", "K", "a", "b", "f(a)", "f(b)");
        fmt::print("--------------------------------------------------------------------------------------\n");

        if(mf(a) > mf(b)) {
            std::swap(a,b);
        }

        double x = b;
        double lastX = a;
        double fx = mf(b);
        double lastFx = mf(a);

        fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), x, lastX, fx, lastFx);

        while(fabs(fx) >= epsilon()) {
            const double x_tmp = calculateX(x, lastX, fx, lastFx);

            lastFx = fx;
            lastX = x;
            x = x_tmp;

            fx = mf(x);

            fmt::print("{:<5}|{:<20.15f}|{:<20.15f}|{:<20.15f}|{:<20.15f}\n", incrementNumberOfIterations(), x, lastX, fx, lastFx);
        }

        fmt::print("\n");

        return x;
    }

private:
    static double calculateX(double x, double lastX, double fx, double lastFx) {
        const double functionDifference = fx - lastFx;
        assert(fabs(functionDifference) >= std::numeric_limits<double>::min());

        return x - fx*(x-lastX)/functionDifference;
    }

    const std::function<double (double)> &mf;
};