/*=======================================================================
  ch11_newton.cpp

  Find one root of f(x) = 0 by Newton's method.

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

  Example:
    f(x) = x + cos(x), starting from x0 = -1.
    Reference root = -0.73908513321516... .

  The program reports the iteration count, root estimate, and f(root).
  Newton's method is local: a different initial guess may converge to
  another root or fail. A small residual should be interpreted in the
  context of the function and the physical problem.
=======================================================================*/

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

void Newton1(double (*fun)(double, int), double x0, double eps,
             double& root, int& flag);
double f(double x, int k);

int main()
{
    // Initial guess and requested residual tolerance
    double x0 = -1.0;
    double eps = 1.0e-7;
    double root = 0.0;
    int flag = 0;

    std::cout << "A single root of f(x)\n"
              << "Method - Newton\n"
              << "Initial guess = " << std::fixed
              << std::setprecision(3) << x0 << '\n';

    Newton1(f, x0, eps, root, flag);

    // Report invalid input or a zero derivative before printing a root
    if (flag == 0) {
        std::cout << "invalid tolerance or zero derivative at an iterate\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, 0) << '\n';
    return 0;
}

void Newton1(double (*fun)(double, int), double x0, double eps,
             double& root, int& flag)
/*=======================================================================
  Locate one root near x0 using Newton's update x_new=x-f(x)/f'(x).

  Input:
    fun   - function evaluating f(x) for k=0 and f'(x) for k=1
    x0    - initial guess
    eps   - requested maximum residual |fun(root,0)| (positive)

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

  Notes:
    This is the ordinary Newton method without damping or bracketing.
    The maximum number of iterations is 200. If the initial guess
    already satisfies the residual test, flag=1 denotes success even
    though no Newton update was needed.
=======================================================================*/
{
    const int iter = 200;
    flag = 0;
    if (eps <= 0.0) return;
    double x = x0;

    // A sufficiently accurate initial guess needs no Newton update.
    double fx = fun(x, 0);
    if (std::abs(fx) <= eps) {
        root = x;
        flag = 1;
        return;
    }

    // Repeatedly use the value and derivative at the current estimate.
    for (int i = 1; i <= iter; ++i) {
        double dfx = fun(x, 1);
        if (dfx == 0.0) return;

        double xnew = x - fx / dfx;
        fx = fun(xnew, 0);

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

    root = x;
    flag = -iter;
}

double f(double x, int k)
/*-----------------------------------------------------------------------
  Evaluate the example function for k=0 and its derivative for k=1.
  To use another function, change both expressions consistently.
-----------------------------------------------------------------------*/
{
    if (k == 0) return x + std::cos(x);
    if (k == 1) return 1.0 - std::sin(x);
    return 0.0;  // The solver calls only k=0 and k=1.
}