/*******************************************************
*        Find minimum of a real function Y=F(X)        *
* ---------------------------------------------------- *
* REFERENCE: "Numerical Recipes, The Art of Scientific *
*             Computing By W.H. Press, B.P. Flannery,  *
*             S.A. Teukolsky and W.T. Vetterling,      *
*             Cambridge University Press, 1986"        *
*             [BIBLI 08].                              *
* ---------------------------------------------------- *
* Sample run:                                          *
*                                                      *
* Find a minimum of function X*SIN(X)-2*COS(X)         *
* between X=4 and X=6.                                 *
*                                                      *
* X1=4  X2=5  X3=6  TOL=1e-6                           *
*                                                      *
* Function minimum is: -5.534529                       *
* for X= 5.232937                                      *
*                                                      *
*Press any key to continue                             *
*                                                      *
* ---------------------------------------------------- *
*                          C++ version by J-P Moreau.  * 
*                              (www.jpmoreau.fr)       *
*******************************************************/
#include <stdio.h>
#include <math.h>

#define ITMAX  100
#define CGOLD  0.3819660
#define ZEPS   1.0E-10
//Maximum allowed number of iterations; golden ration; and a small
//number which protects against trying to achieve fractional accuracy
//for a minimum that happens to be exactly zero.

double  x1,x2,x3,xmini,ymini,tol;

//Function to be analyzed
double F(double x) {
  return x*sin(x)-2.0*cos(x);
}

double MIN(double a, double b) {
  if (a<=b) return a; else return b; 
}

double MAX(double a, double b) { 
  if (a>=b) return a; else return b;
}

double Sign(double a, double b) {
  if (b>=0) return fabs(a);
  else return -fabs(a);
}

double BRENT(double ax, double bx, double cx, double tol, double *xmin) {
// Labels: e1,e2,e3
  double a,b,d,e,fu,fv,fx,fw,u,v,w,x;
  double etemp,p,q,r,tol1,tol2,xm;
  int    iter;
/*Given a function F, and a bracketing triplet of abscissas
  AX,BX,CX (such that BX is between AX and CX, and F(BX) 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 in XMIN, and the minimum
  function value is returned as BRENT, the returned function value.*/
  a=MIN(ax,cx);
  b=MAX(ax,cx);
  v=bx; w=v; x=v; e=0;
  fx=F(x); fv=fx; fw=fx;
  d=0;
  for (iter=1; iter<=ITMAX; iter++) {                        // main loop
    xm=0.5*(a+b);
    tol1=tol*fabs(x)+ZEPS;
    tol2=2.0*tol1;
    if (fabs(x-xm)<=tol2-0.5*(b-a)) goto e3;         //Test for done here
    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);        //bug corrected 07/24/2006 (0.2 instead of 2.0)
      if (q>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)) goto e1;
//  The above conditions determine the acceptability of the 
//  parabolic fit. Here it is o.k.:      
      d=p/q;
      u=x+d;
      if (u-a<tol2 || b-u<tol2)  d=Sign(tol1,xm-x);
      goto e2;
    }
e1: if (x>=xm)  e=a-x;
    else  e=b-x;
    d=CGOLD*e;
e2: if (fabs(d)>=tol1)  u=x+d;
    else  u=x+Sign(tol1,d);
    fu=F(u);    // This the one function evaluation per iteration
    if (fu<=fx) {
      if (u>=x)  a=x;
      else  b=x;
      v=w; fv=fw;
      w=x; fw=fx;
      x=u; fx=fu;
    }
    else {
      if (u<x)  a=u;
      else  b=u;
      if (fu<=fw || w==x) {
        v=w; fv=fw;
        w=u; fw=fu;
      }
      else if (fu<=fv || v==x || v==w) {
        v=u; fv=fu;
	  }
    }
  }
  printf(" Brent exceed maximum iterations.\n");
e3:*xmin=x;   //exit section
  return fx;
} //BRENT()


void main()  {

  x1=4.0; x2=5.0; x3=6.0;
  tol=1e-6;

  ymini = BRENT(x1,x2,x3,tol,&xmini);

  printf("\n Function minimum is %f\n\n",ymini);
  printf(" for X= %f\n\n",xmini);

}

  
// end of file Brent.cpp