/***************************************************************

   FUNCTION: brent

   This function is translated from the Pascal version in
   Numerical Recipes, Press et al., 1986.
   It returns the minimum value of a user-supplied function,
   in a specified interval known to contain a minimum.

*/

#include <stdio.h>
#include <math.h>    /* fabs() */
#include "dbhost.h"  /* PROTOTYPE_ALLOWED */

#if PROTOTYPE_ALLOWED
double sign(double a, double b);
double (*func)(double x);
#else
double sign();
double (*func)();
#endif

#if PROTOTYPE_ALLOWED
double sign(double a, double b)
#else
double sign(a,b)
double a,b;
#endif
{
   return (b > 0.0) ? fabs(a) : -fabs(a);
}

/* Programs using routine BRENT must supply an external function
   func(x:real):real whose minimum is to be found. */

#define ITMAX 50
#define CGOLD 0.3819660
#define ZEPS  1.0e-10

#if PROTOTYPE_ALLOWED
double brent(double ax, double bx, double cx, double (*func)(), double tol,
  double *xmin)
#else
double brent(ax,bx,cx,func,tol,xmin)
double ax, bx, cx, tol, *xmin;
double (*func) ();
#endif
{
   double  a,b,d=0.0,e,etemp;
   double  fu,fv,fw,fx;
   int     iter;
   double  p,q,r,tol1,tol2;
   double  u,v,w,x,xm;

   a = (ax < cx) ? ax : cx;
   b = (ax > cx) ? ax : cx;
   v = bx;
   w = v;
   x = v;
   e = 0.0;
   fx = (*func)(x);
   fv = fx;
   fw = fx;
   for(iter = 1; iter <= ITMAX; iter++)
   {
      xm = 0.5*(a+b);
      tol1 = tol*fabs(x)+ZEPS;
      tol2 = 2.0*tol1;
      if( (fabs(x-xm) <= (tol2-0.5*(b-a))) ) {
         goto step3;
         }
      if(fabs(e) > tol1)
      {
         r = (x-w)*(fx-fv);
         q = (x-v)*(fx-fw);
         p = (x-v)*q-(x-w)*r;
         q = 2.0*(q-r);
         if(q) 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 step1;
         d = p/q;
         u = x+d;
         if( ((u-a) < tol2) || ((b-u) < tol2) )
             d = sign(tol1,xm-x);
         goto step2;
      }

   step1:
      e = (x >= xm) ? (a-x) : (b-x);
      d = CGOLD*e;

   step2:
      if(fabs(d) >= tol1)  u = x+d;
      else u = x+sign(tol1,d);
      fu = (*func)(u);
      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 == 2) )
          {
             v = u;
             fv = fu;
          }
        }
      }
   }
   printf("Warning: max (%d) iterations in brent (function minimizer)\n",ITMAX);

  step3:
     *xmin = x;
     return(fx);
}

#if 0
/* Test routine: */

#if PROTOTYPE_ALLOWED
double polyfunc(double x)
#else
double polyfunc(x)
double x;
#endif
{
    return((x)*(x) + 2*(x)-3);
}

#if PROTOTYPE_ALLOWED
void main(void)
#else
void main()
#endif
{
   double y, xmin;
   y =  brent(-2.5, 0.0, 1.0, polyfunc, 0.001, &xmin);

   printf(" y = %f  xmin = %f\n",y,xmin);
}
#endif