program ALG033;
{     HERMITE INTERPOLATION

      TO OBTAIN THE COEFFICIENTS OF THE HERMITE INTEPOLATING
      POLYNOMIAL H ON THE (N+1) DISTINCT NUMBERS X(0), ..., X(N)
      FOR THE FUNCTION F:

      INPUT:   NUMBERS X(0), X(1), ..., X(N); VALUES F(X(0)), F(X(1)),
               ..., F(X(N)) AND F'(X(0)), F'(X(1)), ..., F'(X(N)).

      OUTPUT:  NUMBERS Q(0,0), Q(1,1), ..., Q(2N + 1,2N + 1) WHERE

               H(X) = Q(0,0) + Q(1,1) * ( X - X(0) ) + Q(2,2) *
                      ( X - X(0) )**2 + Q(3,3) * ( X - X(0) )**2 *
                      ( X - X(1) ) + Q(4,4) * ( X - X(0) )**2 *
                      ( X - X(1) )**2 + ... + Q(2N + 1,2N + 1) *
                      ( X - X(0) )**2 * ( X - X(1) )**2 * ... *
                      ( X - X(N - 1) )**2 * (X - X(N) ).
                                                                               }
type
   FILENAME = string [ 14 ];
var
   Q : array [ 0..25 , 0..25 ] of real;
   X : array [ 0..12 ] of real;
   Z : array [ 0..25 ] of real;
   XX, S : real;
   I,J,K,N,FLAG : integer;
   INP,OUP : text;
   NAME : FILENAME;
   OK : boolean;
   A : char;
function F ( X : real ) : real;
   begin
      F := 1.0
   end;
function FP ( X : real ) : real;
   begin
      FP := 0.0
   end;
procedure OUTPUT;
      begin
         writeln ('Choice of output method: ');
         writeln ('1. Output to screen ');
         writeln ('2. Output to text file ');
         writeln ('Please enter 1 or 2 ');
         readln ( FLAG );
         if ( FLAG = 2 ) then
            begin
               writeln ('Input the file name in the form - drive:name.ext ');
               readln ( NAME );
               assign ( OUP , NAME )
            end
         else  assign ( OUP , 'CON' );
         rewrite ( OUP );
         writeln(OUP,'HERMITE INTERPOLATING POLYNOMIAL');
         writeln(OUP)
      end;
