/*******************************************************
*    Bracketing a minimum of a real function Y=F(X)    *
*             using MNBRAK subroutine                  *
* ---------------------------------------------------- *
* 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 3 points of function X*SIN(X)-2*COS(X)          *
* bracketing a minimum, given two initial points.      *
*                                                      *
* X1=4.0  X2=5.0                                       *
*                                                      *
* The three points are:                                *
*                                                      *
* X1=  4.000000  X2=  5.000000  X3=  6.618034          *
*                                                      *
* Corresponding function values:                       *
*                                                      *
* F1= -1.719923  F2= -5.361946  F3=  0.285940          *
*                                                      *
*                          C++ version by J-P Moreau.  *
*                              (www.jpmoreau.fr)       *
*******************************************************/
#include <stdio.h>
#include <math.h>

/*The first parameter is the default ratio by which successive intervals
  are magnified; the second is the maximum magnification allowed for
  a parabolic-fit step.     */
#define  GOLD    1.618034
#define  GLIMIT  100
#define  TINY    1e-20

double  x1, x2, x3, f1, f2, f3;


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

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);
}


void MNBRAK(double ax, double bx, double *cx, double *fa, 
			double *fb, double *fc)  {
/*Given a function FUNC(X), and given distinct initial 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 which bracket a minimum of the function. Also returned
  are the function values at the three points, fa, fb and fc. */
// LABEL: e1
  double dum,fu,q,r,u,ulim;
  *fa=FUNC(ax);
  *fb=FUNC(bx);
  if (*fb > *fa) {
    dum=ax;
    ax=bx;
    bx=dum;
    dum=*fb;
    *fb=*fa;
    *fa=dum;
  }
  *cx=bx+GOLD*(bx-ax);
  *fc=FUNC(*cx);
e1:if (*fb>=*fc) {
    r=(bx-ax)*(*fb-*fc);
    q=(bx-*cx)*(*fb-*fa);
    u=bx-((bx-*cx)*q-(bx-ax)*r)/(2.*Sign(Max(fabs(q-r),TINY),q-r));
    ulim=bx+GLIMIT*(*cx-bx);
    if ((bx-u)*(u-*cx)>0) {
      fu=FUNC(u);
      if (fu < *fc) {
        ax=bx;
        *fa=*fb;
        bx=u;
        *fb=fu;
        goto e1;
      }
      else if (fu > *fb) {
        *cx=u;
        *fc=fu;
        goto e1;
      }
      u=*cx+GOLD*(*cx-bx);
      fu=FUNC(u);
    }
    else if ((*cx-u)*(u-ulim) > 0) {
      fu=FUNC(u);
      if (fu < *fc) {
        bx=*cx;
        *cx=u;
        u=*cx+GOLD*(*cx-bx);
        *fb=*fc;
        *fc=fu;
        fu=FUNC(u);
      }
    }
    else if ((u-ulim)*(ulim-*cx)>=0) {
      u=ulim;
      fu=FUNC(u);
	}
    else {
      u=*cx+GOLD*(*cx-bx);
      fu=FUNC(u);
    }
    ax=bx;
    bx=*cx;
    *cx=u;
    *fa=*fb;
    *fb=*fc;
    *fc=fu;
    goto e1;
  }
}


void main() {

  x1=4.0; x2=5.0;

  MNBRAK(x1,x2,&x3,&f1,&f2,&f3);

  printf("\n The three points are:\n\n");
  printf(" X1=%10.6f  X2=%10.6f  X3=%10.6f\n",x1,x2,x3);
  printf("\n Corresponding function values:\n\n");
  printf(" F1=%10.6f  F2=%10.6f  F3=%10.6f\n\n\n",f1,f2,f3);

}

// end of file mnbrak.cpp