!*************************************************************************
!*     EVALUATE A POLYNOMIAL AND ITS DERIVATIVES BY HORNER'S METHOD      *
!* --------------------------------------------------------------------- *
!* SAMPLE RUN:                                                           *
!*                                                                       *
!*    EVALUATE A POLYNOMIAL AND ITS DERIVATIVES BY HORNER'S METHOD       *
!*                                                                       *
!*            Example:  P(X) = X^4 + 2X^3 +3X^2 + 4X + 5                 *
!*                                                                       *
!*      X          P(X)          P'(X)         P"(X)         P"'(X)      *
!* --------------------------------------------------------------------- *
!*  0.00000E+00   0.50000E+01   0.40000E+01   0.60000E+01   0.12000E+02  *
!*  0.20000E+00   0.59376E+01   0.54720E+01   0.88800E+01   0.16800E+02  *
!*  0.40000E+00   0.72336E+01   0.76160E+01   0.12720E+02   0.21600E+02  *
!*  0.60000E+00   0.90416E+01   0.10624E+02   0.17520E+02   0.26400E+02  *
!*  0.80000E+00   0.11554E+02   0.14688E+02   0.23280E+02   0.31200E+02  *
!*  0.10000E+01   0.15000E+02   0.20000E+02   0.30000E+02   0.36000E+02  *
!*  0.12000E+01   0.19650E+02   0.26752E+02   0.37680E+02   0.40800E+02  *
!*  0.14000E+01   0.25810E+02   0.35136E+02   0.46320E+02   0.45600E+02  *
!*  0.16000E+01   0.33826E+02   0.45344E+02   0.55920E+02   0.50400E+02  *
!*  0.18000E+01   0.44082E+02   0.57568E+02   0.66480E+02   0.55200E+02  *
!*  0.20000E+01   0.57000E+02   0.72000E+02   0.78000E+02   0.60000E+02  *
!* --------------------------------------------------------------------- *
!*                                                                       *
!* Press any key to continue                                             *
!*                                                                       *
!* --------------------------------------------------------------------- *
!* Reference: From Numath Library By Tuan Dang Trong in Fortran 77       *
!*            [BIBLI 18].                                                *
!*                                                                       *
!*                                F90 Release 1.0 By J-P Moreau, Paris   *
!*                                         (www.jpmoreau.fr)             *
!*************************************************************************  
PROGRAM TEST_HORNER

REAL*8 A(25), RES(25)
REAL*8 DX,X,X1,X2
LOGICAL FLG
  
  N = 4       !ORDER OF POLYNOMIAL
  NP = 11     !NUMBER OF POINTS

  !define coefficients
  A(1) = 1.d0
  A(2) = 2.d0
  A(3) = 3.d0
  A(4) = 4.d0
  A(5) = 5.d0

  X1 = 0.d0   !X BEGIN
  X2 = 2.d0   !X END
  K=3         !NUMBER OF DERIVATIVES 
  FLG=.TRUE.  !EVALUATE ONLY

  DX = (X2-X1) / (NP-1)

  X = X1 - DX

  print *,' '
  print *,'     EVALUATE A POLYNOMIAL AND ITS DERIVATIVES BY HORNER''S METHOD'
  print *,' ' 
  print *,'             Example:  P(X) = X^4 + 2X^3 +3X^2 + 4X + 5'
  print *,' '
  print *,'       X          P(X)          P''(X)         P"(X)         P"''(X)'
  print *,' ---------------------------------------------------------------------'

  DO I=1, NP
    X = X + DX
    CALL HORNER(N,A,X,K,FLG,RES)
    write(*,10) X, (RES(J), J=1,K+1)
  END DO
  print *,' ---------------------------------------------------------------------'
  print *,' '  	  

10 format(5E14.5)

END


      SUBROUTINE HORNER (N,A,X0,K,B,R)
!---------------------------------------------------------------------
!     EVALUATE A POLYNOMIAL AND ITS DERIVATIVES BY HORNER'S METHOD
!
!     INPUTS:
!     N       ORDER OF POLYNOMIAL
!     A       VECTOR OF SIZE N+1 STORING THE COEFFICIENTS OF
!             POLYNOMIAL IN DECREASING ORDERS OF POWERS
!             I.E. P(X) = A(1)*X**N+A(2)*X**(N-1)+...+A(N)*X+A(N+1)
!     X0      GIVEN ARGUMENT
!     K       MAXIMUM ORDER OF DERIVATIVES ASKED FOR
!     B       FLAG
!             = .TRUE.  EVALUATION ONLY
!             = .FALSE. DETERMINATION OF POLYNOMIAL COEFFICIENTS
!                NEAR X0
!     R       VECTOR OF SIZE K+1 CONTAINS:
!             -  VALUES P(X0), P'(X0), P''(X0),..PK(X0),
!                IF B IS TRUE AND IF K < N
!             -  THE N+1 COEFFICIENTS OF POLYNOMIAL
!                Q(X-X0) = R(1)+R(2)*(X-X0)+...+R(N+1)*(X-X0)**(N),
!                IF B IS FALSE AND IF K = N.
!
!     REFERENCE:
!     ALGORITHM 337, COLLECTED ALGORITHMS FROM CACM, W.PANKIEWICKZ
!---------------------------------------------------------------------
      REAL*8 A(*),R(*),X0
      LOGICAL B
      RR = A(1)
      DO 10 I = 1,K+1
   10 R(I) = RR
      DO 20 J = 2,N+1
      R(1) = R(1)*X0+A(J)
      NMJ = N-J+1
      IF (NMJ.GT.K) THEN
      L = K
      ELSE
      L = NMJ
      END IF
      DO 15 I = 2,L+1
      R(I) = R(I)*X0+R(I-1)
   15 CONTINUE
   20 CONTINUE
      IF (B) THEN
      L = 1
      DO 25 I = 2,K+1
      L = L*(I-1)
   25 R(I) = R(I)*L
      ELSE
      END IF
      RETURN
      END

!end of file thorner.f90