#include <iostream>
#include <fstream>
#include <cmath>
#include "testfun.h"
#include <iomanip>

using namespace std;

void bisect(double xa, double xb, int nmax, double &xr, double &errfun)
{
	double fr;
	double fa = FF(xa);
	double fb = FF(xb);
	
	if(fabs(fa) < errfun)
	{
		xr     = xa;
		errfun = fabs(fa);
		return ;
	}
	else if(fabs(fb) < errfun)
	{
		xr     = xb;
		errfun = fabs(fb);
		return ;
	}
	else if(fa*fb > 0.e0)
	{
		cout<<fa<<"  "<<fb<<endl;
		cout<<"error data (BISECT): f(xa)*f(xb) > 0"<<endl;
		return ;
	}
//
//-- Convergence history recorded in file CONV
	ofstream oconv("conv");
	oconv<<xa<<"  "<<fa<<endl;
	oconv<<xb<<"  "<<fb<<endl;
//
//___ start iterating _______________
//
	bool maxIter = false;
	cout.precision(10);
	cout<<fixed<<scientific<<endl;
	for(int niter=1;niter<=nmax;niter++)
	{
		xr=0.5e0*(xa+xb);
		fr=FF(xr);
		oconv<<xr<<"  "<<fr<<endl;
		cout<<setprecision(5)<<niter<<"  "<<setprecision(15)<<xr<<"  "<<setw(5)<<fr<<endl;	//screen message
//-- Is XR a root?
		if(fabs(fr) < errfun)
		{
			errfun=fabs(fr);
			maxIter = false;
			break;
		}
		else
		{ 
//-- define the new interval
			if(fa*fr < 0.e0)  // continue with the LEFT sub-interval
			{
				xb=xr;	
				fb=fr;
			}
			else              // continue with the RIGHT sub-interval
			{
				xa=xr;	
				fa=fr;
			}
		}
		maxIter = true;
	}
	if(maxIter)
        	cout<<"Max. Number of Iterations reached WITHOUT CONVERGENCE\n";

	oconv.close();
}