C      ALGORITHM 748, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 21, NO. 3, September, 1995, P.  327-344.
C
C This file contains 5 files separated by lines of the form
C         C*** filename
C
C The filenames in this file are:
C
C driver.f             enclofx.f            exdrive.f           
C testdata             testout                                  
C
C*** driver.f
C        TEST PROBLEM SET ATTACHED TO THE ALGORITHM "ENCLOFX"
C        WITH THE COMPANION PAPER "ENCLOSING ZEROS OF CONTINUOUS
C        FUNCTIONS" BY
C        G.E.ALEFELD, F.A.POTRA, AND YIXUN SHI.
C
C        THE FOLLOWING RUNS THE SET OF TEST PROBLEMS LISTED IN THE
C        COMPANION PAPER. THE RESULT OF THE TEST RUN IS ALSO LISTED
C        IN THAT PAPER.
C
C ******************************************************************
         INTEGER NPROB,NEPS,ISIGN
	 INTEGER N
         DOUBLE PRECISION A,B,EPS,ROOT
         EXTERNAL ISIGN
	 COMMON /USER/ N
C THIS PROGRAM CALCULATES A ROOT OF A CONTINUOUS FUNCTION F(X)         
C IN AN INTERVAL I=[A, B] PROVIDED THAT F(A)F(B)<0. THE FUNCTION F(X)
C AND THE INITIAL INTERVAL [A, B] ARE TO BE SUPPLIED BY THE USER IN 
C THE SUBROUTINES "FUNC" AND "INIT". THE OUTPUT "ROOT" EITHER SATISFIES
C THAT F(ROOT)=0 OR IS AN APPROXIMATE SOLUTION OF THE EQUATION F(X)=0
C SUCH THAT "ROOT" IS INCLUDED IN AN INTERVAL [AC, BC] WITH
C      F(AC)F(BC)<0,
C AND
C      BC-AC <= TOL = 2*TOLE(AC,BC).
C PRECISION CHOSEN AS THE RELATIVE MACHINE PRECISION,
C AND "UC" IS EQUAL TO EITHER "AC" OR "BC" WITH THE CONDITION
C      |F(UC)| = MIN{ |F(AC)|, |F(BC)| }.
C
C INPUT OF THE PROGRAM:
C  NPROB -- INTEGER. POSITIVE INTEGER STANDING FOR "NUMBER OF PROBLEM",
C           INDICATING THE PROBLEM TO BE SOLVED.
C  N     -- PROBLEM DEPENDENT PARAMETER
C OUTPUT OF THE PROGRAM:
C  ROOT  -- DOUBLE PRECISION. EXACT OR APPROXIMATE SOLUTION OF THE
C           EQUATION F(X)=0.

10       READ(*,*,END=20) NPROB, N
C
C USE MACHINE PRECISION
C
	 CALL RMP(EPS)
	 NEPS = 1000
C
C CALL SUBROUTINE "INIT" TO GET THE INITIAL INTERVAL.
C 
         CALL INIT(NPROB, A, B)
C
C CALL SUBROUTINE "RROOT" TO HAVE THE PROBLEM SOLVED.
C
         CALL RROOT(NPROB,NEPS,EPS,A,B,ROOT)
         CONTINUE
C
C PRINT OUT THE RESULT ON THE SCREEN.
C
         PRINT*, 'PROBLEM NUMBER ', NPROB, 
     +           ' WITH PROBLEM PARAMETER ', N
         PRINT*, 'COMPUTED ROOT =', ROOT
	 PRINT*, ' '
	 GOTO 10
20       STOP
         END

         SUBROUTINE INIT(NPROB, A, B)
         INTEGER NPROB
         DOUBLE PRECISION A,B,PI

         PI=3.1416D0
         GOTO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,
     *    22,23,24,25,26,27,28), NPROB

1        A=PI/2.0D0
         B=PI
         RETURN

2        A=1.0D0+1.0D-9
         B=(2.0D0)*(2.0D0)-1.0D-9
         RETURN

3        A=(2.0D0)*(2.0D0)+1.0D-9
         B=(3.0D0)*(3.0D0)-1.0D-9
         RETURN

4        A=(3.0D0)*(3.0D0)+1.0D-9
         B=(4.0D0)*(4.0D0)-1.0D-9
         RETURN

5        A=(4.0D0)*(4.0D0)+1.0D-9
         B=(5.0D0)*(5.0D0)-1.0D-9
         RETURN

6        A=(5.0D0)*(5.0D0)+1.0D-9
         B=(6.0D0)*(6.0D0)-1.0D-9
         RETURN

7        A=(6.0D0)*(6.0D0)+1.0D-9
         B=(7.0D0)*(7.0D0)-1.0D-9
         RETURN

8        A=(7.0D0)*(7.0D0)+1.0D-9
         B=(8.0D0)*(8.0D0)-1.0D-9
         RETURN

9        A=(8.0D0)*(8.0D0)+1.0D-9
         B=(9.0D0)*(9.0D0)-1.0D-9
         RETURN

10       A=(9.0D0)*(9.0D0)+1.0D-9
         B=(10.0D0)*(10.0D0)-1.0D-9
         RETURN

11       A=(10.0D0)*(10.0D0)+1.0D-9
         B=(11.0D0)*(11.0D0)-1.0D-9
         RETURN

12       CONTINUE

13       CONTINUE

14       A=-9.0D0
         B=31.0D0
         RETURN

15       CONTINUE

16       A=0.0D0
         B=5.0D0
         RETURN

17      A=-0.95D0
        B=4.05D0
        RETURN

18       A=0.0D0
         B=1.5D0
         RETURN

19       CONTINUE

20       CONTINUE

21       CONTINUE

22       CONTINUE

23       A=0.0D0
         B=1.0D0
         RETURN

24       A=1.0D-2
         B=1.0D0
         RETURN

25       A=1.0D0
         B=100.0D0
         RETURN

