/*=======================================================================
  ch11_secant.cpp

  Find one root of f(x) = 0 by the secant method.

  Origin:
    C++ translation of ch11_secant.f90 by Alexander Godunov.
    The secant update and residual stopping test match that program.

  Example:
    f(x) = x + cos(x), starting from x1 = -1 and x2 = -0.9.
    Reference root = -0.73908513321516... .

  The program reports the iteration count, root estimate, and f(root).
  The secant method does not preserve a bracket. Starting values can
  affect convergence and may lead to a different root or failure.
=======================================================================*/

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

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

int main()
{
    // Two initial guesses and requested residual tolerance
    double x1 = -1.0;
    double x2 = -0.9;
    double eps = 1.0e-7;
    double root = 0.0;
    int flag = 0;

    std::cout << "A single root of f(x)\n"
              << "Method - Secant\n"
              << "Initial guesses = " << std::fixed
              << std::setprecision(3) << x1 << ", " << x2 << '\n';

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

    // Report invalid input or a zero secant slope before printing a root
    if (flag == 0) {
        std::cout << "invalid tolerance/guesses or zero secant slope\n";
        return 1;
    }

    std::cout << "residual tolerance = " << std::scientific
              << std::setprecision(5) << eps << '\n'
              << "iterations         = " << std::abs(flag) << '\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 Secant(double (*fun)(double), double x1, double x2,
            double eps, double& root, int& flag)
/*=======================================================================
  Locate one root using secant lines through successive iterates.

  Input:
    fun   - function evaluating f(x)
    x1    - first initial guess
    x2    - second, distinct initial guess
    eps   - requested maximum residual |fun(root)| (positive)

  Output:
    root  - most recent secant estimate when flag is nonzero
    flag  - convergence indicator and iteration count
            flag > 0: residual <= eps; flag is the iteration count
            flag = 0: invalid input or zero secant slope
            flag < 0: iteration limit reached; abs(flag) is the count

  Notes:
    This is the ordinary secant method without damping or bracketing.
    The maximum number of iterations is 200. An initial guess that
    already satisfies the residual test is returned with flag=1 to
    follow the positive-success convention, without an update.
=======================================================================*/
{
    const int iter = 200;
    flag = 0;
    if (eps <= 0.0 || x1 == x2) return;

    double a = x1;
    double b = x2;
    double fa = fun(a);
    double fb = fun(b);

    // An accurate starting guess needs no secant update.
    if (std::abs(fa) <= eps) {
        root = a;
        flag = 1;
        return;
    }
    if (std::abs(fb) <= eps) {
        root = b;
        flag = 1;
        return;
    }

    // Form the secant slope from the two most recent function values.
    double c = b;
    for (int i = 1; i <= iter; ++i) {
        double df = (fb - fa) / (b - a);
        if (df == 0.0) return;
        c = b - fb / df;
        double fc = fun(c);

        // Stop on the residual, as in the original Fortran program.
        if (std::abs(fc) <= eps) {
            root = c;
            flag = i;
            return;
        }

        a = b;
        fa = fb;
        b = c;
        fb = fc;
        if (a == b) return;
    }

    root = c;
    flag = -iter;
}

double f(double x)
/*-----------------------------------------------------------------------
  Example function. To use another problem, change this expression.
-----------------------------------------------------------------------*/
{
    return x + std::cos(x);
}