/*
  A single root of an equation f(x)=0 in [a,b] interval
  Methods: Bisectional, False Position, Secant, Newton's
*/
#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;

double root_bs(double(*)(double), double, double, double, int&);
double fzero(double);
double root_fp(double(*)(double), double, double, double, int&);
double fzero(double);
double root_sm(double(*)(double), double, double, double, int&);
double fzero(double);
double root_mr(double(*)(double), double, double, double, double[], int, int&);
double fzero(double);
double root_nm(void(*)(double, double&, double&), double, double, int&);
void   fzeron(double, double&, double&);


int main()
{
    double a,b,root,eps,x1,x2,x;
    double fx, fxp;
    int flag;
    cout.setf(ios::fixed | ios::showpoint);
    cout.precision(5);

    a =   -10.0;                  // left endpoint
    b =    10.0;                  // right endpoint
    eps = 1.0e-6;                 // desired uncertainity of the root
    
    cout << "*******      a = " << setw(7) << a << setw(7) << " b = " << b << endl;

    cout << "                   iterations" << "      root"<<"        f(root)" << endl;
/* Bisection methods */
    root = root_bs(fzero, a ,b ,eps, flag);
    cout << "Bisection method   ";
    if (flag == 0 ) cout << " no root " << endl;
    else
    { cout << setw(7) << flag << setw(14) << root << setw(14) << fzero(root) << endl;}

/* False position method */
    root = root_fp(fzero, a ,b ,eps, flag);
    cout << "False position     ";   
     if (flag == 0 ) cout << " no root " << endl;
    else
    {cout << setw(7) << flag << setw(14) << root<< setw(14) << fzero(root) << endl;}    

/* Secant method */
    x1 = (b+a)/2.0;               // starting point
    x2 = a - 0.01;                // second point
    root = root_sm(fzero, x1, x2, eps, flag);
    cout << "Secant method      ";   
     if (flag == 0 ) cout << " no root for 1000 iterations " << endl;
    else
    {cout << setw(7) << flag << setw(14) << root<< setw(14) << fzero(root) << endl;}

/* Newton method */
    x =   (b+a)/2.0;                       // starting poin
    root = root_nm(fzeron, x, eps, flag);
    cout << "Newton method      ";   
     if (flag == 0 ) cout << " no root, try new x initial  " << endl;
    else
    {   fzeron(root,fx,fxp);
        cout << setw(7) << flag << setw(14) << root<< setw(14) << fx << endl;}


/* Multiple roots */
	const int intervals = 1000;         // number of subintervals to to look for roots
    double roots[intervals];
    int i, nroots;
    root_mr(fzero, a ,b ,eps, roots, intervals ,nroots);
    cout << endl << "Brute force method " << endl;
    cout << " root" << "      value"<<"        f(root)" << endl;    
    if (nroots == 0 ) cout << " no roots found" << endl;
    for (i = 1; i <= nroots; i = i+1)
    cout << setw(5) << i << setw(12) << roots[i-1]<<setw(14) << fzero(root) << endl;
    
    return 0;
}
/*
 *  Function f(x) for the equation f(x)=0
*/
    double fzero(double x)
{
    double y;

    y = x+cos(x);
//    y = exp(x) - x*sin(3.1415926*x/2.0);
//    y = x*x*x - 2.0*x*x -2.0*x +1;
//    y = x*x*x*x -6.0*x*x*x +12.0*x*x -36.0*x-18.0; 
// physics
//   y = (10000.0*1000.0 + 1703.3/(x*x))*(x-0.001169) - 461.5*800.0; 

    return y;
}

/*
 *  Function fzeron(x) Newton method. (includes the first derivative
*/
    void fzeron(double x, double& fx, double& fpx)
{
    fx = x+cos(x);
    fpx = 1.0 - sin(x);
//    fx = exp(x) - x*sin(3.1415926*x/2.0);
//    fpx= exp(x) - sin(3.1415926*x/2.0) - 0.5*x*cos(3.1415926*x/2.0)*3.1415926;
//    fx = x*x*x - 2.0*x*x -2.0*x +1;
//    fpx= 3.0*x*x-4.0*x-2.0;
//    fx = x*x*x*x -6.0*x*x*x +12.0*x*x -36.0*x-18.0;
//    fpx= 4.0*x*x*x - 18.0*x*x +24.0*x -36.0;

}

    double root_bs(double(*f)(double), double a, double b, double eps, int& flag)
/*
====================================================================
 Program to find a single root of an equation f(x)=0
 using the Bisectional method
 written by:    Alex Godunov
 Last revision: February 2022
--------------------------------------------------------------------
 input ...
 f   - function which evaluates f(x) for any x in the interval a,b
 a   - left endpoint of initial interval
 b   - right endpoint of initial interval
 eps - desired uncertainity of the root

 output ...
 x0    - root of equation f(x)=0
 flag - indicator of success
        0 - no solutions for the Bisectional method
        i - a single root found after i iterations

 Limitations: a function f(x) must change sign between a and b
 Max number of allowed iterations is 1000 (accuracy (b-a)/2**1000)
====================================================================
*/
{
    double xl,x0,xr;
    int i, iter=1000;           /* Max number of iterations */

/* check the bisection condition */
    if(f(a)*f(b) > 0.0 )
    { flag = 0;
      return 0.0;}

/* finding a single root */
    i = 0;
    xl = a;
    xr = b;

    while (fabs(xr - xl) >= eps)
    {
     i = i + 1;
     x0 = (xr + xl)/2.0;
     if((f(xl) * f(x0)) <= 0.0 )
        xr = x0;
        else
        xl = x0;
     if(i >= iter) break;
     }
    flag = i;
    return x0;
}


    double root_fp(double(*f)(double), double a, double b, double eps, int& flag)