26       A=-1.0D0
         B=4.0D0
         RETURN

27       A=-10000
         B=PI/2.0D0
         RETURN

28       A=-10000
         B=1.0D-4
         RETURN

         END

C        *****************

         SUBROUTINE FUNC(NPROB, X, FX)
         INTEGER NPROB, N,I
         DOUBLE PRECISION X, FX, DN, DFLOAT, DI
	 COMMON /USER/ N

C         N --- AN INTEGER USED IN THE FUNCTION. 
C

         DN=DFLOAT(N)
         
         GOTO (10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,
     *         160,170,180,190,200,210,220,230,240,250,260,
     *         270,280), NPROB

10       FX=DSIN(X)-X/2.0D0
         RETURN

20       CONTINUE
30       CONTINUE
40       CONTINUE
50       CONTINUE
60       CONTINUE
70       CONTINUE
80       CONTINUE
90       CONTINUE
100      CONTINUE
110      CONTINUE
         FX=0.0D0
         DO 115 I=1,20
           DI=DFLOAT(I)
           FX=FX+((2.0D0*DI-5.0D0)**2)/(X-DI*DI)**3
115      CONTINUE
         FX=-2.0D0*FX
         RETURN
120      FX=-40.0D0*X*DEXP(-1.0D0*X)
         RETURN

130      FX=-100.0D0*X*DEXP(-2.0D0*X)
         RETURN

140      FX=-200.0D0*X*DEXP(-3.0D0*X)
         RETURN

150      FX=X**N-0.2D0
         RETURN

160      CONTINUE
170      FX=X**N-1.0D0
         RETURN

180      FX=DSIN(X)-0.5D0
         RETURN

190      FX=2.0D0*X*DEXP(-DN)-2.0D0*DEXP(-DN*X)+1.0D0
         RETURN

200      FX=(1.0D0+(1.0D0-DN)**2)*X - (1.0D0-DN*X)**2
         RETURN

210      FX=X**2-(1.0D0-X)**N
         RETURN

220      FX=(1.0D0+(1.0D0-DN)**4)*X - (1.0D0-DN*X)**4
         RETURN

230      FX=(X-1.0D0)*DEXP(-DN*X)+X**N
         RETURN

240      FX=(DN*X-1.0D0)/((DN-1.0D0)*X)
         RETURN

250      FX=X**(1.0D0/DN)-DN**(1.0D0/DN)
         RETURN

260      CONTINUE
         IF(X .EQ. 0.0D0)THEN
           FX=0.0D0
         ELSE
           FX=X/DEXP(1.0D0/(X*X))
         ENDIF
         RETURN

270      CONTINUE
         IF(X .GE. 0.0D0)THEN
           FX=(X/1.5D0 + DSIN(X) -1.0D0)*DN/20.0D0
         ELSE
           FX=(-1.0D0*DN)/20.0D0
         ENDIF
         RETURN

280      CONTINUE
         IF(X .GE. (1.0D-3)*2.0D0/(DN+1.0D0)) THEN
           FX=DEXP(1.0D0)-1.859D0
         ELSEIF (X .GE. 0.0D0)THEN
            FX=DEXP((DN+1.0D0)*0.5D0*X*1.0D3)-1.859D0
         ELSE
           FX=-0.859D0
         ENDIF
         RETURN

         END
C*** enclofx.f
         SUBROUTINE RROOT(NPROB,NEPS,EPS,A,B,ROOT)
C FINDS EITHER AN EXACT SOLUTION OR AN APPROXIMATE SOLUTION OF THE
C EQUATION F(X)=0 IN THE INTERVAL [A,B]. AT THE BEGINING OF EACH
C ITERATION, THE CURRENT ENCLOSING INTERVAL IS RECORDED AS [A0,B0].
C THE FIRST ITERATION IS SIMPLY A SECANT STEP. STARTING WITH THE
C SECOND ITERATION, THREE STEPS ARE TAKEN IN EACH ITERATION. FIRST
C TWO STEPS ARE EITHER QUADRATIC INTERPOLATION OR CUBIC INVERSE
C INTERPOLATION. THE THIRD STEP IS A DOUBLE-SIZE SECANT STEP. IF THE
C DIAMETER OF THE ENCLOSING INTERVAL OBTAINED AFTER THOSE THREE STEPS
C IS LARGER THAN 0.5*(B0-A0), THEN AN ADDITIONAL BISECTION STEP WILL
C BE TAKEN.
C  NPROB -- INTEGER. INDICATING THE PROBLEM TO BE SOLVED;
C  NEPS  -- INTEGER. USED TO DETERMINE THE TERMINATION CRITERION;
C  EPS   -- DOUBLE PRECISION. USED IN THE TERMINATION CRITERION;
C  A,B   -- DOUBLE PRECISION. INPUT AS THE INITIAL INTERVAL AND 
C           OUTPUT AS THE ENCLOSING INTERVAL AT THE TERMINATION;
C  ROOT  -- DOUBLE PRECISION. OUTPUT SOLUTION OF THE EQUATION.
         INTEGER NPROB,ITNUM,ISIGN,NEPS
         DOUBLE PRECISION A,B,FA,FB,C,U,FU,MU,A0,B0,TOL,D,FD
         DOUBLE PRECISION PROF,E,FE,EPS,ROOT
         EXTERNAL ISIGN
         PARAMETER (MU=0.5D0)
C
C INITIALIZATION. SET THE NUMBER OF ITERATION AS 0. CALL SUBROUTINE
C "FUNC" TO OBTAIN THE INITIAL FUNCTION VALUES F(A) AND F(B). SET
C DUMB VALUES FOR THE VARIABLES "E" AND "FE".
C
         ITNUM=0
         CALL FUNC(NPROB,A,FA)
         CALL FUNC(NPROB,B,FB)
         E=1.0D5
         FE=1.0D5
