/*********************************************************************
   Finds the root of a function, f(x)=0, using Brent's method. The 
   root, returned as zbrent, will be refined until is accuracy is tol.
   C.A. Bertulani        May/25/2000
*********************************************************************/
typedef double Number;

#include <math.h>
#include<iostream.h>
#include "butil.h"

int main(){
	Number func(Number );
	int zbrac(Number (*)(Number), Number *, Number *);
	Number zbrent(Number (*)(Number), Number, Number, Number );
	Number a, b, tol, res, eps, x;

	tol = 1.e-6;
	eps = 0.001;
	/* Find minimum of sin(x) between pi and 2*pi */	 
	cout << "\nI will find the zero of sin(x) somewhere around a\n";
	cout <<   "-------------------------------------------------\n\n";
	cout << "Enter a\n";
	cin >> a;
	b = a + eps;
	int c = zbrac(&func,&a,&b);  /* First bracket the interval. */
	if(c == 1){
		cout <<"\nCongratulations!. Bracketing successful";
		cout << "\nHere is a good bracketed interval to start search\n";
		cout << a << "   " << b << endl;
		cout << "\nNow let's get that root.";
	}
	else{
		cout <<"\nIf you've got this message, don't waste your time playing the lottery";
		cout <<"Sorry. Bracketing failed!";
	}
	x = zbrent(&func,a,b,tol);  /* Find the root and where it is. */
	cout << "\n\nAnd the winner of the root search is...   x = " <<x;
	res = sin(x);
	cout << "\nThe function there is ... sin(x) = " << res;
	cout << "\n\nThis is too easy! Test another function.\n";
	return 0;
}
/*********************************************************************
	Function to find minimum.
*********************************************************************/
Number func(Number x){
		return sin(x);
}
/*********************************************************************
   Finds the root of a function, f(x)=0, in an interval x1, x2, using 
   Brent's method. The root, returned as zbrent, will be refined until
   is accuracy is tol.
   This is a robust program, which should be used for general purpose.
   C.A. Bertulani        May/25/2000
*********************************************************************/
#define ITMAX 100			/* Maximum allowed number of iterations */
#define EPS 3.0e-8			/* Machine Numbering point precision. 
								You can change it for double precision.*/

Number zbrent(Number (*func)(Number), Number x1, Number x2, Number tol)
{
	int iter;
	Number a=x1,b=x2,c=x2,d,e,min1,min2;
	Number fa=(*func)(a),fb=(*func)(b),fc,p,q,r,s,tol1,xm;

	if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0))
		cerr << "Root must be bracketed in zbrent";
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++) {
		if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
			c=a;		/* rename a, b, c, and adjust bounding interval d. */
			fc=fa;
			e=d=b-a;
		}
		if (fabs(fc) < fabs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*fabs(b)+0.5*tol;		/* Convergence check. */
		xm=0.5*(c-b);
		if (fabs(xm) <= tol1 || fb == 0.0) return b;
		if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
			s=fb/fa;		/* Attempt inverse quadratic interpolation. */
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0) q = -q;		/* Check whether in bounds. */
			p=fabs(p);
			min1=3.0*xm*q-fabs(tol1*q);
			min2=fabs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;		/* Accept interpolation. */
				d=p/q;
			} else {
				d=xm;		/* Interpolation failed, use bisection. */
				e=d;
			}
		} else {		/* Bounds decreasing too slowly use bisection. */
			d=xm;
			e=d;
		}
		a=b;		/* Move last guess to a. */
		fa=fb;
		if (fabs(d) > tol1)		/* Evaluate new trial root. */
			b += d;
		else
			b += SIGN(tol1,xm);
		fb=(*func)(b);
	}
	cerr << "Maximum number of iterations exceeded in zbrent\n";
	return 0.0;		/* Never get here. */
}
#undef ITMAX
#undef EPS
/*********************************************************************
   Given a function func and its initial guessed range x1 to x2, the 
   routine expands the range geometrically until a root is bracketed
   by the returned values x1 and x2 (in which case zbrac returns 1)
   or until the range becomes unacceptably large (in which case zbrac
   returns 0).
   C.A. Bertulani        May/25/2000
*********************************************************************/
#define FACTOR 1.6
#define NTRY 50

int zbrac(Number (*func)(Number), Number *x1, Number *x2)
{
	int j;
	Number f1,f2;

	if (*x1 == *x2) cerr << "Bad initial range in zbrac";
	f1=(*func)(*x1);
	f2=(*func)(*x2);
	for (j=1;j<=NTRY;j++) {
		if (f1*f2 < 0.0) return 1;
		if (fabs(f1) < fabs(f2))
			f1=(*func)(*x1 += FACTOR*(*x1-*x2));
		else
			f2=(*func)(*x2 += FACTOR*(*x2-*x1));
	}
	return 0;
}
#undef FACTOR
#undef NTRY
/*********************************************************************/