/*=======================================================================
  ch11_multiple_roots.cpp

  Find several distinct roots by scanning an interval and using
  bisection on each subinterval with a sign change.

  Origin:
    C++ translation of ch11_multiple_roots.f90 by Alexander Godunov.
    Based on the default scan plus bisection method in BForce.f90.

  Example:
    f(x) = exp(x) - x*sin(pi*x/2) on [-10.5,2].
    The chapter example has six distinct roots in this interval.

  The tolerance applies to the final width of each root bracket.
  A sign-change scan can miss even-multiplicity roots, pairs of roots
  within one subinterval, and roots near a discontinuity. Inspect a
  plot and repeat with more subintervals when counting roots matters.
=======================================================================*/

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

void MultipleRoots(double (*fun)(double), double x1, double x2, int n,
                   double eps, double roots[], int& nroots, int& nfail);
void BisectBracket(double (*fun)(double), double x1, double x2,
                   double eps, double& root, int& flag);
double f(double x);

int main()
{
    const int n = 1000;
    double roots[n + 1];

    // Search interval, number of scan subintervals, and bracket-width tolerance
    double x1 = -10.5;
    double x2 = 2.0;
    double eps = 1.0e-7;
    int nroots = 0;
    int nfail = 0;

    std::cout << "Roots of f(x) by interval scanning and bisection\n"
              << "interval             = [" << std::fixed
              << std::setprecision(3) << x1 << ", " << x2 << "]\n"
              << "scan subintervals    = " << n << '\n'
              << "bracket-width tol.   = " << std::scientific
              << std::setprecision(5) << eps << '\n';

    MultipleRoots(f, x1, x2, n, eps, roots, nroots, nfail);

    if (nroots == 0) {
        std::cout << "no roots detected by the scan\n";
    } else {
        std::cout << "number       root                  f(root)\n";
        for (int i = 0; i < nroots; ++i) {
            std::cout << std::setw(5) << i + 1 << "   " << std::fixed
                      << std::setprecision(10) << std::setw(18) << roots[i]
                      << "   " << std::scientific << std::setprecision(8)
                      << std::setw(16) << f(roots[i]) << '\n';
        }
    }
    if (nfail > 0)
        std::cout << "brackets not refined to the requested width: "
                  << nfail << '\n';
    return 0;
}

void MultipleRoots(double (*fun)(double), double x1, double x2, int n,
                   double eps, double roots[], int& nroots, int& nfail)
/*=======================================================================
  Scan [x1,x2] and bisect subintervals with opposite endpoint signs.

  Input:
    fun     - function continuous on the search interval
    x1,x2   - endpoints of the search interval, with x1 < x2
    n       - number of equal scan subintervals (positive)
    eps     - maximum final bracket width (positive)

  Output:
    roots   - detected roots, in increasing order; capacity n+1
    nroots  - number of roots stored
    nfail   - number of sign-change brackets that reached the limit

  Notes:
    A zero at a grid point is recorded once, including at x2.
    A sign change only certifies a root when fun is continuous there.
    The scan does not certify that all roots have been found.
=======================================================================*/
{
    nroots = 0;
    nfail = 0;
    if (n <= 0 || x1 >= x2 || eps <= 0.0) return;

    double dx = (x2 - x1) / n;
    double a = x1;
    double fa = fun(a);

    // Each grid point is checked once; a sign change brackets another root.
    for (int j = 1; j <= n; ++j) {
        double b = x1 + j * dx;
        double fb = fun(b);

        if (fa == 0.0) {
            roots[nroots] = a;
            ++nroots;
        } else if (fa * fb < 0.0) {
            double root = 0.0;
            int flag = 0;
            BisectBracket(fun, a, b, eps, root, flag);
            if (flag > 0) {
                roots[nroots] = root;
                ++nroots;
            } else {
                ++nfail;
            }
        }

        a = b;
        fa = fb;
    }

    // The rightmost grid point is not a left endpoint in the loop.
    if (fa == 0.0) {
        roots[nroots] = a;
        ++nroots;
    }
}

void BisectBracket(double (*fun)(double), double x1, double x2,
                   double eps, double& root, int& flag)
/*=======================================================================
  Refine one sign-change bracket by bisection.

  Input:
    fun     - function continuous on [x1,x2]
    x1,x2   - endpoints with opposite signs of fun
    eps     - requested maximum final bracket width

  Output:
    root    - midpoint of the final bracket
    flag    - positive iteration count on convergence; -200 on limit

  Notes:
    The stopping rule is bracket width, as in ch11_bisection.cpp.
=======================================================================*/
{
    const int iter = 200;
    double a = x1;
    double b = x2;
    double fa = fun(a);
    flag = -iter;

    for (int i = 1; i <= iter; ++i) {
        double c = a + (b - a) / 2.0;
        double fc = fun(c);
        if (fa * fc <= 0.0) {
            b = c;
        } else {
            a = c;
            fa = fc;
        }
        if (std::abs(b - a) <= eps) {
            flag = i;
            break;
        }
    }

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

double f(double x)
/*-----------------------------------------------------------------------
  Example from Chapter 11. Replace this expression for another problem.
  Choose a scan interval on which the new function is continuous.
-----------------------------------------------------------------------*/
{
    const double pi = 3.141592653589793;
    return std::exp(x) - x * std::sin(pi * x / 2.0);
}