C
C ITERATION STARTS. THE ENCLOSING INTERVAL BEFORE EXECUTING THE
C ITERATION IS RECORDED AS [A0, B0].
C
 10      A0=A
         B0=B
C
C UPDATES THE NUMBER OF ITERATION.
C
         ITNUM=ITNUM+1
C
C CALCULATES THE TERMINATION CRITERION. STOPS THE PROCEDURE IF THE
C CRITERION IS SATISFIED.
C
         IF(DABS(FB) .LE. DABS(FA)) THEN
          CALL TOLE(B,TOL,NEPS,EPS)
         ELSE
          CALL TOLE(A,TOL,NEPS,EPS)
         ENDIF
         IF((B-A).LE.TOL)GOTO 400
C
C FOR THE FIRST ITERATION, SECANT STEP IS TAKEN.
C
         IF(ITNUM .EQ. 1)THEN
          C=A-(FA/(FB-FA))*(B-A)
C
C CALL SUBROUTINE "BRACKT" TO GET A SHRINKED ENCLOSING INTERVAL AS
C WELL AS TO UPDATE THE TERMINATION CRITERION. STOP THE PROCEDURE
C IF THE CRITERION IS SATISFIED OR THE EXACT SOLUTION IS OBTAINED.
C
          CALL BRACKT(NPROB,A,B,C,FA,FB,TOL,NEPS,EPS,D,FD)
          IF((FA.EQ.0.0D0).OR.((B-A).LE.TOL))GOTO 400
          GOTO 10
         ENDIF
C
C STARTING WITH THE SECOND ITERATION, IN THE FIRST TWO STEPS, EITHER
C QUADRATIC INTERPOLATION IS USED BY CALLING THE SUBROUTINE "NEWQUA"
C OR THE CUBIC INVERSE INTERPOLATION IS USED BY CALLING THE SUBROUTINE
C "PZERO". IN THE FOLLOWING, IF "PROF" IS NOT EQUAL TO 0, THEN THE
C FOUR FUNCTION VALUES "FA", "FB", "FD", AND "FE" ARE DISTINCT, AND
C HENCE "PZERO" WILL BE CALLED.
C
         PROF=(FA-FB)*(FA-FD)*(FA-FE)*(FB-FD)*(FB-FE)*(FD-FE)          
         IF((ITNUM .EQ. 2) .OR. (PROF .EQ. 0.0D0)) THEN
          CALL NEWQUA(A,B,D,FA,FB,FD,C,2)
         ELSE
          CALL PZERO(A,B,D,E,FA,FB,FD,FE,C)
          IF((C-A)*(C-B) .GE. 0.0D0)THEN
           CALL NEWQUA(A,B,D,FA,FB,FD,C,2)
          ENDIF
         ENDIF
         E=D
         FE=FD
C
C CALL SUBROUTINE "BRACKT" TO GET A SHRINKED ENCLOSING INTERVAL AS
C WELL AS TO UPDATE THE TERMINATION CRITERION. STOP THE PROCEDURE
C IF THE CRITERION IS SATISFIED OR THE EXACT SOLUTION IS OBTAINED.
C
         CALL BRACKT(NPROB,A,B,C,FA,FB,TOL,NEPS,EPS,D,FD)
         IF((FA.EQ.0.0D0).OR.((B-A).LE.TOL))GOTO 400
         PROF=(FA-FB)*(FA-FD)*(FA-FE)*(FB-FD)*(FB-FE)*(FD-FE)          
         IF(PROF .EQ. 0.0D0) THEN
          CALL NEWQUA(A,B,D,FA,FB,FD,C,3)
         ELSE
          CALL PZERO(A,B,D,E,FA,FB,FD,FE,C)
          IF((C-A)*(C-B) .GE. 0.0D0)THEN
           CALL NEWQUA(A,B,D,FA,FB,FD,C,3)
          ENDIF
         ENDIF
C
C CALL SUBROUTINE "BRACKT" TO GET A SHRINKED ENCLOSING INTERVAL AS
C WELL AS TO UPDATE THE TERMINATION CRITERION. STOP THE PROCEDURE
C IF THE CRITERION IS SATISFIED OR THE EXACT SOLUTION IS OBTAINED.
C
         CALL BRACKT(NPROB,A,B,C,FA,FB,TOL,NEPS,EPS,D,FD)
         IF((FA.EQ.0.0D0).OR.((B-A).LE.TOL))GOTO 400
         E=D
         FE=FD
C
C TAKES THE DOUBLE-SIZE SECANT STEP.
C
         IF (DABS(FA) .LT. DABS(FB))THEN
          U=A
          FU=FA
         ELSE
          U=B
          FU=FB
         ENDIF
         C=U-2.0D0*(FU/(FB-FA))*(B-A)
         IF(DABS(C-U) .GT. (0.5D0*(B-A)))THEN
          C=A+0.5D0*(B-A)
         ENDIF
C
C CALL SUBROUTINE "BRACKT" TO GET A SHRINKED ENCLOSING INTERVAL AS
C WELL AS TO UPDATE THE TERMINATION CRITERION. STOP THE PROCEDURE
C IF THE CRITERION IS SATISFIED OR THE EXACT SOLUTION IS OBTAINED.
C
         CALL BRACKT(NPROB,A,B,C,FA,FB,TOL,NEPS,EPS,D,FD)
         IF((FA.EQ.0.0D0).OR.((B-A).LE.TOL))GOTO 400
C
C DETERMINES WHETHER AN ADDITIONAL BISECTION STEP IS NEEDED. AND TAKES
C IT IF NECESSARY.
C
         IF((B-A) .LT. (MU*(B0-A0)))THEN
           GOTO 10
         ENDIF
         E=D
         FE=FD
