      PROGRAM NEWTON
************************************************************************
* Program to find an approximate complex root of a polynomial P using  *
* Newton's method.  Identifiers used are:                              *
*     DEGLIM : (parameter) limit on the polynomial's degree            *
*     DEGREE : actual degree of the polynomial                         *
*     COEFF  : array of coefficients of the polynomial                 *
*     XOLD   : previous approximation (initially the first one)        *
*     PPOLD  : value of the derivative of P at XOLD                    *
*     XNEW   : the new approximation                                   *
*     PVALUE : value of P at an approximation                          *
*     EPSIL  : repetition stops when ABS(PVALUE) is less than EPSIL    *
*     NUMITS : limit on number of iterations                           *
*     N      : number of iterations                                    *
*     I      : index                                                   *
*                                                                      *
* Input:  An initial approximation                                     *
* Output: A sequence of approximations to the root or an error message *
*         indicating that the method fails                             *
************************************************************************

      INTEGER DEGLIM, DEGREE, NUMITS, I, N
      PARAMETER (DEGLIM = 20)
      REAL EPSIL
      COMPLEX COEFF(0:DEGLIM), XOLD, PPOLD, XNEW,  P, PPRIME, PVALUE

* Get degree and coefficients of the polynomial
      PRINT *, 'ENTER THE DEGREE OF THE POLYNOMIAL:'
      READ *, DEGREE
      PRINT *, 'ENTER THE COMPLEX COEFFICIENTS OF THE POLYONIMAL'
      PRINT *, '(IN DESCENDING ORDER OF THE POWERS OF X):'
      READ *, (COEFF(I), I = DEGREE, 0, -1)

* Get termination values EPSIL & NUMITS and initial approximation
      PRINT *, 'ENTER EPSILON, LIMIT ON # OF ITERATIONS, AND'
      PRINT *, 'THE INITIAL APPROXIMATION:'
      READ *, EPSIL, NUMITS, XOLD
      PVALUE = P(COEFF, DEGREE, XOLD)
      N = 0
      PRINT *, '  N             X(N)           ABS(P(X(N)))'
      PRINT *, '============================================'
      PRINT 10, 0, XOLD, ABS(PVALUE)
 10   FORMAT(1X, I3, F10.5, ' +', F10.5, 'I', EG14.5)



* Iterate using Newton's method while ABS(PVALUE) is greater
* than or equal to EPSIL and N has not reached NUMITS.
* Terminate if the derivative is 0 at some approximation.

20    IF ((ABS(PVALUE) .GE. EPSIL) .AND. (N .LT. NUMITS)) THEN
          N = N + 1
          PPOLD = PPRIME(COEFF, DEGREE, XOLD)
          IF (PPOLD .EQ. (0,0)) THEN
              PRINT *, 'NEWTON''S METHOD FAILS -- DERIVATIVE = 0'
*             Force repetition to terminate
              N = NUMITS
          ELSE
              XNEW = XOLD - (PVALUE / PPOLD)
              PVALUE = P(COEFF, DEGREE, XNEW)
              PRINT 10, N, XNEW, ABS(PVALUE)
              XOLD = XNEW
          END IF
      GO TO 20
      END IF
      END


**P(A,N,X)**************************************************************
* For the polynomial with N coefficients stored in array A, compute    *
* its value at X, using Horner's method.  Local variable used:         *
*     I : index                                                        *
************************************************************************

      FUNCTION P(A, N, X)

      INTEGER N, I
      COMPLEX A(0:N), X, P

      P = A(N)
      DO 10 I = N - 1, 0, -1
         P = A(I)  +  X * P
10    CONTINUE
      END


**PPRIME(A,N,X)*********************************************************
* For the polynomial with N coefficients stored in array A, compute    *
* the value of its derivative at X, using Horner's method.  Local      *
* variable used:                                                       *
*     I : index                                                        *
************************************************************************

      FUNCTION PPRIME(A, N, X)

      INTEGER N, I
      COMPLEX A(0:N), X, PPRIME

      PPRIME = N * A(N)
      DO 10 I = N - 1, 1, -1
         PPRIME = I * A(I)  +  X * PPRIME
10    CONTINUE
      END   

