% ========================================================================
% quadratic.m -- Coefficients are read in from keyboard
%                Roots of Quadratic are printed to screen.
%                                                                      
% Note : Naive implementation of quadratic equation solver. This algorithm
%        does not take into account possible loss of accuracy when two
%        floating point numbers of almost equal size are subtracted.
%
% Written By : Mark Austin                                       July 1997
% ========================================================================

% Print Welcome Message 

  disp('Welcome to the Quadratic Equation Solver (Version 1)');
  disp('====================================================');

% Prompt User for Coefficients of Quadratic Equation 

  disp('Please enter coefficients for equation a.x^2 + b.x + c');
  A = input ('Enter coefficient a : ');
  B = input ('Enter coefficient b : ');
  C = input ('Enter coefficient c : ');

% Print Quadratic Equation to Screen

  fprintf('Equation you have entered is : %g.x^2 + %g.x + %g\n', ...
          A, B, C);

% Compute Roots of simplified equations : A equals zero 

  RootsFound = 0;
  if A == 0 & B == 0,
     fprintf('Cannot solve extremely degenerate equation' );
     fprintf('%14.8g = 0.0\n', C );
     RootsFound = 1;
  end;

  if A == 0 & B ~= 0 & RootsFound == 0,
     Root1 = - C/B;
     fprintf('Degenerate root : Root = %14.8g\n', Root1 );
     RootsFound = 1;
  end;

% Compute Roots of Quadratic Equation : A not equal to zero 

  if RootsFound == 0,
     Discriminant = discriminant(A,B,C);   % Compute discriminant of
                                           % quadratic equation.
     if Discriminant >= 0,                 % Case for two real roots 

        Root1 = -B/2.0/A - sqrt( Discriminant )/2.0/A;
        Root2 = -B/2.0/A + sqrt( Discriminant )/2.0/A;

        fprintf('Two real roots : Root1 = %14.8g\n', Root1 );
        fprintf('               : Root2 = %14.8g\n', Root2 );
     else                             % Case for complex roots 

        fprintf('Two complex roots : Root1 = %14.8g + %14.8g i\n', ...
                -B/2.0/A,  sqrt( -Discriminant )/2.0/A);
        fprintf('                  : Root2 = %14.8g + %14.8g i\n', ....
                -B/2.0/A, -sqrt( -Discriminant )/2.0/A);
     end;
  end;