C
C CALL SUBROUTINE "BRACKT" TO GET A SHRINKED ENCLOSING INTERVAL AS
C WELL AS TO UPDATE THE TERMINATION CRITERION. STOP THE PROCEDURE
C IF THE CRITERION IS SATISFIED OR THE EXACT SOLUTION IS OBTAINED.
C
         CALL BRACKT(NPROB,A,B,A+0.5D0*(B-A),FA,FB,TOL,NEPS,EPS,D,FD)         
         IF((FA .EQ. 0.0D0).OR.((B-A).LE.TOL))GOTO 400
         GOTO 10
C
C TERMINATES THE PROCEDURE AND RETURN THE "ROOT".
C
 400     CONTINUE
         ROOT=A
         RETURN
         END

         SUBROUTINE BRACKT(NPROB,A,B,C,FA,FB,TOL,NEPS,EPS,D,FD)
C GIVEN CURRENT ENCLOSING INTERVAL [A,B] AND A NUMBER C IN (A,B), IF
C F(C)=0 THEN SETS THE OUTPUT A=C. OTHERWISE DETERMINES THE NEW
C ENCLOSING INTERVAL: [A,B]=[A,C] OR [A,B]=[C,B]. ALSO UPDATES THE
C TERMINATION CRITERION CORRESPONDING TO THE NEW ENCLOSING INTERVAL.
C  NPROB   -- INTEGER. INDICATING THE PROBLEM TO BE SOLVED;
C  A,B     -- DOUBLE PRECISION. [A,B] IS INPUT AS THE CURRENT 
C             ENCLOSING INTERVAL AND OUTPUT AS THE SHRINKED NEW 
C             ENCLOSING INTERVAL;
C  C       -- DOUBLE PRECISION. USED TO DETERMINE THE NEW ENCLOSING
C             INTERVAL;
C  D       -- DOUBLE PRECISION. OUTPUT: IF THE NEW ENCLOSING INTERVAL
C             IS [A,C] THEN D=B, OTHERWISE D=A;
C  FA,FB,FD-- DOUBLE PRECISION. FA=F(A), FB=F(B), AND FD=F(D);
C  TOL     -- DOUBLE PRECISION. INPUT AS THE CURRENT TERMINATION 
C             CRITERION AND OUTPUT AS THE UPDATED TERMINATION 
C             CRITERION ACCORDING TO THE NEW ENCLOSING INTERVAL;
C  NEPS    -- INTEGER. USED TO DETERMINE THE TERMINATION CRITERION;
C  EPS     -- DOUBLE PRECISION. USED IN THE TERMINATION CRITERION.
         INTEGER NPROB,ISIGN,NEPS
         DOUBLE PRECISION A,B,C,FA,FB,FC,D,FD,TOL,EPS
         EXTERNAL ISIGN
C
C ADJUST C IF (B-A) IS VERY SMALL OR IF C IS VERY CLOSE TO A OR B.
C
         TOL=0.7D0*TOL
         IF((B-A) .LE. 2.0D0*TOL)THEN
          C=A+0.5D0*(B-A)
         ELSEIF(C .LE. A+TOL)THEN
          C=A+TOL
         ELSE
          IF(C .GE. B-TOL)THEN
           C=B-TOL
          ENDIF
         ENDIF
C
C CALL SUBROUTINE "FUNC" TO OBTAIN F(C)
C
         CALL FUNC(NPROB, C, FC)
C
C IF F(C)=0, THEN SET A=C AND RETURN. THIS WILL TERMINATE THE 
C PROCEDURE IN SUBROUTINE "RROOT" AND GIVE THE EXACT SOLUTION OF 
C THE EQUATION F(X)=0.
C
         IF(FC .EQ. 0.0D0)THEN
          A=C
          FA=0.0D0
          D=0.0D0
          FD=0.0D0
          RETURN
         ENDIF
C
C IF F(C) IS NOT ZERO, THEN DETERMINE THE NEW ENCLOSING INTERVAL.
C
         IF((ISIGN(FA)*ISIGN(FC)) .LT. 0)THEN
          D=B
          FD=FB
          B=C
          FB=FC
         ELSE
          D=A
          FD=FA
          A=C
          FA=FC
         ENDIF
C
C UPDATE THE TERMINATION CRITERION ACCORDING TO THE NEW ENCLOSING
C INTERVAL.
C
         IF(DABS(FB) .LE. DABS(FA)) THEN
          CALL TOLE(B,TOL,NEPS,EPS)
         ELSE
          CALL TOLE(A,TOL,NEPS,EPS)
         ENDIF
C
C END OF THE SUBROUTINE.
C
         RETURN
         END

         INTEGER FUNCTION ISIGN(X)
C INDICATES THE SIGN OF THE VARIABLE "X".
C  X     -- DOUBLE PRECISION.
C  ISIGN -- INTEGER.
         DOUBLE PRECISION X
         IF(X .GT. 0.0D0)THEN
          ISIGN=1
         ELSEIF(X .EQ. 0.0D0)THEN
          ISIGN=0
         ELSE
          ISIGN=-1
         ENDIF
         RETURN
         END

         SUBROUTINE TOLE(B,TOL,NEPS,EPS)
C DETERMINES THE TERMINATION CRITERION.
C  B    -- DOUBLE PRECISION.
C  NEPS -- INTEGER.
C  EPS  -- DOUBLE PRECISION.
C  TOL  -- DOUBLE PRECISION. OUTPUT AS THE TERMINATION CRITERION.
C           TOL =2*(2*EPS*|B| + 10D-{NEPS}),  IF NEPS IS NOT 1000;
C    AND    TOL =2*(2*EPS*|B|),               IF NEPS = 1000. 
         INTEGER NEPS,I
         DOUBLE PRECISION B,TOL,EPS
         IF(NEPS .EQ. 1000)THEN
          TOL=0.0D0
         ELSE
          TOL=1.0D0
          DO 10 I=1,NEPS
           TOL=TOL/10.0D0
 10       CONTINUE
         ENDIF
         TOL=TOL+2.0D0*DABS(B)*EPS
         TOL=2.0D0*TOL
         RETURN
         END

         SUBROUTINE NEWQUA(A,B,D,FA,FB,FD,C,K)
