SUBROUTINE FFUNP(N,NUMT,MMAXT,KDEG,COEF,CL,X,
     $ XX,TRM,DTRM,CLX,DXNP1,F,DF)
C
C FFUNP  EVALUATES THE SYSTEM "F(X)=0" AND ITS PARTIAL
C DERIVATIVES, USING THE "TABLEAU" INPUT: N,NUMT,KDEG,COEF.
C
C FFUNP  CAN BE MADE MORE EFFICIENT BY CUSTOMIZING IT TO
C PARTICULAR SYSTEM TYPES.  FOR EXAMPLE,
C IF X(1)**2 AND X(1)**3 ARE USED IN SEVERAL
C EQUATIONS, THE CURRENT  FFUNP  RECOMPUTES BOTH OF THESE FOR
C EACH EQUATION.  BUT (OF COURSE) WE CAN COMPUTE
C X1SQ=X(1)**2 AND X1CU=XSQ(1)*X(1), AND
C USE THESE IN EACH OF THE EQUATIONS.
C
C THE PART OF THE CODE BELOW LABELED "BLOCK A" CAN BE
C CUSTOMIZED IN THIS WAY.   (THE CODE OUTSIDE OF
C BLOCK A CONCERNS THE PROJECTIVE TRANSFORMATION AND NEED NOT
C BE CHANGED.)  HOWEVER, BLOCK A REQUIRES THE HOMOGENEOUS FORM
C OF THE POLYNOMIALS RATHER THAN THE STANDARD FORM.  FURTHER,
C THE PARTIAL DERIVATIVES WITH RESPECT TO ALL N+1 PROJECTIVE
C VARIABLES MUST BE COMPUTED.  MORE EXPLICITLY,
C THE ORIGINAL SYSTEM, F(X)=0, IS GIVEN IN "NON-HOMOGENEOUS FORM" AS
C DESCRIBED IN SUBROUTINE POLSYS.  F(X)  IS
C REPRESENTED IN "HOMOGENEOUS FORM" AS FOLLOWS:
C
C              NUMT(J)
C
C    F(J) =     SUM   TRM(J,K)
C
C               K=1
C
C WHERE  TRM(J,K)=COEF(J,K) * XX(J,1,K)*XX(J,2,K)* ... *XX(J,N+1,K)
C
C WITH XX(J,L,K) = X(L)**KDEG(J,L,K) FOR J=1 TO N, L=1 TO N, AND
C K=1 TO NUMT(J), AND WITH XX(J,N+1,K) = XNP1**KDEG(J,N+1,K) FOR J=1 TO
C N AND K=1 TO NUMT(J), WHERE  XNP1  IS THE "HOMOGENEOUS COORDINATE,"
C KDEG(J,N+1,K)=IDEG(J)-(KDEG(J,1,K)+ ... + KDEG(J,N,K)),
C AND IDEG(J) THE DEGREE OF THE J-TH EQUATION.   XNP1  IS GENERATED
C FROM  X  AND  CL  BEFORE BLOCK A.
C
C IN THIS DISCUSSION WE HAVE OMITTED, FOR SIMPLICITY OF
C EXPOSITION, THE LEADING INDEX, WHICH DIFFERENTIATES THE
C REAL AND IMAGINARY PARTS.  HOWEVER, THIS INDEX MUST NOT BE
C OMITTED IN THE CODE.
C
C WE COMPLETE THE EXPOSITION OF "REPLACING BLOCK A WITH MORE EFFICIENT
C CODE" WITH AN EXPLICIT EXAMPLE.  FIRST, THE SYSTEM IS DESCRIBED.
C THEN THE CODE THAT SHOULD BE USED IS GIVEN (COMMENTED OUT).
C IN TESTS  POLSYS  WITH THE MORE EFFICIENT  FFUNP  RAN ABOUT TWICE AS
C FAST AS WITH THE GENERIC  FFUNP .
C
C HERE IS THE SYSTEM TO BE SOLVED:
C
C     F(1) = COEF(1,1) * X(1)**4
C    &     + COEF(1,2) * X(1)**3 * X(2)
C    &     + COEF(1,3) * X(1)**3
C    &     + COEF(1,4) * X(1)
C    &     + COEF(1,5)
C     F(2) = COEF(2,1) * X(1)     * X(2)**2
C    &     + COEF(2,2)              X(2)**2
C    &     + COEF(2,3)
C
C THE REPLACEMENT CODE REQUIRES THE FOLLOWING DECLARATIONS:
C     DOUBLE PRECISION X1SQ,X1CU,X2SQ,X3SQ,X3CU,
C    &  TEMPA,TEMPB,TEMPC,TEMPD,TEMPE,TEMPF
C     DIMENSION X1SQ(2),X1CU(2),X2SQ(2),X3SQ(2),X3CU(2),
C    &  TEMPA(2),TEMPB(2),TEMPC(2),TEMPD(2),TEMPE(2),TEMPF(2)
C
C HERE IS CODE TO REPLACE BLOCK A:
C
C******************  BEGIN BLOCK A  *******************
C
C     CALL MULP(X(1,1),X(1,1),X1SQ)
C     CALL MULP(X1SQ  ,X(1,1),X1CU)
C     CALL MULP(X(1,2),X(1,2),X2SQ)
C     CALL MULP(XNP1,  XNP1,  X3SQ)
C     CALL MULP(X3SQ  ,XNP1,  X3CU)
C
C     DO 1 I=1,2
C       TEMPA(I)=   COEF(1,1) * X(I,1)
C    &            + COEF(1,2) * X(I,2)
C    &            + COEF(1,3) * XNP1(I)
C       TEMPB(I)=   COEF(1,4) * X(I,1)
C    &            + COEF(1,5) * XNP1(I)
C 1   CONTINUE
C
C     CALL MULP(X1SQ,  TEMPA,TEMPC)
C     CALL MULP(X(1,1),TEMPC,TEMPD)
C     CALL MULP(X3SQ,  TEMPB,TEMPE)
C     CALL MULP(XNP1,  TEMPE,TEMPF)
C
C     DO 2 I=1,2
C       F(I,1)=TEMPD(I) + TEMPF(I)
C       DF(I,1,1)= 3. *TEMPC(I) + COEF(1,1)*X1CU(I) + COEF(1,4)*X3CU(I)
C       DF(I,1,2)= COEF(1,2) * X1CU(I)
C       DF(I,1,3)= COEF(1,3)*X1CU(I) + 3. *TEMPE(I) + COEF(1,5)*X3CU(I)
C
C       TEMPA(I) = COEF(2,1) * X(I,1) + COEF(2,2) * XNP1(I)
C  2  CONTINUE
C
C     CALL MULP(TEMPA,X(1,2),TEMPB)
C     CALL MULP(TEMPB,X(1,2),TEMPC)
C
C     DO 3 I=1,2
C       F(I,2) = TEMPC(I) + COEF(2,3) * X3CU(I)
C       DF(I,2,1) = COEF(2,1) * X2SQ(I)
C       DF(I,2,2) = 2. * TEMPB(I)
C       DF(I,2,3) = COEF(2,2) * X2SQ(I) + COEF(2,3) * 3. * X3SQ(I)
C  3  CONTINUE
C******************  END OF BLOCK A  *******************
C
C ON INPUT:
C
C N  IS THE NUMBER OF EQUATIONS AND VARIABLES.
C
C NUMT(J)  IS THE NUMBER OF TERMS IN THE JTH EQUATION.
C
C MMAXT  IS AN UPPER BOUND ON NUMT(J) FOR J=1 TO N.
C
C KDEG(J,L,K)  IS THE DEGREE OF THE L-TH VARIABLE IN THE K-TH TERM
C   OF THE J-TH EQUATION.
C
C COEF(J,K)  IS THE K-TH COEFFICIENT OF THE J-TH EQUATION.
C
C CL  IS USED TO DEFINE THE PROJECTIVE TRANSFORMATION.  IF
C   THE PROJECTIVE TRANSFORMATION IS NOT SPECIFIED, THEN  CL
C   CONTAINS DUMMY VALUES.
C
C X(1,J), X(2,J)  ARE THE REAL AND IMAGINARY PARTS RESPECTIVELY OF
C   THE J-TH INDEPENDENT VARIABLE.
C
C XX, TRM, DTRM, CLX, DXNP1  ARE WORKSPACE VARIABLES.
C
C ON OUTPUT:
C
C F(1,J), F(2,J)  ARE THE REAL AND IMAGINARY PARTS RESPECTIVELY OF
C   THE J-TH EQUATION.
C
C DF(1,J,K), DF(2,J,K)  ARE THE REAL AND IMAGINARY PARTS RESPECTIVELY
C   OF THE K-TH PARTIAL DERIVATIVE OF THE J-TH EQUATION.
C
C
C VARIABLES: XNP1,TEMP1,TEMP2.
C
C NOTE:  XNP1(1), XNP1(2)  ARE THE REAL AND IMAGINARY PARTS,
C   RESPECTIVELY, OF THE PROJECTIVE VARIABLE.  XNP1  IS UNITY
C   IF THE PROJECTIVE TRANSFORMATION IS NOT SPECIFIED.
C
C  SUBROUTINES: MULP,POWP,DIVP.
C
C
C DECLARATION OF INPUT AND OUTPUT:
      INTEGER N,NUMT,MMAXT,KDEG
      DOUBLE PRECISION COEF,CL,X,XX,TRM,DTRM,CLX,DXNP1,F,DF
      DIMENSION NUMT(N),KDEG(N,N+1,MMAXT),
     $  COEF(N,MMAXT),CL(2,N+1),X(2,N),
     $  XX(2,N,N+1,MMAXT),TRM(2,N,MMAXT),DTRM(2,N,N+1,MMAXT),
     $  CLX(2,N),DXNP1(2,N),F(2,N),DF(2,N,N+1)