procedure INPUT;
   begin
      writeln('This is Hermite interpolation.');
      OK := false;
      while ( not OK ) do
         begin
            writeln ('Choice of input method: ');
            writeln ('1. Input entry by entry from keyboard ');
            writeln ('2. Input data from a text file ');
            writeln ('3. Generate data using a function F ');
            writeln ('Choose 1, 2, or 3 please ');
            readln ( FLAG );
            if ( FLAG = 1 ) or ( FLAG = 2 ) or ( FLAG = 3 ) then OK := true
         end;
      case FLAG of
         1 : begin
                OK := false;
                while ( not OK ) do
                   begin
                      writeln ('Input the number of data points minus 1 ');
                      readln ( N );
                      if ( N > 0 ) then
                         begin
                            OK := true;
                            for I := 0 to N do
                               begin
                                  write ('Input X(',I,'), F(X(',I,')), and ');
                                  write ('F''(X(',I,')) separated by spaces ');
                                  writeln ('- include decimal point ');
                                  readln ( X[I], Q[2 * I,0], Q[2 * I + 1,1] )
                               end
                         end
                      else writeln ('Number must be a positve integer ')
                   end
             end;
         2 : begin
                write ('Has a text file been created with the data in three ');
                writeln ('columns? ');
                writeln ('Enter Y or N ');
                readln ( A );
                if ( A = 'Y' ) or ( A = 'y' ) then
                   begin
                      write ('Input the file name in the form - ');
                      writeln ('drive:name.ext ');
                      writeln ('as example:   A:DATA.DTA ');
                      readln ( NAME );
                      assign ( INP, NAME );
                      reset ( INP );
                      OK := false;
                      while ( not OK ) do
                         begin
                            writeln ('Input the number of data points minus 1 ');
                            readln ( N );
                            if ( N > 0 ) then
                               begin
                                  for I := 0 to N do
                                     readln ( INP, X[I], Q[2 * I,0],
                                              Q[2 * I + 1,1] );
                                  close ( INP );
                                  OK := true
                               end
                            else writeln ('Number must be a positive integer ')
                         end
                   end
                else
                   begin
                      write ('Please create the input file in three column ');
                      writeln ('form with the X values, F(X), and ');
                      writeln ('F''(X) values in the corresponding columns. ');
                      write ('The program will end so the input file can ');
                      writeln ('be created. ');
                      OK := false;
                   end
             end;
         3 : begin
                write ('Have the functions F and FP been created in ');
                writeln ('the program immediately preceding ');
                writeln ('the INPUT procedure? ');
                writeln ('Enter Y or N ');
                readln ( A );
                if ( A = 'Y' ) or ( A = 'y' ) then
                   begin
                      OK := false;
                      while ( not OK ) do
                         begin
                            writeln ('Input the number of data points minus 1 ');
                            readln ( N );
                            if ( N > 0 ) then
                               begin
                                  for I := 0 to N do
                                     begin
                                        write ('Input X(',I,') - include ');
                                        writeln ('decimal point ');
                                        readln ( X[I] );
                                        Q[2 * I,0] := F( X[I] );
                                        Q[2 * I + 1,1] := FP( X[I] )
                                     end;
                                  OK := true
                               end
                            else writeln ('Number must be a positive integer ')
                         end
                   end
                else
                   begin
                      write ('The program will end so that the functions F ');
                      writeln ('and FP can be created. ');
                      OK := false
                   end
             end
      end
   end;
procedure OUTPUT2;
   begin
      writeln (OUP,'The input data follows: ');
      writeln (OUP,'  X, F(X), F''(X) ');
      for I := 0 to N do
         writeln (OUP, X[I], Q[2 * I,0], Q[2 * I + 1,1] );
      writeln(OUP);
      writeln (OUP,'The Coefficients of the Hermite Interpolation Polynomial ');
      writeln (OUP,'in order of increasing exponent follow: '); writeln(OUP);
      for I := 0 to K do writeln (OUP, Q[I,I] )
   end;
   begin
      INPUT;
      if ( OK ) then
         begin
{           STEP 1                                                             }
            for I := 0 to N do
               begin
{                 STEP 2                                                       }
                  Z[2 * I] := X[I];
                  Z[2 * I + 1] := X[I];
                  Q[2 * I + 1,0] := Q[2 * I,0];
{                 STEP 3                                                       }
                  if ( I <> 0 ) then
                     Q[2 * I,1] := ( Q[2 * I,0] - Q[2 * I - 1, 0] )
                                    / ( Z[2 * I] - Z[2 * I - 1] );
               end;
{           STEP 4                                                             }
            K := 2 * N + 1;
            for I := 2 to K do
               for J := 2 to I do
                 Q[I,J] := ( Q[I,J - 1] - Q[I - 1,J - 1] )
                            / ( Z[I] - Z[I - J] );
{           STEP 5                                                             }
            OUTPUT;
            OUTPUT2;
            writeln('Do you wish to evaluate this polynomial?');
            writeln('Enter Y or N');
            readln(A);
            if (A = 'Y') or (A = 'y') then
               begin
                  writeln('Enter a point at which to evaluate');
                  readln(XX);
                  S := Q[K,K]*(XX-Z[K-1]);
                  for I := 2 to K do
                     begin
                        J := K - I + 1;
                        S := (S + Q[J,J])*(XX-Z[J-1])
                     end;
                  S := S + Q[0,0];
                  writeln(OUP,'x-value and interpolated-value');
                  writeln(OUP,XX,S)
               end
         end;
         close(OUP)
   end.