C USES K NEWTON STEPS TO APPROXIMATE THE ZERO IN (A,B) OF THE
C QUADRATIC POLYNOMIAL INTERPOLATING F(X) AT A, B, AND D. SAFEGUARD
C IS USED TO AVOID OVERFLOW.
C  A,B,D,FA,FB,FD -- DOUBLE PRECISION. D LIES OUTSIDE THE INTERVAL
C                    [A,B]. FA=F(A), FB=F(B), AND FD=F(D). F(A)F(B)<0.        
C  C              -- DOUBLE PRECISION. OUTPUT AS THE APPROXIMATE ZERO
C                    IN (A,B) OF THE QUADRATIC POLYNOMIAL.
C  K              -- INTEGER. INPUT INDICATING THE NUMBER OF NEWTON 
C                    STEPS TO TAKE.
         INTEGER K,IERROR,ISIGN,I
         DOUBLE PRECISION A,B,D,FA,FB,FD,C,A0,A1,A2,PC,PDC
         EXTERNAL ISIGN
C
C INITIALIZATION. FIND THE COEFFICIENTS OF THE QUADRATIC POLYNOMIAL.
C
         IERROR=0
         A0=FA
         A1=(FB-FA)/(B-A)
         A2=((FD-FB)/(D-B)-A1)/(D-A)
C
C SAFEGUARD TO AVOID OVERFLOW.
C
  10     IF((A2 .EQ. 0.0D0).OR.(IERROR .EQ. 1))THEN
          C=A-A0/A1
          RETURN
         ENDIF
C
C DETERMINE THE STARTING POINT OF NEWTON STEPS.
C
         IF(ISIGN(A2)*ISIGN(FA) .GT. 0)THEN
          C=A
         ELSE
          C=B
         ENDIF
C
C START THE SAFEGUARDED NEWTON STEPS.
C
         DO 20 I=1,K
         IF(IERROR .EQ. 0)THEN
          PC=A0+(A1+A2*(C-B))*(C-A)
          PDC=A1+A2*((2.0D0*C)-(A+B))
          IF(PDC .EQ. 0.0D0)THEN
           IERROR=1
          ELSE
           C=C-PC/PDC
          ENDIF
         ENDIF
  20     CONTINUE
         IF(IERROR .EQ. 1)GOTO 10
         RETURN
         END

         SUBROUTINE PZERO(A,B,D,E,FA,FB,FD,FE,C)
C USES CUBIC INVERSE INTERPOLATION OF F(X) AT A, B, D, AND E TO
C GET AN APPROXIMATE ROOT OF F(X). THIS PROCEDURE IS A SLIGHT
C MODIFICATION OF AITKEN-NEVILLE ALGORITHM FOR INTERPOLATION
C DESCRIBED BY STOER AND BULIRSCH IN "INTRO. TO NUMERICAL ANALYSIS"
C SPRINGER-VERLAG. NEW YORK (1980).
C  A,B,D,E,FA,FB,FD,FE -- DOUBLE PRECISION. D AND E LIE OUTSIDE
C                         THE INTERVAL [A,B]. FA=F(A), FB=F(B),
C                         FD=F(D), AND FE=F(E).
C  C                   -- DOUBLE PRECISION. OUTPUT OF THE SUBROUTINE.        
         DOUBLE PRECISION A,B,D,E,FA,FB,FD,FE,C,Q11,Q21,Q31,D21,D31         
         DOUBLE PRECISION Q22,Q32,D32,Q33
C
         Q11=(D-E)*FD/(FE-FD)
         Q21=(B-D)*FB/(FD-FB)
         Q31=(A-B)*FA/(FB-FA)
         D21=(B-D)*FD/(FD-FB)
         D31=(A-B)*FB/(FB-FA)
         Q22=(D21-Q11)*FB/(FE-FB)
         Q32=(D31-Q21)*FA/(FD-FA)
         D32=(D31-Q21)*FD/(FD-FA)
         Q33=(D32-Q22)*FA/(FE-FA)
C
C CALCULATE THE OUTPUT C.
C
         C=Q31+Q32+Q33
         C=A+C
         RETURN
         END

         SUBROUTINE RMP(REL)
C CALCULATES THE RELATIVE MACHINE PRECISION (RMP).
C  REL -- DOUBLE PRECISION. OUTPUT OF RMP.
C
         DOUBLE PRECISION REL,BETA,A,B
         BETA=2.0D0
         A=1.0D0
 10      B=1.0D0+A
         IF(B .GT. 1.0D0)THEN
          A=A/BETA
          GOTO 10
         ENDIF
         REL=A*BETA
         RETURN
         END
C*** exdrive.f
C ALGORITHM ENCLOFX
C
C G.E.ALEFELD, F.A.POTRA, AND YIXUN SHI --- JULY, 1994
C
C COMPANION PAPER: ENCLOSING ZEROS OF CONTINUOUS FUNCTIONS
C
C THIS IMPLEMENTS THE ALGORITHM 4.2 OF THE COMPANION PAPER.
C
C FORTRAN/77 PROGRAM.
C BUILT-IN FORTRAN/77 FUNCTION DABS(X) IS USED, WHERE
C    DABS(X)=|X| WITH DOUBLE PRECISION.
C ******************************************************************
         INTEGER NPROB,NEPS,ISIGN,INDI
         DOUBLE PRECISION A,B,EPS,ROOT
         EXTERNAL ISIGN
