'********************************************************
'*        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" [8].   *
'* ---------------------------------------------------- *
'* 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                             *
'*                                                      *
'* ---------------------------------------------------- *
'*                         BASIC version by J-P Moreau. *
'*                               (www.jpmoreau.fr)      *
'********************************************************
defdbl a-h,o-z

ITMAX=100
CGOLD=0.3819660
ZEPS=1E-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.

ax=4.0 : bx=5.0 : cx=6.0
tol=1e-6

gosub 2000  'call Brent method

cls
print
print using " Function minimum is ##.######"; fx
print
print using " for X= ##.######"; x
print : print

end

'Function to be analyzed
1000 F=xx*sin(xx)-2*cos(xx)
return

'double MIN(double a, double b)
1100 if a<=b then
    xx=a
  else
    xx=b
  end if
return

'double MAX(double a, double b)
1200  if a>=b then
    xx=a
  else
    xx=b
  end if
return

'double Sign(double a, double b)
1300 if b>=0 then
    xx=abs(a)
  else
    xx=-abs(a)
  end if
return

'double BRENT(double ax, double bx, double cx, double tol, double *xmin)
'------------------------------------------------------------------
' Labels: 1,2
' 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.
'------------------------------------------------------------------
2000:
  a=ax : b = cx : gosub 1100 : a=xx  'a=MIN(ax,cx)
  gosub 1200 : b=xx                  'b=MAX(ax,cx)
  v=bx: w=v: x=v: e=0
  xx=x : gosub 1000 : fx=F           'fx=F(x)
  fv=fx : fw=fx
  d=0
  for iter=1 to ITMAX                        'main loop
    xm=0.5*(a+b)
    tol1=tol*abs(x)+ZEPS
    tol2=2.0*tol1
    if abs(x-xm)<=tol2-0.5*(b-a) then return         'Test for done here
    if abs(e)>tol1 then                 'Construct a trial parabolic fit
      r=(x-w)*(fx-fv)
      q=(x-v)*(fx-fw)
      p=(x-v)*q-(x-w)*r
      q=2#*(q-r)          'bug corrected 07/24/2006 (0.2 instead of 2.0)
      if q>0 then  p=-p
      q=abs(q)
      etemp=e
      e=d
      if abs(p)>=abs(0.5*q*etemp) or p<=q*(a-x) or p>=q*(b-x) then goto 1
'   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 or b-u<tol2 then
        a=tol1:b=xm-x:gosub 1300:d=xx    'd=Sign(tol1,xm-x)
      end if
      goto 2
    end if
1   if x>=xm then
      e=a-x
    else
      e=b-x
    end if
    d=CGOLD*e
2   if abs(d)>=tol1 then
      u=x+d
    else
      a=tol1:b=d:gosub 1300: u=x+xx      'u=x+Sign(tol1,d)
    end if
    xx=u : gosub 1000 : fu=F             'fu=F(u)
    ' This the one function evaluation per iteration
    if fu<=fx then
      if u>=x then
        a=x
      else
        b=x
      end if
      v=w: fv=fw
      w=x: fw=fx
      x=u: fx=fu
    else
      if u<x then
        a=u
      else
        b=u
      end if
      if fu<=fw or w=x then
        v=w: fv=fw
        w=u: fw=fu
      else
        if fu<=fv or v=x or v=w then
          v=u: fv=fu
        end if
      end if
    end if
  next iter
  print " Brent exceed 100 iterations."
  print : print


' end of file Brent.bas