      SUBROUTINE ADCOMP(F)
C
C        ===============================================================
C
C **  THIS PROGRAM COMPUTES THE PIECEWISE POLYNOMIAL APPROXIMATION ON
C     THE CURRENT INTERVAL. IT ALSO ESTIMATES THE ERROR
C
      DOUBLE PRECISION A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP,
     * DSCTOL, ERROR, ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM,
     * XBREAK, XDD, XINTRP, XLEFT, XRIGHT, BUFFER
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50)
      DIMENSION DDTEMP(20,20), FACTOR(20), FINTRP(18), FLEFT(10),
     * FRIGHT(10), XDD(20), XINTRP(18)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, RIGHT, RIGHTX, SMOOTH
      LOGICAL DISCRD
      DOUBLE PRECISION ABSC, AJ, DX, FDERVL, FDERVR, FDUMB, R, WGTS
      DIMENSION ABSC(4), WGTS(4), FDERVL(9), FDERVR(9), FDUMB(9)
      DOUBLE PRECISION ERRINT, F, POLYDD
      EXTERNAL F, POLYDD
C
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR, MAXSTK,
     * NPAR, NSTACK, RIGHT, DISCRD, BUFFER
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
      EQUIVALENCE (FLEFT(2),FDERVL(1)), (FRIGHT(2),FDERVR(1))
      DATA R/1.5D0/
C
C     THIRTY DIGIT VALUES FOR THE GAUSS INTEGRATION CONSTANTS
C                    .861136311594052575223946488893D0
C                    .339981043584856264802665759103D0
C                    .347854845137453857373063949222D0
C                    .652145154862546142626936050778D0
C
C     ***** THE ABSISSAE AND WEIGHTS ARE GIVEN BELOW TO 30 DIGITS.
C           THE PARAMETER EPS0 IN ADSET SPECIFIES THE ACCURACY OF
C           THESE CONSTANTS. IF THE ACCURACY IS CHANGED TO K DECIMAL
C           DIGITS THEN SET EPS0 = 10**(-K).
C
      DATA ABSC(1) /-.861136311594052575223946488893D0 /
      DATA ABSC(2) /-.339981043584856264802665759103D0 /
      DATA ABSC(3) / .339981043584856264802665759103D0 /
      DATA ABSC(4) / .861136311594052575223946488893D0 /
      DATA WGTS(1) / .347854845137453857373063949222D0 /
      DATA WGTS(2) / .652145154862546142626936050778D0 /
      DATA WGTS(3) / .652145154862546142626936050778D0 /
      DATA WGTS(4) / .347854845137453857373063949222D0 /
C
C             COMPUTE INTERPOLATION INFORMATION
      NINTRP = DEGREE - 2*SMOOTH - 1
C
C          INCREASE NUMBER OF INTERPOLATION POINTS IF BREAK POINTS ARE
C          SPECIFIED WITH FEWER DERIVATIVES THAN SMOOTH
      IF (BREAK.EQ.LEFT .OR. BREAK.EQ.RIGHT) NINTRP = NINTRP + SMOOTH -
     * DBREAK(IBREAK)
      IF (BREAK.EQ.BOTH) NINTRP = NINTRP + 2*SMOOTH - DBREAK(IBREAK) -
     * DBREAK(IBREAK+1)
      IF (NINTRP.EQ.0) GO TO 20
C
C             GENERATE EQUAL SPACED INTERPOLATION POINTS
      AJ = NINTRP+1
      DX = (XRIGHT(NSTACK)-XLEFT(NSTACK))/AJ
      DO 10 J=1,NINTRP
        AJ = J
        XINTRP(J) = XLEFT(NSTACK) + AJ*DX
   10 CONTINUE
C
C             GET LEFT AND RIGHT F-VALUES, PUT F-VALUE IN FIRST ELEMENT
C             OF ARRAYS FLEFT AND FRIGHT.  GET DERIVATIVES BACK AS
C             OTHER ELEMENTS VIA THE SUBARRAYS FDERVL AND FDERVR.
   20 FLEFT(1) = F(XLEFT(NSTACK),FDERVL)
      FRIGHT(1) = F(XRIGHT(NSTACK),FDERVR)
      LEFTX = SMOOTH + 1
      RIGHTX = LEFTX
C            GET F-VALUES AT OTHER INTERPOLATION POINTS, IF ANY
      IF (NINTRP.EQ.0) GO TO 40
      DO 30 J=1,NINTRP
        FINTRP(J) = F(XINTRP(J),FDUMB)
   30 CONTINUE
C
C          CHECK FOR BREAK POINTS, MODIFY VALUES IF NECESSARY
   40 CONTINUE
      IF (BREAK.NE.LEFT) GO TO 50
      LEFTX = DBREAK(IBREAK) + 1
      FLEFT(LEFTX) = BRIGHT(IBREAK)
   50 IF (BREAK.NE.RIGHT) GO TO 60
      RIGHTX = DBREAK(IBREAK) + 1
      FRIGHT(RIGHTX) = BLEFT(IBREAK)
   60 IF (BREAK.NE.BOTH) GO TO 70
      LEFTX = DBREAK(IBREAK) + 1
      RIGHTX = DBREAK(IBREAK+1) + 1
      FLEFT(LEFTX) = BRIGHT(IBREAK)
      FRIGHT(RIGHTX) = BLEFT(IBREAK+1)
   70 CONTINUE
C
C           COMPUTE DIVIDED DIFFERENCES, NEWTON FORM OF POLYNOMIAL
      CALL NEWTON(LEFTX, RIGHTX, NINTRP)
C
C         COMPUTE NORM OF ERROR OF THIS APPROMIMATION USING FOUR PTS
C         ADD 50 PERCENT FUDGE FACTOR
      ERRORI = ERRINT(F,POLYDD,XLEFT(NSTACK),XRIGHT(NSTACK),ABSC,WGTS)
      ERRORI = R*ERRORI
      RETURN
      END
