/*=======================================================================
  ch11_bisection.cpp

  Find one root of f(x) = 0 by bisection on an initial interval.

  Origin:
    C++ translation of ch11_bisection.f90 by Alexander Godunov.
    The bisection steps and convergence criterion match that program.

  Example:
    f(x) = sqrt(x) - cos(x) on [0,1].
    Reference root = 0.64171437087288... .

  The program reports the iteration count, root estimate, and f(root).
  The requested tolerance applies to the width of the final bracket.
  A small bracket alone does not establish continuity or exclude a
  singularity, so inspect the residual and the function's domain.
=======================================================================*/

#include <cmath>
#include <iomanip>
#include <iostream>

void bisection(double (*fun)(double), double x1, double x2,
               double eps, double& root, int& flag);
double f(double x);

int main()
{
    // Initial bracket and requested maximum bracket width
    double x1 = 0.0;
    double x2 = 1.0;
    double eps = 1.0e-6;
    double root = 0.0;
    int flag = 0;

    std::cout << "A single root of f(x)\n"
              << "Method - Bisection\n";

    bisection(f, x1, x2, eps, root, flag);

    // Report a missing bracket or invalid input before printing a root
    if (flag == 0) {
        std::cout << "invalid interval or no initial sign change\n";
        return 1;
    }

    std::cout << "iterations  = " << std::abs(flag) << '\n'
              << "tolerance   = " << std::scientific
              << std::setprecision(5) << eps << '\n';
    if (flag < 0)
        std::cout << "iteration limit reached before the requested width\n";

    std::cout << "root        = " << std::scientific
              << std::setprecision(5) << root << '\n'
              << "f(root)     = " << f(root) << '\n';
    return 0;
}

void bisection(double (*fun)(double), double x1, double x2,
               double eps, double& root, int& flag)
/*=======================================================================
  Locate one root of fun(x) on [x1,x2] by repeatedly halving the bracket.

  Input:
    fun   - function to evaluate, continuous on the initial bracket
    x1    - left endpoint of the initial bracket
    x2    - right endpoint of the initial bracket
    eps   - requested maximum width of the final bracket (positive)

  Output:
    root  - midpoint of the final bracket when flag is nonzero
    flag  - convergence indicator and iteration count
            flag > 0: width <= eps; flag is the iteration count
            flag = 0: invalid interval/tolerance or no sign change
            flag < 0: iteration limit reached; abs(flag) is the count

  Notes:
    A sign change guarantees a root only when fun is continuous.
    The maximum number of iterations is 200. If the requested width
    is not reached by then, the estimate is not marked converged.
    The algorithm uses bracket width, not |fun(root)|, to stop.
=======================================================================*/
{
    const int iter = 200;

    // Check that the inputs describe a usable initial bracket
    if (x1 >= x2 || eps <= 0.0) {
        flag = 0;
        return;
    }
    if (fun(x1) * fun(x2) > 0.0) {
        flag = 0;
        return;
    }

    // Initialize the current bracket
    double a = x1;
    double b = x2;

    // Repeatedly retain the half whose endpoints have opposite signs
    for (int i = 1; i <= iter; ++i) {
        // This midpoint form avoids overflow in a+b when both are large.
        double c = a + (b - a) / 2.0;
        if (fun(a) * fun(c) <= 0.0)
            b = c;
        else
            a = c;

        // Stop when the current bracket is sufficiently narrow.
        if (std::abs(b - a) <= eps) {
            root = (b + a) / 2.0;
            flag = i;
            return;
        }
    }

    root = (b + a) / 2.0;
    flag = -iter;
}

// Example function. The chosen bracket must stay in the domain x >= 0.
double f(double x)
{
    // To try the alternative from the Fortran file, use [-1,0].
    // return x + std::cos(x);
    return std::sqrt(x) - std::cos(x);
}