C
C DECLARATION OF VARIABLES:
      INTEGER I,IERR,J,K,L,M,NNNN,NP1
      DOUBLE PRECISION TEMP1,TEMP2,XNP1
      DIMENSION TEMP1(2),TEMP2(2),XNP1(2)
C
      NP1=N+1
C
C GENERATE XNP1, THE PROJECTIVE COORDINATE, AND ITS DERIVATIVES.
      DO 40 J=1,N
         CALL MULP(CL(1,J),X(1,J),CLX(1,J))
 40   CONTINUE
C
      DO 60 I=1,2
          XNP1(I)=CL(I,NP1)
          DO 50 J=1,N
              XNP1(I) = XNP1(I) + CLX(I,J)
              DXNP1(I,J)=CL(I,J)
 50       CONTINUE
 60   CONTINUE
C
C******************  BEGIN BLOCK A  *******************
C
C "BLOCK A" TAKES  X  AND  XNP1  AS INPUT AND RETURNS  F
C AND  DF  AS OUTPUT.   F  IS THE HOMOGENEOUS FORM OF THE
C ORIGINAL  F, AND  DF  CONSISTS OF THE PARTIAL
C DERIVATIVES OF THE HOMOGENEOUS FORM OF  F  WITH RESPECT
C TO THE N+1 VARIABLES X(1), ... ,X(N), XNP1.
C
C BEGIN "COMPUTE F"
C
      DO 100 J=1,N
        DO 100 K=1,NUMT(J)
          CALL POWP(KDEG(J,NP1,K),XNP1, XX(1,J,NP1,K))
          DO 100 L=1,N
            CALL POWP(KDEG(J, L,K),X(1,L),XX(1,J,  L,K))
 100  CONTINUE
      DO 200 J=1,N
        DO 200 K=1,NUMT(J)
          TRM(1,J,K)=COEF(J,K)
          TRM(2,J,K)=0.0
          DO 120 L=1,NP1
            CALL MULP(XX(1,J,L,K), TRM(1,J,K),TEMP1)
            TRM(1,J,K  )=TEMP1(1)
            TRM(2,J,K  )=TEMP1(2)
 120      CONTINUE
 200  CONTINUE
      DO 300 J=1,N
          F(1,J)=0.0
          F(2,J)=0.0
          DO 220 I=1,2
          DO 220 K=1,NUMT(J)
              F(I,J)= F(I,J) + TRM(I,J,K)
 220      CONTINUE
 300  CONTINUE
