/*=======================================================================
  ch11_false_position.cpp

  Find one root of f(x) = 0 by the false position (regula falsi) method.

  Origin:
    C++ translation of ch11_false_position.f90 by Alexander Godunov.
    The interpolation formula and convergence criterion match that file.

  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 residual |f(root)|. Unlike
  bisection, ordinary false position may leave one endpoint fixed for
  many iterations, so the bracket width is not the stopping test.
=======================================================================*/

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

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

int main()
{
    // Initial bracket and requested residual tolerance
    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 - False position\n";

    false_position(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'
              << "residual tolerance = " << std::scientific
              << std::setprecision(5) << eps << '\n';
    if (flag < 0)
        std::cout << "iteration limit reached before the requested residual\n";

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

void false_position(double (*fun)(double), double x1, double x2,
                    double eps, double& root, int& flag)
/*=======================================================================
  Locate one root of fun(x) on [x1,x2] by ordinary false position.

  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 residual |fun(root)| (positive)

  Output:
    root  - most recent false-position estimate when flag is nonzero
    flag  - convergence indicator and iteration count
            flag > 0: |fun(root)| <= 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. One bracket endpoint
    can remain nearly fixed, so a small interval width is not assumed.
    An endpoint root is returned immediately with a positive flag.
=======================================================================*/
{
    const int iter = 200;
    flag = 0;

    // Check that the inputs describe a usable initial bracket
    if (x1 >= x2 || eps <= 0.0) return;

    double a = x1;
    double b = x2;
    double fa = fun(a);
    double fb = fun(b);
    if (fa * fb > 0.0) return;

    // Handle an exact zero at an endpoint before forming a secant line.
    // Flag 1 follows the existing convention that positive means success.
    if (fa == 0.0) {
        root = a;
        flag = 1;
        return;
    }
    if (fb == 0.0) {
        root = b;
        flag = 1;
        return;
    }

    // Interpolate between the current endpoints and retain the sign change.
    double c = a;
    for (int i = 1; i <= iter; ++i) {
        c = b - fb * (b - a) / (fb - fa);
        double fc = fun(c);

        // False position uses the residual, not the bracket width, to stop.
        if (std::abs(fc) <= eps) {
            root = c;
            flag = i;
            return;
        }

        if (fa * fc <= 0.0) {
            b = c;
            fb = fc;
        } else {
            a = c;
            fa = fc;
        }
    }

    root = c;
    flag = -iter;
}

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