/*********************************************************************
   Finds the minimum of a function in an interval ax, bx, using the
   Golden search method.
   C.A. Bertulani        May/18/2000
*********************************************************************/
typedef double Number;

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

int main(){
	Number func(Number );
	void mnbrak(Number*, Number*, Number*, Number*, Number*, Number*,Number (*)(Number));
	Number golden(Number, Number, Number, Number (*)(Number), Number, Number *);
	Number brent(Number, Number, Number, Number (*)(Number), Number, Number *);
	Number a, b, c, fa, fb, fc, tol, x, res, pi, eps;

	tol = 1.e-8;
	pi = acos(-1.);
	eps = 0.01;
	/* Find minimum of sin(x) between pi and 2*pi */	 
	cout << "\nI will find the minimum of sin(x) somewhere around a\n";
	cout << "----------------------------------------------------\n\n";
	cout <<"Enter a\n";
	cin >> a;
	b=a+eps;
	mnbrak(&a,&b,&c,&fa,&fb,&fc,&func);  /* First backet the interval. */
	cout << "\nHere the optmized bracketed interval\n";
	cout << a << "   " << b << "  " << c;
//	a=3.1;
//	b=3.3;
//	c=6.2;
	cout << "\nThe above routine sucks for functions with more than one minimum.\n"; 
	cout << "Sometimes results are very bad. Ex: Just try to use a=3.\n";
	cout << "Solution for such functions: bracket it yourself. Golden will do the rest.\n";
	res = golden(a,b,c,&func,tol,&x);  /* Find the minimum and where it is. */
	cout << "\nAnd the winner of the golden search is...   x = " << x;
	cout << "\nThe function there is ... sin(x) = " << res;
	cout << "\n\nNow, use Brent's method: our preferred one\n";
	res = brent(a,b,c,&func,tol,&x);  /* Find the minimum and where it is. */
	cout << "\nAnd the winner of the brent's search is...   x = " << x;
	cout << "\nThe function there is ... sin(x) = " << res << endl;
	cout << "\n\n But, this is too easy! Test another function.\n";
	return 0;
}
/*********************************************************************
	Function to find minimum.
*********************************************************************/
Number func(Number x){
		return sin(x);
}
/*********************************************************************
   Given a function f, and given a bracketing triplet of abcissas ax,
   bx, and cx (such that bx is between ax and cx, and f(bx) is less
   than both f(ax) and f(cx)), this routine performs a golden search
   for the minimum, isolating it to a fractional precision of about tol.
   The abcissa of the minimum,  is returned as xmin, and the minimum
   function value is returned as golden, the returned function value.
   tol should be set to the square root of the machine precision.
   C.A. Bertulani        May/18/2000
*********************************************************************/
#define R 0.61803399  /* The golden ratios.  */
#define C (1.0-R)
#define SHFT2(a,b,c) (a)=(b);(b)=(c);
#define SHFT3(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);

Number golden(Number ax, Number bx, Number cx, Number (*f)(Number), Number tol,
	Number *xmin)
{
	Number f1,f2,x0,x1,x2,x3;

	x0=ax;			/* At any given time we will keep track of four points */
	x3=cx;			/* x0,x1,x2,x3. */
	if (fabs(cx-bx) > fabs(bx-ax)) {		/* Make x0 to x1 the smaller segment, */
		x1=bx;
		x2=bx+C*(cx-bx);				/* and fill in the new point to be tried. */
	} else {
		x2=bx;
		x1=bx-C*(bx-ax);
	}
	f1=(*f)(x1);			/* The initial function evaluations. Note that */
	f2=(*f)(x2);			/* we never need to evaluate the function   */
	while (fabs(x3-x0) > tol*(fabs(x1)+fabs(x2))) {	/* at the original endpoints. */
		if (f2 < f1) {		/* One possible outcome,  */
			SHFT3(x0,x1,x2,R*x1+C*x3)	/* its housekeeping, */
			SHFT2(f1,f2,(*f)(x2))		/* and a new function evaluation. */
		} else {						/* The other outcome,  */
			SHFT3(x3,x2,x1,R*x2+C*x0)
			SHFT2(f2,f1,(*f)(x1))		/* and its new function evaluation.  */
		}
	}
	if (f1 < f2) {				/* Back to see if we are done.  */
		*xmin=x1;				/* We are done. Output the best of the two   */
		return f1;				/* current values.  */
	} else {
		*xmin=x2;
		return f2;
	}
}
#undef C
#undef R
#undef SHFT2
#undef SHFT3
/****************************************************************************
   Given a function func, and given distinct points ax and bx, this 
   routine searches in the downhill direction (defined by the function 
   as evaluated at the initial points) and returns new points ax,bx,cx
   that bracket a minimum of the function. Also returned are the function
   values at the three points, fa, fb, and fc.
   C.A. Bertulani        May/18/2000
*****************************************************************************/
#define GOLD 1.618034	/* Default ratio by which successive intervals are magnified. */
#define GLIMIT 100.0	/* The maximum magnification allowed for a parabolic-fit step. */
#define TINY 1.0e-20
#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);