C
C END OF "COMPUTE F"
C
C BEGIN "COMPUTE DF"
C
      DO 400 J=1,N
        DO 400 K=1,NUMT(J)
        DO 400 M=1,NP1
C
C IF TERM DOES NOT INCLUDE X(M), SET PARTIAL DERIVATIVE OF TERM
C   EQUAL TO ZERO.
          IF(KDEG(J,M,K) .EQ. 0) THEN
            DTRM(1,J,M,K)=0.0
            DTRM(2,J,M,K)=0.0
          ELSE
C
C IF TERM DOES INCLUDE X(M), TRY COMPUTING THE PARTIAL BY DIVIDING
C   THE TERM BY X(M).
            IF(M.LE.N) CALL DIVP(TRM(1,J,K),X(1,M),DTRM(1,J,M,K),IERR)
            IF(M.EQ.NP1) CALL DIVP(TRM(1,J,K),XNP1,DTRM(1,J,M,K),IERR)
            IF (IERR .EQ. 0) THEN
              DTRM(1,J,M,K)=KDEG(J,M,K)*DTRM(1,J,M,K)
              DTRM(2,J,M,K)=KDEG(J,M,K)*DTRM(2,J,M,K)
            ELSE
C
C IF DIVISION WOULD CAUSE OVERFLOW, GENERATE THE PARTIAL BY
C   THE POLYNOMIAL FORMULA.
              DTRM(1,J,M,K)=COEF(J,K)
              DTRM(2,J,M,K)=0.0
              DO 320 L=1,NP1
                IF (L .EQ. M) GOTO 320
                CALL MULP(XX(1,J,L,K),DTRM(1,J,M,K),TEMP1)
                DTRM(1,J,M,K)=TEMP1(1)
                DTRM(2,J,M,K)=TEMP1(2)
 320          CONTINUE
              NNNN=KDEG(J,M,K)-1
              IF (M .LE. N) CALL POWP(NNNN,X(1,M),TEMP2)
              IF (M .EQ. NP1) CALL POWP(NNNN,XNP1 ,TEMP2)
              CALL MULP(TEMP2,TEMP1,DTRM(1,J,M,K))
              DTRM(1,J,M,K)=KDEG(J,M,K)*DTRM(1,J,M,K)
              DTRM(2,J,M,K)=KDEG(J,M,K)*DTRM(2,J,M,K)
            END IF
          END IF
 400  CONTINUE
      DO 600 J=1,N
      DO 600 M=1,NP1
          DF(1,J,M)=0.0
          DF(2,J,M)=0.0
          DO 420 I=1,2
            DO 420 K=1,NUMT(J)
              DF(I,J,M)= DF(I,J,M) + DTRM(I,J,M,K)
 420      CONTINUE
 600  CONTINUE
C
C END OF "COMPUTE DF"
C*******************  END BLOCK A  ********************
C
C CONVERT  DF  TO BE PARTIALS WITH RESPECT TO  X(1), ... ,X(N),
C BY APPLYING THE CHAIN RULE WITH  XNP1  CONSIDERED A FUNCTION OF
C OF  X(1), ... ,X(N).
C
      DO 700 J=1,N
        DO 700 K=1,N
          CALL MULP(DF(1,J,NP1),DXNP1(1,K),TEMP1)
          DO 700 I=1,2
            DF(I,J,K)=DF(I,J,K)+TEMP1(I)
 700  CONTINUE
      RETURN
      END

.