C THIS PROGRAM CALCULATES A ROOT OF A CONTINUOUS FUNCTION F(X)         
C IN AN INTERVAL I=[A, B] PROVIDED THAT F(A)F(B)<0. THE FUNCTION F(X)
C AND THE INITIAL INTERVAL [A, B] ARE TO BE SUPPLIED BY THE USER IN 
C THE SUBROUTINES "FUNC" AND "INIT". THE OUTPUT "ROOT" EITHER SATISFIES
C THAT F(ROOT)=0 OR IS AN APPROXIMATE SOLUTION OF THE EQUATION F(X)=0
C SUCH THAT "ROOT" IS INCLUDED IN AN INTERVAL [AC, BC] WITH
C      F(AC)F(BC)<0,
C AND
C      BC-AC <= TOL = 2*TOLE(AC,BC).
C IN THE ABOVE TERMINATION CRITERION,
C      TOLE(AC,BC)=2*EPS*|UC|+10D-{NEPS}, IF "NEPS" IS NOT 1000,        
C AND
C      TOLE(AC,BC)=2*EPS*|UC|,            IF "NEPS" IS SET AS 1000,
C WHERE "NEPS" AND "EPS" ARE THE INPUTS OF THE PROGRAM DETERMINED BY
C USER WITH "EPS" BEING LARGER THAN OR EQUAL TO THE RELATIVE MACHINE
C PRECISION (IT IS USUALLY CHOSEN AS THE RELATIVE MACHINE PRECISION),        
C AND "UC" IS EQUAL TO EITHER "AC" OR "BC" WITH THE CONDITION
C      |F(UC)| = MIN{ |F(AC)|, |F(BC)| }.
C
C INPUT OF THE PROGRAM:
C  INDI  -- INTEGER. INDICATES WHETHER THE USER WANTS TO USE THE
C           RELATIVE MACHINE PRECISION AS "EPS". IF INPUT 0, THEN
C           "EPS" WILL BE AUTOMATICALLY SET AS THE RELATIVE MACHINE
C           PRECISION, OTHERWISE THE USER WILL BE ASKED TO ENTER
C           "EPS";
C  EPS   -- DOUBLE PRECISION. USED IN TERMINATION CRITERION. LARGER
C           THAN OR EQUAL TO THE RELATIVE MACHINE PRECISION. IT IS
C           USUALLY CHOSEN AS THE RELATIVE MACHINE PRECISION;
C  NEPS  -- INTEGER. POSITIVE INTEGER USED IN DETERMINING THE
C           TERMINATION CRITERION;
C  NPROB -- INTEGER. POSITIVE INTEGER STANDING FOR "NUMBER OF PROBLEM",
C           INDICATING THE PROBLEM TO BE SOLVED.
C OUTPUT OF THE PROGRAM:
C  ROOT  -- DOUBLE PRECISION. EXACT OR APPROXIMATE SOLUTION OF THE
C           EQUATION F(X)=0.

C
C INDICATE WHETHER TO USE THE RELATIVE MACHINE PRECISION AS "EPS"
C OR TO ENTER A USER-CHOSEN "EPS". THE VALUE OF "EPS" SHOULD BE
C LARGER THAN OR EQUAL TO THE RELATIVE MACHINE PRECISION.
C
         PRINT*, 'ENTER 0 TO TAKE THE RELATIVE MACHINE PRECISION'
         PRINT*, 'AS "EPS", ENTER 1 TO TAKE OTHER VALUE FOR "EPS".'          
         READ*, INDI
         IF (INDI .EQ. 0)THEN
          CALL RMP(EPS)
         ELSE
          PRINT*, 'ENTER EPS.'
          READ*, EPS
         ENDIF
C
C ENTER THE INPUTS "NEPS", AND "NPROB".
C
 10      PRINT*, 'ENTER NEPS.'
         READ*, NEPS
 20      PRINT*, 'ENTER THE NUMBER OF PROBLEM.'
         READ*, NPROB
C
C CALL SUBROUTINE "INIT" TO GET THE INITIAL INTERVAL.
C 
         CALL INIT(NPROB, A, B)
C
C CALL SUBROUTINE "RROOT" TO HAVE THE PROBLEM SOLVED.
C
         CALL RROOT(NPROB,NEPS,EPS,A,B,ROOT)
         CONTINUE
C
C PRINT OUT THE RESULT ON THE SCREEN.
C
         PRINT*, 'THE NUMBER OF PROBLEM IS', NPROB
         PRINT*, 'THE ROOT IS', ROOT
         STOP
         END
         SUBROUTINE INIT(NPROB, A, B)
C  NPROB -- INTEGER. INDICATING THE NUMBER OF PROBLEM.
C  A, B  -- DOUBLE PRECISION. GIVE THE INTERVAL [A,B].
         INTEGER NPROB
         DOUBLE PRECISION A,B
         GOTO (10,20), NPROB
 10      CONTINUE
         A=0.0D0
         B=1.0D0
         RETURN
 20      CONTINUE
         A=1.0D0
         B=100.0D0
         RETURN
         END

         SUBROUTINE FUNC(NPROB, X, FX)
C GIVE THE CONTINUOUS FUNCTION F(X) OF THE PROBLEM.
C  NPROB -- INTEGER. INDICATING THE NUMBER OF PROBLEM.
C  X, FX -- DOUBLE PRECISION. FX=F(X).
         INTEGER NPROB
         DOUBLE PRECISION X, FX
         GOTO (10,20), NPROB
 10      CONTINUE
         FX=X**2-(1.0D0-X)**5
         RETURN
 20      CONTINUE
         FX=X**(1.0D0/5.0D0)-5.0D0**(1.0D0/5.0D0)
         RETURN
         END
