/*******************************************************
*        Find minimum of a real function Y=F(X)        *
*        using routine for Golden Section Search       *
* ---------------------------------------------------- *
* 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                                      *
*                                                      *
* ---------------------------------------------------- *
*                           C++ version by J-P Moreau. *
*                               (www.jpmoreau.fr)      *
*******************************************************/
#include <stdio.h>
#include <math.h>

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

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


double GOLDEN(double ax,double bx,double cx,double tol,double *xmin) {
/*------------------------------------------------------------------
  Given a function F, and given a bracketing triplet of fabscissas 
  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 performs a golden section 
  search for the minimum, isolating it to a fractional precision of 
  about tol. The fabscissa of the minimum is returned as XMIN, and the minimum
  function value is returned as GOLDEN, the returned function value.
--------------------------------------------------------------------*/
  double f0,f1,f2,f3,x0,x1,x2,x3;
  double R,C;
  R=0.61803399; C=1.0-R;    // golden ratios
  x0=ax;  //At any given time we will keep trace of 4 points: x0,x1,x2,x3
  x3=cx;
  if (fabs(cx-bx) > fabs(bx-ax)) {
    x1=bx; x2=bx+C*(cx-bx);
  }
  else {
    x2=bx; x1=bx-C*(bx-ax);
  };
  f1=F(x1); f2=F(x2);    //Initial function evaluations
  //main loop
  while (fabs(x3-x0) > tol*(fabs(x1)+fabs(x2))) {
    if (f2 < f1) {
      x0=x1; x1=x2;
      x2=R*x1+C*x3;
      f0=f1; f1=f2;
      f2=F(x2);
    }
    else {
      x3=x2; x2=x1;
      x1=R*x2+C*x0;
      f3=f2; f2=f1;
      f1=F(x1);
    }
  }

  if (f1 < f2) {
	*xmin=x1;
    return f1;
  }
  else {
    *xmin=x2;
	return f2;
  }
} // Function Golden()


void main() {

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

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

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

}

  
//end of file Golden.cpp