      SUBROUTINE ZBRENTX(X,A,B,T,ERRMSG)

C MICHAEL B. PORTER  8/84

C FORTRAN CONVERSION OF ALGOL PROGRAM PUBLISHED IN
C    THE COMPUTER JOURNAL 14(4):422-425 (1971)
C    BY R. P. BRENT

C RETURNS A ZERO X OF THE FUNCTION F IN THE GIVEN INTERVAL
C    [A,B], TO WITHIN A TOLERANCE 6*MACHEP*ABS(X)+2*T, WHERE
C    MACHEP IS THE RELATIVE MACHINE PRECISION AND T IS A POSITIVE
C    TOLERANCE.  THE PROCEDURE ASSUMES THAT FUNCT(A) AND FUNCT(B) HAVE
C    DIFFERENT SIGNS.

C    THIS IS THE EXTENDED RANGE VERSION WHICH EXPECT THE FUNCTION
C    TO HAVE THE FORM
C         SUBROUTINE FUNCT(X,G,IPOW)
C    WHERE G*10**IPOW GIVES FUNCT(X)

      IMPLICIT REAL*8 (A-H, O-Z)
      REAL*8 MACHEP,M
      CHARACTER*80 ERRMSG

      ERRMSG = ' '
      MACHEP = 1.0E-16
      TEN = 10.0

      CALL FUNCT(A,FA,IEXPA)
      CALL FUNCT(B,FB,IEXPB)

      IF ( ( (FA .GT. 0.0) .AND. (FB .GT. 0.0) ) .OR.
     &     ( (FA .LT. 0.0) .AND. (FB .LT. 0.0) ) ) THEN
         ERRMSG = ' *** ZBRENT ERROR: FUNCT SGN SAME AT INTRVL ENDPTS'
         RETURN
      ENDIF

C     INTERNAL ROOT

 2000 C = A
      FC = FA
      IEXPC = IEXPA
      E = B-A
      D = E

C     EXTERNAL ROOT

      IF (IEXPA .LT. IEXPB) THEN
         F1 = FC*TEN**(IEXPC-IEXPB)
         F2 = FB
      ELSE
         F1 = FC
         F2 = FB*TEN**(IEXPB-IEXPC)
      ENDIF
 3000 IF (ABS(F1) .LT. ABS(F2)) THEN
         A = B
         B = C
         C = A
         FA = FB
         IEXPA = IEXPB
         FB = FC
         IEXPB = IEXPC
         FC = FA
         IEXPC = IEXPA
      ENDIF
      TOL = 2.0*MACHEP*ABS(B)+T
      M = 0.5*(C-B)
      IF ((ABS(M) .GT. TOL) .AND. (FB .NE. 0.0)) THEN

C     ------ SEE IF A BISECTION IS FORCED
      IF (IEXPA .LT. IEXPB) THEN
         F1 = FA*TEN**(IEXPA-IEXPB)
         F2 = FB
      ELSE
         F1 = FA
         F2 = FB*TEN**(IEXPB-IEXPA)
      ENDIF

         IF ( (ABS(E) .LT. TOL) .OR. 
     &       (ABS( F1 ) .LE. ABS(F2)) ) THEN
            E = M
            D = E
         ELSE
            S = FB/FA*TEN**(IEXPB-IEXPA)
            IF (A .EQ. C) THEN
C              ------ LINEAR INTERPOLATION
               P = 2.0*M*S
               Q = 1.0-S
            ELSE
C              ------ INVERSE QUADRATIC INTERPOLATION
               Q = FA/FC*TEN**(IEXPA-IEXPC)
               R = FB/FC*TEN**(IEXPB-IEXPC)
               P = S*(2.0*M*Q*(Q-R)-(B-A)*(R-1.0))
               Q = (Q-1.0)*(R-1.0)*(S-1.0)
            ENDIF
            IF (P .GT. 0.0) THEN
               Q = -Q
            ELSE
               P = -P
            ENDIF
            S = E
            E = D
            IF ((2.0*P .LT. 3.0*M*Q-ABS(TOL*Q)) .AND.
     &          (P .LT. ABS(0.5*S*Q))) THEN
               D = P/Q
            ELSE
               E = M
               D = E
            ENDIF
         ENDIF

         A = B
         FA = FB
         IEXPA = IEXPB

         IF (ABS(D) .GT. TOL) THEN
            B = B+D
         ELSE
            IF (M .GT. 0.0) THEN
               B = B+TOL
            ELSE
               B = B-TOL
            ENDIF
         ENDIF

         CALL FUNCT(B,FB,IEXPB)
         IF ((FB .GT. 0.0) .EQV. (FC .GT. 0.0)) GOTO 2000
         GOTO 3000
      ENDIF

      X = B

      RETURN
      END
