program alg028;
{  MULLER'S ALGORITHM 2.8

   To find a solution to f(x) = 0 given three approximations x0, x1
   and x2:

   INPUT:  x0,x1,x2; tolerance TOL; maximum number of iterations NO.

   OUTPUT: approximate solution p or message of failure.

   This implementation allows for a switch to complex arithmetic.
   The coefficients are stored in the vector A, so the dimension
   of A may have to be changed.                                        }
const
   ZERO = 1.0e-20;
var
   A : array[1..50] of real;
   ZR,ZI,GR,GI,F,X : array[1..4] of real;
   CH1R,CH1I,H : array [1..3] of real;
   CDEL1R,CDEL1I,DEL1 : array[1..2] of real;
   CDELR,CDELI,CBR,CBI,CDR,CDI,CER,CEI : real;
   DEL,B,D,E,TOL : real;
   QR,QI,ER,EI,FR,FI,P : real;
   FLAG,N,M,I,J,K,ISW,KK : integer;
   OK : boolean;
   NAME : string [ 30 ];
   INP, OUP : text;
   AA : char;
procedure CADD ( A,B,C,D : real; var E,F : real );
{ procedure to perform complex addition :
   (A + Bi) + (C + Di) -> E + Fi  }
   begin
      E := A + C;
      F := B + D
   end;
procedure CMULT ( A,B,C,D : real; var E,F : real );
{ procedure to perform complex multiplication :
   (A + Bi) * (C + Di) -> E + Fi  }
   begin
      E := ( A * C ) - ( B * D );
      F := A * D + B * C
   end;
procedure CDIV ( A,B,C,D : real; var E,F : real );
{ procedure to perform complex division :
   (A + Bi) / (C + Di) -> E + Fi  }
   var
      G : real;
   begin
      G := C * C + D * D;
      E := ( A * C + B * D ) / G;
      F := ( B * C - A * D ) / G
   end;
function CABS ( A,B : real ) : real;
{ function to compute complex absolute value :
   | A + Bi | = sqrt(A*A + B*B)  }
   var
      C : real;
   begin
      C := sqrt(A * A + B * B);
      CABS := C
   end;
 procedure CSQRT ( A,B : real; var C,D : real );
 {procedure to compute complex square root:
    sqrt( A + Bi) -> C + Di  }
   const
      ZERO = 1.0E-20;
   var
      G,R,T,HP : real;
   begin
      HP := 0.5*pi;
      if ( abs( A ) <= ZERO ) then
         begin
            if ( abs( B ) <= ZERO ) then
               begin
                  R := 0.0;
                  T := 0.0
               end
            else
               begin
                  T := HP;
                  if ( B < 0.0 ) then T := -T;
                  R := abs( B )
               end
         end
      else
         begin
            R :=  sqrt(A * A + B * B) ;
            if (abs(B) < ZERO) then
               begin
                  T := 0.0;
                  if (A < 0.0) then T := pi
               end
            else
               begin
                  T := arctan( B / A );
                  if (A < 0.0) then T := T + pi
               end
         end;
      G := sqrt( R );
      C := G * cos( 0.5 * T );
      D := G * sin( 0.5 * T )
   end;
procedure CSUB ( A,B,C,D : real; var E,F : real );
{ procedure to perform complex subtraction :
   (A + Bi) - (C + Di) -> E + Fi  }
   begin
      E := A - C;
      F := B - D
   end;
procedure INPUT;
   begin
      writeln('This is Mullers Method');
      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 ('Choose 1 or 2 please ');
            readln ( FLAG );
            if ( FLAG = 1 ) or ( FLAG = 2 ) then OK := true
         end;
      case FLAG of
         1 : begin
                OK := false;
                while ( not OK ) do
                   begin
                      writeln('Input the degree n of the polynomial');
                      readln(N);
                      if ( N > 0 ) then
                         begin
                            OK := true;
                            write('Input the coefficients of the');
                            writeln(' polynomial in ascending order');
                            writeln('of exponent at the prompt');
                            N := N+1;
                            for I := 1 to N do
                               begin
                                  J := I-1;
                                  writeln('Input A( ',J,' )');
                                  readln(A[I])
                               end
                         end
                      else writeln(' n must be a positive integer.');
                   end;
             end;
         2 : begin
                writeln('Is there a text file containing the coefficients');
                writeln('in ascending order of exponent?');
                writeln ('Enter Y or N ');
                readln ( AA );
                if ( AA = 'Y' ) or ( AA = 'y' ) then
                   begin
                      OK := true;
                      write ('Input the file name in the form - ');
                      writeln ('drive:name.ext ');
                      writeln ('for example:   A:DATA.DTA ');
                      readln ( NAME );
                      assign ( INP, NAME );
                      reset ( INP );
                      OK := false;
                      while ( not OK ) do
                         begin
                            writeln ('Input n');
                            readln ( N );
                            if ( N > 0 ) then
                               begin
                                  OK := true;
                                  N := N+1;
                                  for I := 1 to N do
                                     read ( INP,A[I]);
                                  close ( INP )
                               end
                            else writeln ('Number must be a positive integer ')
                         end
                   end
                else
                   begin
                      writeln ('Please create the input file.');
                      writeln ('The program will end so the input file can ');
                      writeln ('be created. ');
                      OK := false
                   end
             end
      end;
      if (A[N] = 0) and OK then
         begin
            writeln('Leading coefficient is 0 - error in input');
            OK := false
         end;
      if (N = 2) and OK then
         begin
            P := -A[1]/A[2];
            writeln('Polynomial is linear:  root is ',P);
            OK := false;
         end;
      if OK then
         begin
            OK := false;
            while ( not OK ) do
               begin
                  writeln ('Input tolerance ');
                  readln ( TOL );
                  if (TOL <= 0.0) then
                     writeln ('Tolerance must be positive ')
                  else OK := true
               end;
            OK := false;
            while ( not OK ) do
               begin
                  write('Input maximum number of iterations ');
                  writeln('- no decimal point ');
                  readln ( M );
                  if ( M <= 0 ) then
                     writeln ('Must be positive integer ')
                  else OK := true
               end;
            writeln('Input the first of three starting values');
            readln(X[1]);
            writeln('Input the second of three starting values');
            readln(X[2]);
            writeln('Input the third starting value');
            readln(X[3]);
         end;
   end;
procedure OUTPUT;
   begin
      writeln ('Select output destination ');
      writeln ('1. Screen ');
      writeln ('2. Text file ');
      writeln ('Enter 1 or 2 ');
      readln ( FLAG );
      if ( FLAG = 2 ) then
         begin
            write ('Input the file name in the form - ');
            writeln ('drive:name.ext ');
            writeln ('for example:   A:OUTPUT.DTA ');
            readln ( NAME );
            assign ( OUP, NAME )
         end
      else assign ( OUP, 'CON');
      rewrite ( OUP );
      writeln(OUP,'MULLERS METHOD');
      writeln(OUP,'The input polynomial:');