'*****************************************************
'*      Program to demonstrate the real domain       *
'*               Zbrent subroutine                   *
'* ------------------------------------------------- *
'* Reference:  BORLAND MATHEMATICAL LIBRARY          *
'*                                                   *
'*              Basic version by J-P Moreau, Paris.  *
'*                       (www.jpmoreau.fr)           *
'* ------------------------------------------------- *
'* Example:    Find a real root of f(x)=(x+1)^5      *
'*                                                   *
'* SAMPLE RUN:                                       *
'*                                                   *
'*  Input interval (X1,X2):                          *
'*                                                   *
'*        X1 = -2                                    *
'*        X2 =  0                                    *
'*                                                   *
'*  Convergence criterion: 1e-10                     *
'*  Maximum number of iterations: 10                 *
'*                                                   *
'*  The estimated root is:                           *
'*                                                   *
'*        X = -1                                     *
'*                                                   *
'*  The associated Y value is Y = 0                  *
'*                                                   *
'*  The number of iterations was: 2                  *
'*  The error code is: 0                             *
'*                                                   *
'*****************************************************
DEFINT I-N
DEFDBL A-H, O-Z

  CLS
  PRINT
  PRINT " Input interval (X1,X2):"
  PRINT
  INPUT "       X1 = ", x1
  INPUT "       X2 = ", x2
  PRINT
  INPUT " Convergence criterion: ", e
  INPUT " Maximum number of iterations: ", maxiter

  GOSUB 2000  'Call BrentRoots subroutine

  PRINT
  PRINT " The estimated root is:"
  PRINT
  PRINT "       X = "; result
  PRINT
  PRINT " The associated Y value is Y = "; yr
  PRINT
  PRINT " The number of iterations was: "; niter
  PRINT " The error code is: "; ierror
  PRINT

END


1000 'Test function for Brent method
  AFunction = (x + 1)^5
RETURN

1200 'TRUE if FA*FB negative
  IF FA * FB < 0 THEN
    IRootBracketed = 1
  ELSE
    IRootBracketed = 0
  END IF
RETURN

1400 'returns the minimum of two real numbers xx1 and xx2
  IF xx1 <= xx2 THEN
    xminimum = xx1
  ELSE
    xminimum = xx2
  END IF
RETURN


'*****************************************************
'*              Brent Method Function                *
'* ------------------------------------------------- *
'* The purpose is to find a real root of a real      *
'* function f(x) using Brent method.                 *
'*                                                   *
'* INPUTS:  x1,x2     : interval of root             *
'*          Tolerance : desired accuracy for root    *
'*          maxIter   : maximum number of iterations *
'*                                                   *
'* OUTPUTS: The function returns the root value      *
'*          ValueAtRoot : value of f(root)           *
'*          niter    : number of done iterations     *
'*          ierror   : =0, all OK                    *
'*                   : =1, no root found in interval *
'*                   : =2, no more iterations !      *
'*****************************************************
2000 FPP = 1E-11
  xnearzero = 1E-20
  i = 0: ierror = 0
  'Get values of test function for x=x1 and x=x2
  x = x1: GOSUB 1000: FA = AFunction
  x = x2: GOSUB 1000: FB = AFunction
  'Test if a root exists in the interval (x1,x2)
  GOSUB 1200 'f(x1)*f(x2) not negative?
  IF IRootBracketed = 0 THEN
    ierror = 1: RETURN
  END IF
  FC = FB : AA = x1 : BB = x2
2100 'main loop
  If FB*FC > 0 THEN
    CC=AA : FC=FA : DD = BB - AA : EE = DD
  END IF
  IF ABS(FC) < ABS(FB) THEN
    AA = BB: BB = CC: CC = AA
    FA = FB: FB = FC: FC = FA
  END IF
  Tol1 = 2# * FPP * ABS(BB) + .5# * e
  xm = .5# * (CC - BB)
  IF (ABS(xm) <= Tol1) OR (ABS(FA) < xnearzero) THEN
    result = BB
    x = result: GOSUB 1000: yr = AFunction
    GOTO 2500
  END IF
  IF (ABS(EE) >= Tol1) AND (ABS(FA) > ABS(FB)) THEN
    SS = FB / FA
    IF ABS(AA - CC) < xnearzero THEN
      PP = 2# * xm * SS
      QQ = 1# - SS
    ELSE
      QQ = FA / FC
      RR = FB / FC
      PP = SS * (2# * xm * QQ * (QQ - RR) - (BB - AA) * (RR - 1#))
      QQ = (QQ - 1#) * (RR - 1#) * (SS - 1#)
    END IF
    IF PP > xnearzero THEN QQ = -QQ
    PP = ABS(PP)
    xx1 = 3# * xm * QQ - ABS(Tol1 * QQ): xx2 = ABS(EE * QQ): GOSUB 1400
    IF 2# * PP < xminimum THEN
      EE = DD:   DD = PP / QQ
    ELSE
      DD = xm: EE = DD
    END IF
  ELSE
    DD = xm
    EE = DD
  END IF
  AA = BB
  FA = FB
  IF ABS(DD) > Tol1 THEN
    BB = BB + DD
  ELSE
    IF xm > 0 THEN
      BB = BB + ABS(Tol1)
    ELSE
     BB = BB - ABS(Tol1)
    END IF
  END IF
  x = BB: GOSUB 1000: FB = AFunction
  i = i+1
  IF i >= maxiter THEN
    ierror = 2: GOTO 2500
  END IF
  'Continue
  GOTO 2100
2500 'A real root has been found
  niter = i
RETURN

'End of file zbrent.bas