libcruft-util/roots/bisection.hpp

/*
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
 *
 * Copyright 2016 Danny Robson <danny@nerdcruft.net>
 */

#pragma once

#include "../maths.hpp"

#include <cmath>

namespace cruft::roots {
    /// find a root of a function using the bisection method
    ///
    /// the user is responsible for ensuring there is in fact a root. there
    /// is no iteration limit currently, so poorly considered arguments could
    /// result in an infinite loop.
    ///
    /// a: lower bound
    /// b: upper bound
    /// f: real function
    /// tolerance: minimum range termination condition
    template <typename T>
    T
    bisection (T a, T b, T (*f)(T), T tolerance)
    {
        // to simplify the exit conditions we assume that a..b is an
        // increasing range.
        if (a > b)
            std::swap (a, b);

        while (b - a > tolerance) {
            // calculate the midpoint
            auto   c = (a + b) / T{2};
            auto f_c = f (c);

            // return early if we happened across the root
            if (exactly_zero (f_c))
                return c;

            // check which interval we fell into, and divide for the next
            // iteration
            auto f_a = f (a);
            if (samesign (f_a, f_c))
                a = c;
            else
                b = c;
        }

        // we didn't find it exactly, but it's somewhere in this range so take
        // a punt and return the middle of the range. not sure if this is
        // valid mathematically.
        return (a + b) / T{2};
    }
}