!****************************************************
!* Program to demonstrate Horner_Shift() Subroutine *
!* ------------------------------------------------ *
!* Reference: BASIC Scientific Subroutines, Vol. II *
!* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].*
!*                                                  *
!*              F90 Version by J.-P. Moreau, Paris. *
!*                       (www.jpmoreau.fr)          *
!* ------------------------------------------------ *
!* SAMPLE RUN:                                      *
!*                                                  *
!* Input the five coefficients:                     *
!*                                                  *
!*    A( 0) =  1                                    *
!*    A( 1) =  4                                    *
!*    A( 2) =  6                                    *
!*    A( 3) =  4                                    *
!*    A( 4) =  1                                    *
!*                                                  *
!* What is the expansion point: 1                   *
!*                                                  *
!* The shifted coefficients are:                    *
!*                                                  *
!*    B( 0) =  16.                                  *
!*    B( 1) =  32.                                  *
!*    B( 2) =  24.                                  *
!*    B( 3) =   8.                                  *
!*    B( 4) =   1.                                  *
!*                                                  *
!****************************************************
PROGRAM DEMO_HORNER

real*8, allocatable :: C(:,:), A(:), B(:)
real*8  x0
integer i

  allocate(C(0:4,0:5),A(0:5),B(0:5)) 
  print *,' '
  print *,'Input the five coefficients:'
  print *,' '
  do i = 0, 4
    write(*,25,advance='no')  i
	read *, A(i)
  end do
  print *,' '
  write(*,"(' What is the expansion point: ')",advance='no')
  read *, x0

  call Horner_Shift(x0,A,B,C);

  print *,' '
  print *,'The shifted coefficients are:'
  print *,' '
  do i = 0, 4 
    write(*,50)  i, B(i) 
  end do
  print *,' '
  
  stop
25 format('  A(',i2,') = ')
50 format('  B(',i2,') = ',f6.0)
end

!***************************************************
!*       Horner's shifting rule subroutine         *
!* ----------------------------------------------- *
!* This routine takes a given quartic polynomial   *
!* and converts it to a Taylor expansion.          *
!* The input series coefficients are A(i), the     *
!* expansion point is x0. The shifted coefficients *
!* are returned in B(i).                           *
!***************************************************
Subroutine Horner_Shift(x0,A,B,C)  
  integer i, j
  real*8  C(0:4,0:5), A(0:5), B(0:5)  
  real*8 x0
  do j = 0, 4
    C(j,0) = A(4-j)
  end do
  i=0
  do while (i<5) 
    C(0,i+1) = C(0,i)
    j = 1
    do while (j <= 4-i) 
      C(j,i+1) = x0*C(j-1,i+1) + C(j,i)
      j = j + 1
    end do
    i = i + 1
  end do
  do i = 0, 4  
    B(4-i) = C(i,4-i+1)
  end do 
  return
end

!End of file horner.f90