!**************************************************** !* 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