{***************************************************
*        Program to demonstrate the series         *
*              reversion procedure                 *
* ------------------------------------------------ *
* Reference: BASIC Scientific Subroutines, Vol. II *
* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].*
*                                                  *
*         Pascal 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;
Uses WinCrt;


VAR
        A, B  : Array[0..10] of DOUBLE;
        i, n : INTEGER;
 

{*****************************************************
*           Series reversion subroutine              *
* -------------------------------------------------- *
* 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.                                      *
*****************************************************}
PROCEDURE Reverse_coeff;
Var
        a1,a2,a3,a4,a5,a6,a7,aa,bb : DOUBLE;
Begin
  a1 := A[1];
  IF a1 = 0.0 THEN
  begin
    Writeln(' Divide zero error (A(1) must be nonzero).');
    exit
  end;
  B[1] := 1.0 / a1;
  aa := 1.0 / a1;
  bb := aa * aa; aa := aa * bb;
  B[2] := -a2 / aa;
  a3 := A[3]; aa := aa * bb;
  B[3] := aa * (2.0 * a2 * a2 - a1 * a3);
  a4 := A[4]; aa := aa * bb;
  B[4] := aa * (5.0 * a1 * a2 * a3 - a1 * a1 * a4 - 5.0 * a2 * a2 * a2);
  a5 := A[5]; aa := aa * bb;
  B[5] := 6.0 * a1 * a1 * a2 * a4 + 3.0 * a1 * a1 * a3 * a3 + 14.0 * a2 * a2 * a2 * a2;
  B[5] := B[5] - a1 * a1 * a1 * a5 - 21.0 * a1 * a2 * a2 * a3;
  B[5] := aa * B[5];
  a6 := A[6]; aa := aa * bb;
  B[6] := 7.0 * a1 * a1 * a1 * a2 * a5 + 7.0 * a1 * a1 * a1 * a3 * a4 + 84 * a1 * a2 * a2 * a2 * a3;
  B[6] := B[6] - a1 * a1 * a1 * a1 * a6 - 28.0 * a1 * a1 * a2 * a2 * a4;
  B[6] := B[6] - 28.0 * a1 * a1 * a2 * a3 * a3 - 42.0 * a2 * a2 * a2 * a2 * a2;
  B[6] := aa * B[6];
  a7 := A[7]; aa := aa * bb;
  B[7] := 8.0 * a1 * a1 * a1 * a1 * a2 * a6 + 8.0 * a1 * a1 * a1 * a1 * a3 * a5;
  B[7] := B[7] + 4.0 * a1 * a1 * a1 * a1 * a4 * a4 + 120.0 * a1 * a1 * a2 * a2 * a2 * a4;
  B[7] := B[7] + 180.0 * a1 * a1 * a2 * a2 * a3 * a3 + 132.0 * a2 * a2 * a2 * a2 * a2 * a2;
  B[7] := B[7] - a1 * a1 * a1 * a1 * a1 * a7 - 36.0 * a1 * a1 * a1 * a2 * a2 * a5;
  B[7] := B[7] - 72.0 * a1 * a1 * a1 * a2 * a3 * a4 - 12.0 * a1 * a1 * a1 * a3 * a3 * a3;
  B[7] := B[7] - 330.0 * a1 * a2 * a2 * a2 * a2 * a3;
  B[7] := aa * B[7];
  B[0] := 0.0;
  aa := A[0];
  FOR i := 1 TO 7 DO
  begin
    B[0] := B[0] - B[i] * aa;
    aa := aa * A[0]
  end
End;


{main program}
BEGIN
  ClrScr;
  Writeln;
  Write(' What is the degree of the input polynomial: '); read(n);
  Writeln;
  Writeln(' Input the coefficients as prompted:');
  Writeln;
  FOR i := 0 TO n DO
  begin
    Write('  A(',i,') = '); read(A[i])
  end;
  Writeln;

  Reverse_coeff;

  IF A[1] <> 0 THEN
  begin
    Writeln(' The reversed polynomial coefficients are:');
    Writeln;
    FOR i := 0 TO 7 DO
    begin
      if abs(B[i]) < 1e-20 then B[i] := 0.0;
      Writeln('  B(',i,') = ',B[i]:3:0);
    end
  end;
  Writeln;
  ReadKey; DoneWinCrt
END.

{End of file reverse.pas}