C*** testdata
1  1
2  1
3  1
4  1
5  1
6  1
7  1
8  1
9  1
10  1
11  1
12  1
13  1
14  1
15  4
15  6
15  8
15  10
15  12
16  4
16  6
16  8
16  10
16  12
17  8
17  10
17  12
17  14
18  1
19  1
19  2
19  3
19  4
19  5
19  20
19  40
19  60
19  80
19  100
20  5
20  10
20  20
21  2
21  5
21  10
21  15
21  20
22  1
22  2
22  4
22  5
22  8
22  15
22  20
23  1
23  5
23  10
23  15
23  20
24  2
24  5
24  15
24  20
25  2
25  3
25  4
25  5
25  6
25  7
25  9
25  11
25  13
25  15
25  17
25  19
25  21
25  23
25  25
25  27
25  29
25  31
25  33
26  1
27  1
27  2
27  3
27  4
27  5
27  6
27  7
27  8
27  9
27  10
27  11
27  12
27  13
27  14
27  15
27  16
27  17
27  18
27  19
27  20
27  21
27  22
27  23
27  24
27  25
27  26
27  27
27  28
27  29
27  30
27  31
27  32
27  33
27  34
27  35
27  36
27  37
27  38
27  39
27  40
28  20
28  21
28  22
28  23
28  24
28  25
28  26
28  27
28  28
28  29
28  30
28  31
28  32
28  33
28  34
28  35
28  36
28  37
28  38
28  39
28  40
28  100
28  200
28  300
28  400
28  500
28  600
28  700
28  800
28  900
28  1000
C*** testout
PROBLEM NUMBER   1 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =    1.8954942670340
 
PROBLEM NUMBER   2 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =    3.0229153472731
 
PROBLEM NUMBER   3 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =    6.6837535608081
 
PROBLEM NUMBER   4 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =    11.238701655002
 
PROBLEM NUMBER   5 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =    19.676000080623
 
PROBLEM NUMBER   6 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =    29.828227326505
 
PROBLEM NUMBER   7 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =    41.906116195289
 
PROBLEM NUMBER   8 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =    55.953595800143
 
PROBLEM NUMBER   9 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =    71.985665586588
 
PROBLEM NUMBER   10 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =    90.008868539167
 
PROBLEM NUMBER   11 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =    110.02653274833
 
PROBLEM NUMBER   12 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =  0.
 
PROBLEM NUMBER   13 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =  0.
 
PROBLEM NUMBER   14 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =  0.
 
PROBLEM NUMBER   15 WITH PROBLEM PARAMETER   4
COMPUTED ROOT =   0.66874030497642
 
PROBLEM NUMBER   15 WITH PROBLEM PARAMETER   6
COMPUTED ROOT =   0.76472449133173
 
PROBLEM NUMBER   15 WITH PROBLEM PARAMETER   8
COMPUTED ROOT =   0.81776543395794
 
PROBLEM NUMBER   15 WITH PROBLEM PARAMETER   10
COMPUTED ROOT =   0.85133992252078
 
PROBLEM NUMBER   15 WITH PROBLEM PARAMETER   12
COMPUTED ROOT =   0.87448527222117
 
PROBLEM NUMBER   16 WITH PROBLEM PARAMETER   4
COMPUTED ROOT =    1.0000000000000
 
PROBLEM NUMBER   16 WITH PROBLEM PARAMETER   6
COMPUTED ROOT =    1.0000000000000
 
PROBLEM NUMBER   16 WITH PROBLEM PARAMETER   8
COMPUTED ROOT =    1.0000000000000
 
PROBLEM NUMBER   16 WITH PROBLEM PARAMETER   10
COMPUTED ROOT =    1.0000000000000
 
PROBLEM NUMBER   16 WITH PROBLEM PARAMETER   12
COMPUTED ROOT =    1.0000000000000
 
PROBLEM NUMBER   17 WITH PROBLEM PARAMETER   8
COMPUTED ROOT =    1.0000000000000
 
PROBLEM NUMBER   17 WITH PROBLEM PARAMETER   10
COMPUTED ROOT =    1.0000000000000
 
PROBLEM NUMBER   17 WITH PROBLEM PARAMETER   12
COMPUTED ROOT =    1.0000000000000
 
PROBLEM NUMBER   17 WITH PROBLEM PARAMETER   14
COMPUTED ROOT =    1.0000000000000
 
PROBLEM NUMBER   18 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =   0.52359877559830
 
PROBLEM NUMBER   19 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =   0.42247770964124
 
PROBLEM NUMBER   19 WITH PROBLEM PARAMETER   2
COMPUTED ROOT =   0.30669941048320
 
PROBLEM NUMBER   19 WITH PROBLEM PARAMETER   3
COMPUTED ROOT =   0.22370545765466
 
PROBLEM NUMBER   19 WITH PROBLEM PARAMETER   4
COMPUTED ROOT =   0.17171914751951
 
PROBLEM NUMBER   19 WITH PROBLEM PARAMETER   5
COMPUTED ROOT =   0.13825715505682
 
PROBLEM NUMBER   19 WITH PROBLEM PARAMETER   20
COMPUTED ROOT =    3.4657359020854D-02
 
PROBLEM NUMBER   19 WITH PROBLEM PARAMETER   40
COMPUTED ROOT =    1.7328679513999D-02
 
PROBLEM NUMBER   19 WITH PROBLEM PARAMETER   60
COMPUTED ROOT =    1.1552453009332D-02
 
PROBLEM NUMBER   19 WITH PROBLEM PARAMETER   80
COMPUTED ROOT =    8.6643397569993D-03
 
PROBLEM NUMBER   19 WITH PROBLEM PARAMETER   100
COMPUTED ROOT =    6.9314718055995D-03
 
PROBLEM NUMBER   20 WITH PROBLEM PARAMETER   5
COMPUTED ROOT =    3.8402551840622D-02
 
PROBLEM NUMBER   20 WITH PROBLEM PARAMETER   10
COMPUTED ROOT =    9.9000099980005D-03
 
PROBLEM NUMBER   20 WITH PROBLEM PARAMETER   20
COMPUTED ROOT =    2.4937500390620D-03
 
PROBLEM NUMBER   21 WITH PROBLEM PARAMETER   2
COMPUTED ROOT =   0.50000000000000
 
PROBLEM NUMBER   21 WITH PROBLEM PARAMETER   5
COMPUTED ROOT =   0.34595481584824
 
PROBLEM NUMBER   21 WITH PROBLEM PARAMETER   10
COMPUTED ROOT =   0.24512233375331
 
