!****************************************************
!*   Program to demonstrate the series reversion    *
!* ------------------------------------------------ *
!* 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: 3    *
!*                                                  *
!* Input the coefficients as prompted:              *
!*                                                  *
!*    A( 0) = 1                                     *
!*    A( 1) = 1                                     *
!*    A( 2) = 1                                     *
!*    A( 3) = 1                                     *
!*                                                  *
!* The reversed polynomial coefficients are:        *
!*                                                  *
!*    B( 0) =  9                                    *
!*    B( 1) =  1                                    *
!*    B( 2) =  0                                    *
!*    B( 3) = -1                                    *
!*    B( 4) =  0                                    *
!*    B( 5) =  3                                    *
!*    B( 6) =  0                                    *
!*    B( 7) = -12                                   *
!*                                                  *
!****************************************************
PROGRAM Reverse

real*8  A(10), B(10)
integer i, n

  print *,' '
  write(*,"(' What is the degree of the input polynomial: ')",advance='no')
  read *, n
  print *, 'Input the coefficients as prompted:'
  print *,' '
  do i = 0, n
    write(*,25,advance='no') i
    read *, A(i)
  end do
  print *,' '

  call Reverse_coeff(A,B)

  if (A(1).ne.0) then
    print *, 'The reversed polynomial coefficients are:'
    print *,' '
    do i = 0, 7
      if (dabs(B(i)) < 1.d-20)  B(i) = 0.d0
      write(*,50) i, B(i) 
    end do
  end if
  print *,' '
25 format('  A(',i2,') = ')
50 format('  B(',i2,') = ',e14.7)
  stop
end

!******************************************************
!* Series reversion - This routine takes a polynomial *
!* Y = A(0) + A(1) * X + ... and returns a polynomial *
!* X = B(0) + B(1) * Y + ... A(1) must be <> 0.       *
!* The degree of reversion is limited to seven.       *
!* -------------------------------------------------- *
!* Reference: CRC Standard Mathematical Tables,       *
!* 24th edition.                                      *
!******************************************************
Subroutine Reverse_coeff(A,B)
  real*8 A(10), B(10)    
  real*8 a1,a2,a3,a4,a5,a6,a7,aa,bb
  a1 = A(1); a2=0.d0
  if (a1.eq.0) then
    print *,'Divide zero error (A(1) must be nonzero).'
    return
  else 
    B(1) = 1.d0 / a1
    aa = 1.d0 / a1
  end if
  bb = aa * aa; aa = aa * bb
  B(2) = -a2 / aa
  a3 = A(3); aa = aa * bb
  B(3) = aa * (2.d0 * a2 * a2 - a1 * a3)
  a4 = A(4); aa = aa * bb
  B(4) = aa * (5.d0 * a1 * a2 * a3 - a1 * a1 * a4 - 5.d0 * a2 * a2 * a2)
  a5 = A(5); aa = aa * bb
  B(5) = 6.d0 * a1 * a1 * a2 * a4 + 3.d0 * a1 * a1 * a3 * a3 + 14.d0 * a2 * a2 * a2 * a2
  B(5) = B(5) - a1 * a1 * a1 * a5 - 21.d0 * a1 * a2 * a2 * a3
  B(5) = aa * B(5)
  a6 = A(6); aa = aa * bb
  B(6) = 7.d0 * a1 * a1 * a1 * a2 * a5 + 7.d0 * a1 * a1 * a1 * a3 * a4 + 84.d0 * a1 * a2 * a2 * a2 * a3
  B(6) = B(6) - a1 * a1 * a1 * a1 * a6 - 28.d0 * a1 * a1 * a2 * a2 * a4
  B(6) = B(6) - 28.d0 * a1 * a1 * a2 * a3 * a3 - 42.d0 * a2 * a2 * a2 * a2 * a2
  B(6) = aa * B(6)
  a7 = A(7); aa = aa * bb
  B(7) = 8.d0 * a1 * a1 * a1 * a1 * a2 * a6 + 8.d0 * a1 * a1 * a1 * a1 * a3 * a5
  B(7) = B(7) + 4.d0 * a1 * a1 * a1 * a1 * a4 * a4 + 120.d0 * a1 * a1 * a2 * a2 * a2 * a4
  B(7) = B(7) + 180.d0 * a1 * a1 * a2 * a2 * a3 * a3 + 132.d0 * a2 * a2 * a2 * a2 * a2 * a2
  B(7) = B(7) - a1 * a1 * a1 * a1 * a1 * a7 - 36.d0 * a1 * a1 * a1 * a2 * a2 * a5
  B(7) = B(7) - 72.d0 * a1 * a1 * a1 * a2 * a3 * a4 - 12.d0 * a1 * a1 * a1 * a3 * a3 * a3
  B(7) = B(7) - 330.d0 * a1 * a2 * a2 * a2 * a2 * a3
  B(7) = aa * B(7)
  B(0) = 0.d0
  aa = A(0)
  do i = 1, 7
    B(0) = B(0) - B(i) * aa
    aa = aa * A(0)
  end do
  return
end


! End of file reverse.f90