void mnbrak(Number *ax, Number *bx, Number *cx, Number *fa, Number *fb, Number *fc,
	Number (*func)(Number))
{
	Number ulim,u,r,q,fu,dum;

	*fa=(*func)(*ax);
	*fb=(*func)(*bx);
	if (*fb > *fa) {		/* Switch roles of a and b so that we can go */
		SHFT(dum,*ax,*bx,dum)		/* downhill in the direction from a to b.  */
		SHFT(dum,*fb,*fa,dum)
	}
	*cx=(*bx)+GOLD*(*bx-*ax);		/* First guess for c. */
	*fc=(*func)(*cx);
	while (*fb > *fc) {		/* Keep returning here until we bracket. */
		r=(*bx-*ax)*(*fb-*fc);	/* Compute u by parabolic extrapolation  from */
		q=(*bx-*cx)*(*fb-*fa);		/* a,b,c. TINY is used to prevent any possible */
		u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/	/* division by zero.  */
			(2.0*SIGN(FMAX(fabs(q-r),TINY),q-r));
		ulim=(*bx)+GLIMIT*(*cx-*bx);
		/* we won't go farther than this. Test various possibilities. */
		if ((*bx-u)*(u-*cx) > 0.0) {	/* Parabolic u is between b and c; try it. */
			fu=(*func)(u);
			if (fu < *fc) {	/* Got a minimum between b anc c. */
				*ax=(*bx);
				*bx=u;
				*fa=(*fb);
				*fb=fu;
				return;
			} else if (fu > *fb) {		/* Got a minimum between a and u. */
				*cx=u;
				*fc=fu;
				return;
			}
			u=(*cx)+GOLD*(*cx-*bx);	/* Parabolic fit was no use. Use default magnification. */
			fu=(*func)(u);
		} else if ((*cx-u)*(u-ulim) > 0.0) {	 /* Parabolic fit is between c and its */
			fu=(*func)(u);				/* allowed limit. */
			if (fu < *fc) {
				SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
				SHFT(*fb,*fc,fu,(*func)(u))
			}
		} else if ((u-ulim)*(ulim-*cx) >= 0.0) {	/* Limit parabolic u to maximum allowed */
			u=ulim;								/* value. */
			fu=(*func)(u);
		} else {			/* Reject parabolic u, use default magnification.  */
			u=(*cx)+GOLD*(*cx-*bx);
			fu=(*func)(u);
		}
		SHFT(*ax,*bx,*cx,u)		/* Eliminate oldest point and continue. */
		SHFT(*fa,*fb,*fc,fu)
	}
}
#undef GOLD
#undef GLIMIT
#undef TINY
#undef SHFT
/****************************************************************************
   Given a function f, and given a bracketing triple ax, bx,  and cx (such 
   that bx is between ax and cx, and f(ax) is less than both f(ax) and f(cx)),
   this routine isolates the minimum to a fractional precision of about tol
   using Brent's method. The abscissa of the minimum is returned as xmin, and
   the minimum function value is returned as brent, the returned function value.
   C.A. Bertulani        May/18/2000
*****************************************************************************/
#define ITMAX 100 /* The maximum allowed number of iterations */
#define CGOLD 0.3819660 /* Golden ratio. */
#define ZEPS 1.0e-10	/* Small number that protects against trying to achieve 
							fractional accuracy for a minimum that happens to be
							exactly zero. */

#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);

Number brent(Number ax, Number bx, Number cx, Number (*f)(Number), Number tol,
	Number *xmin)
{
	int iter;
	Number a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
	Number e=0.0;		/* This will be the distance moved on the step
						  before the last.  */

	a=(ax < cx ? ax : cx);		/* a and b must be in ascending order, */
	b=(ax > cx ? ax : cx);		/* but input abcissas need not be. */
	x=w=v=bx;					/* Initializations... */
	fw=fv=fx=(*f)(x);
	for (iter=1;iter<=ITMAX;iter++) {		/* Main program loop. */
		xm=0.5*(a+b);
		tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
		if (fabs(x-xm) <= (tol2-0.5*(b-a))) {		/* Test for done here. */
			*xmin=x;
			return fx;
		}
		if (fabs(e) > tol1) {		/* Construct a trial parabolic fit. */
			r=(x-w)*(fx-fv);
			q=(x-v)*(fx-fw);
			p=(x-v)*q-(x-w)*r;
			q=2.0*(q-r);
			if (q > 0.0) p = -p;
			q=fabs(q);
			etemp=e;
			e=d;
			if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
				d=CGOLD*(e=(x >= xm ? a-x : b-x));
			/* The above conditions determine the acceptability of the parabolic fit.
				Here we take the golden search step into the larger of the two
				segments. */
			else {
				d=p/q;		/* Take the parabolic step. */
				u=x+d;
				if (u-a < tol2 || b-u < tol2)
					d=SIGN(tol1,xm-x);
			}
		} else {
			d=CGOLD*(e=(x >= xm ? a-x : b-x));
		}
		u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
		fu=(*f)(u);
		/* This is te one function evaluation per iteration. */
		if (fu <= fx) {	/* Now decide what to do with our function evaluation. */
			if (u >= x) a=x; else b=x;
			SHFT(v,w,x,u)		/* Housekeeping follows. */
			SHFT(fv,fw,fx,fu)
		} else {
			if (u < x) a=u; else b=u;
			if (fu <= fw || w == x) {
				v=w;
				w=u;
				fv=fw;
				fw=fu;
			} else if (fu <= fv || v == x || v == w) {
				v=u;
				fv=fu;
			}
		}		/* Done with housekeeping. Back for another iteration. */
	}
	cerr << "Too many iterations in brent";
	*xmin=x;
	return fx;
}
#undef ITMAX
#undef CGOLD
#undef ZEPS
#undef SHFT
/*****************************************************************************************/