/*
====================================================================
 Program to find a single root of an equation f(x)=0
 using the Bisectional method
 written by:    Alex Godunov
 Last revision: February 2022
--------------------------------------------------------------------
 input ...
 f   - function which evaluates f(x) for any x in the interval a,b
 a   - left endpoint of initial interval
 b   - right endpoint of initial interval
 eps - desired uncertainity of the root

 output ...
 x0    - root of equation f(x)=0
 flag - indicator of success
        0 - no solutions for the Bisectional method
        i - a single root found after i iterations

 Limitations: a function f(x) must change sign between a and b
 Max number of allowed iterations is 1000 (accuracy (b-a)/2**1000)
====================================================================
*/
{
    double xl,x0,xr;
    int i, iter=100;           /* Max number of iterations */
    cout.setf(ios::fixed | ios::showpoint);
    cout.precision(5);

/* check the bisection condition */
    if(f(a)*f(b) > 0.0 )
    { flag = 0;
      return 0.0;}

/* finding a single root */
    i = 0;
    xl = a;
    xr = b;
    x0 = (xr + xl)/2.0;
//    cout << "  i    a       f(a)      b       f(b)      c       f(c)  "<< endl;
    while (fabs(f(x0)) >= eps)
    {
     i = i + 1;
//     x0 = (xr + xl)/2.0;
     x0 = xr - f(xr)*(xr - xl)/(f(xr)-f(xl));
//     cout << setw(3) << i << setw(9) << xl << setw(9) << f(xl)  
//                          << setw(9) << xr << setw(9) << f(xr) 
//                          << setw(9) << x0 << setw(9) << f(x0) << endl;
     if((f(xl) * f(x0)) <= 0.0 )
        xr = x0;
        else
        xl = x0;
     if(i >= iter) break;
     }
    flag = i;
    return x0;
}
    double root_sm(double(*f)(double), double x1, double x2, double eps, int& flag)
/*
====================================================================
 Program to find a single root of an equation f(x)=0
 using the method of secants
 written by:    Alex Godunov
 last revision: February 2022
--------------------------------------------------------------------
 input ...
 f  - function which evaluates f(x) for any x in the interval a,b
 x1 - initial point
 x2 - second initial point
 eps    - desired uncertainity of the root

 output ...
 x3    - root of equation f(x)=0
 iflag - indicator of success
         0 - no solutions found for the secant method
         i - a single root found after i iterations

 Max number of allowed iterations is 1000
====================================================================
*/
{
    double x3;
    int i, iter=1000;
    i = 0;
    while (fabs(x2 - x1) >= eps)
    {
     i = i + 1;
     x3 = x2 - (f(x2)*(x2-x1))/(f(x2)-f(x1));
     x1 = x2;
     x2 = x3;
     if(i >= iter) break;
     }

    flag = i;
    if (i == iter) flag = 0;
    return x3;
}

    double root_mr(double(*f)(double), double a, double b, double eps, 
    double roots[], int intervals, int& nroots)
/*
====================================================================
 Program to find multiple roots of an equation f(x)=0
 on an interval [a,b] using the Bisectional method for subintervals
 written by:    Alex Godunov
 Last revision: February 2022
--------------------------------------------------------------------
 input ...
 f   - function which evaluates f(x) for any x in the interval a,b
 a   - left endpoint of initial interval
 b   - right endpoint of initial interval
 eps - desired uncertainity of the root
 intervals - number of subintervals for [a,b]

 output ...
 roots  - roots of equation f(x)=0
 nroots - number of roots 

 Limitations: a function f(x) must change sign between a and b
 Max number of allowed iterations is 1000 (accuracy (b-a)/2**1000)
====================================================================
*/        
{
   double xl, x0, xr, dx;    
   int i, j, n, iter=1000;
   n = 0;                   // n is the counter for roots
   dx = (b-a)/intervals; 
//* start loop over subintervals
   j= 0;
   while (j < intervals)
   { j = j + 1;
	xl = a + (j-1)*dx;
    xr = xl + dx;
    if(f(xl)*f(xr) > 0.0) continue;
//** start bisectional loop over the subinterval j
    i = 0;
    while (fabs(xr - xl) >= eps)
    {
	  i = i + 1;
      x0 = (xr + xl)/2.0;
      if((f(xl) * f(x0)) <= 0.0 ) 
        xr = x0;
        else 
        xl = x0;
      if(i >= iter) break;   
    }
    if (i >= iter) continue;  //do not count roots when max. iterations exceeded
    roots[n] = x0;
    n = n+1;
   }
    nroots = n; 
    return 0;
} 
    double root_nm(void(*f)(double, double&, double&), 
                   double x, double eps, int& flag)
/*
====================================================================
 Program to find a single root of an equation f(x)=0
 using Newton's method 
 written by:    Alex Godunov
 last revision: February 2022
--------------------------------------------------------------------
 input ...
 f - function which evaluates f(x) for any x in the interval a,b
 x - initial point
 eps    - desired uncertainity of the root

 output ...
 xc    - root of equation f(x)=0
 iflag - indicator of success
         0 - no solutions found for the secant method
         i - a single root found after i iterations

 Max number of allowed iterations is 1000
====================================================================
*/
{
    double fx, fpx, xc;
    int i, iter=1000;
    i = 0;
    do {
        i = i + 1;
        f(x,fx,fpx);
        xc = x - fx/fpx;
        x = xc;
        if(i >= iter) break;
        } while (fabs(fx) >= eps);    
    
    flag = i;
    if (i == iter) flag = 0;
    return xc;
}