PROBLEM NUMBER   21 WITH PROBLEM PARAMETER   15
COMPUTED ROOT =   0.19554762353657
 
PROBLEM NUMBER   21 WITH PROBLEM PARAMETER   20
COMPUTED ROOT =   0.16492095727644
 
PROBLEM NUMBER   22 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =   0.27550804099948
 
PROBLEM NUMBER   22 WITH PROBLEM PARAMETER   2
COMPUTED ROOT =   0.13775402049974
 
PROBLEM NUMBER   22 WITH PROBLEM PARAMETER   4
COMPUTED ROOT =    1.0305283778156D-02
 
PROBLEM NUMBER   22 WITH PROBLEM PARAMETER   5
COMPUTED ROOT =    3.6171081789041D-03
 
PROBLEM NUMBER   22 WITH PROBLEM PARAMETER   8
COMPUTED ROOT =    4.1087291849640D-04
 
PROBLEM NUMBER   22 WITH PROBLEM PARAMETER   15
COMPUTED ROOT =    2.5989575892908D-05
 
PROBLEM NUMBER   22 WITH PROBLEM PARAMETER   20
COMPUTED ROOT =    7.6685951221853D-06
 
PROBLEM NUMBER   23 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =   0.40105813754155
 
PROBLEM NUMBER   23 WITH PROBLEM PARAMETER   5
COMPUTED ROOT =   0.51615351875793
 
PROBLEM NUMBER   23 WITH PROBLEM PARAMETER   10
COMPUTED ROOT =   0.53952222690842
 
PROBLEM NUMBER   23 WITH PROBLEM PARAMETER   15
COMPUTED ROOT =   0.54818229434066
 
PROBLEM NUMBER   23 WITH PROBLEM PARAMETER   20
COMPUTED ROOT =   0.55270466667849
 
PROBLEM NUMBER   24 WITH PROBLEM PARAMETER   2
COMPUTED ROOT =   0.50000000000000
 
PROBLEM NUMBER   24 WITH PROBLEM PARAMETER   5
COMPUTED ROOT =   0.20000000000000
 
PROBLEM NUMBER   24 WITH PROBLEM PARAMETER   15
COMPUTED ROOT =    6.6666666666667D-02
 
PROBLEM NUMBER   24 WITH PROBLEM PARAMETER   20
COMPUTED ROOT =    5.0000000000000D-02
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   2
COMPUTED ROOT =    2.0000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   3
COMPUTED ROOT =    3.0000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   4
COMPUTED ROOT =    4.0000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   5
COMPUTED ROOT =    5.0000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   6
COMPUTED ROOT =    6.0000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   7
COMPUTED ROOT =    7.0000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   9
COMPUTED ROOT =    9.0000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   11
COMPUTED ROOT =    11.000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   13
COMPUTED ROOT =    13.000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   15
COMPUTED ROOT =    15.000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   17
COMPUTED ROOT =    17.000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   19
COMPUTED ROOT =    19.000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   21
COMPUTED ROOT =    21.000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   23
COMPUTED ROOT =    23.000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   25
COMPUTED ROOT =    25.000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   27
COMPUTED ROOT =    27.000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   29
COMPUTED ROOT =    29.000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   31
COMPUTED ROOT =    31.000000000000
 
PROBLEM NUMBER   25 WITH PROBLEM PARAMETER   33
COMPUTED ROOT =    33.000000000000
 
PROBLEM NUMBER   26 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =    2.2317679157465D-02
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   1
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   2
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   3
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   4
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   5
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   6
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   7
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   8
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   9
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   10
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   11
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   12
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   13
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   14
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   15
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   16
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   17
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   18
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   19
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   20
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   21
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   22
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   23
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   24
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   25
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   26
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   27
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   28
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   29
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   30
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   31
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   32
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   33
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   34
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   35
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   36
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   37
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   38
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   39
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   27 WITH PROBLEM PARAMETER   40
COMPUTED ROOT =   0.62380651896161
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   20
COMPUTED ROOT =    5.9051305594220D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   21
COMPUTED ROOT =    5.6367155339937D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   22
COMPUTED ROOT =    5.3916409455592D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   23
COMPUTED ROOT =    5.1669892394942D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   24
COMPUTED ROOT =    4.9603096699145D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   25
COMPUTED ROOT =    4.7695285287639D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   26
COMPUTED ROOT =    4.5928793239949D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   27
COMPUTED ROOT =    4.4288479195665D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   28
COMPUTED ROOT =    4.2761290257883D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   29
COMPUTED ROOT =    4.1335913915954D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   30
COMPUTED ROOT =    4.0002497338020D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   31
COMPUTED ROOT =    3.8752419296207D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   32
COMPUTED ROOT =    3.7578103559958D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   33
COMPUTED ROOT =    3.6472865219959D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   34
COMPUTED ROOT =    3.5430783356532D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   35
COMPUTED ROOT =    3.4446594929961D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   36
COMPUTED ROOT =    3.3515605877800D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   37
COMPUTED ROOT =    3.2633616249437D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   38
COMPUTED ROOT =    3.1796856858426D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   39
COMPUTED ROOT =    3.1001935436965D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   40
COMPUTED ROOT =    3.0245790670210D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   100
COMPUTED ROOT =    1.2277994232462D-05
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   200
COMPUTED ROOT =    6.1695393904409D-06
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   300
COMPUTED ROOT =    4.1198585298293D-06
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   400
COMPUTED ROOT =    3.0924623877272D-06
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   500
COMPUTED ROOT =    2.4752044261050D-06
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   600
COMPUTED ROOT =    2.0633567678513D-06
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   700
COMPUTED ROOT =    1.7690120078154D-06
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   800
COMPUTED ROOT =    1.5481615698859D-06
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   900
COMPUTED ROOT =    1.3763345366022D-06
 
PROBLEM NUMBER   28 WITH PROBLEM PARAMETER   1000
COMPUTED ROOT =    1.2388385788997D-06