!****************************************************
!*   Program to demonstrate the series inversion    *
!* ------------------------------------------------ *
!* 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:                                      *
!*                                                  *
!* What is the degree of the input polynomial: 2    *
!*                                                  *
!* What is the degree of the inverted polynomial: 6 *
!*                                                  *
!* Input the polynomial coefficients:               *
!*                                                  *
!*    A( 0) = ? 1                                   *
!*    A( 1) = ? .1                                  *
!*    A( 2) = ? .01                                 *
!*                                                  *
!* The inverted polynomial coefficients are:        *
!*                                                  *
!*    B( 0) =  1.000000                             *
!*    B( 1) = -0.100000                             *
!*    B( 2) =  0.000000                             *
!*    B( 3) =  0.001000                             *
!*    B( 4) = -0.000100                             *
!*    B( 5) =  0.000000                             *
!*    B( 6) =  0.000001                             *
!*                                                  *
!****************************************************
PROGRAM Recipro

real*8, allocatable :: A(:), B(:)
integer i, m, n

  print *,' '
  write(*,"(' What is the degree of the input polynomial: ')",advance='no')
  read *, n
  write(*,"(' What is the degree of the inverted polynomial: ')",advance='no')
  read *, m

  allocate(A(10),B(10))

  print *,'Input the polynomial coefficients:'
  print *,' '

  do i = 0, n
    write(*,25,advance='no') i 
	read *, A(i)
  end do

  call Reciprocal(m,n,A,B)

  print *,' '
  print *,'The inverted polynomial coefficients are:'
  print *,' '
  do i = 0, m
    write(*,50)  i, B(i)
  end do 
  print *,' '
  stop
25 format('   A(',i2,') = ') 
50 format('   B(',i2,') = ',e13.6)  
end


!******************************************************
!*       Reciprocal power series subroutine           *
!* Reference: Computational analysis by Henrici       *
!* -------------------------------------------------- *
!* The input series coefficients are A(i). The output *
!* series coefficients are B(i). The degree of the    *
!* input polynomial is n. The degree of the inverted  *
!* polynomial is m. The routine takes care of norma-  *
!* lization using l = A(0).                           *
!******************************************************
Subroutine Reciprocal(m,n,A,B)  
! Label used: 100
  real*8  A(10), B(10)
  real*8 l
  integer i,j,m,n
  l = A(0)
  do i = 0, n
    A(i) = A(i) / l; B(i) = 0.d0
  end do
  ! Clear arrays
  do i = n+1, m
    A(i) = 0.d0; B(i) = 0.d0
  end do
  ! Calculate the B(i) coefficients
  B(0) = 1.d0
  do i = 1, m
    j = 1
100 B(i) = B(i) - A(j) * B(i - j)
    j = j + 1
    if (j <= i) goto 100
  end do
  ! Un-normalize the A(i) and B(i)
  do i = 0, m
    A(i) = A(i) * l
    B(i) = B(i) / l
  end do
  return
end


! End of file recipro.f90