% Some examples for playing with the fl( x ) and QuadRoot( a, b, c ) functions
%
% Step one: initialize the floating point system.  See the comments in fl.m by typing

help fl

% Notice that using the notation in the handout of floating point arithmetic we need
% to indicate t, U, L, and beta.  beta = 10 is assumed.  Also, the method (rounded or
% truncated) must be indicated.  The other values are set
% as global variables which are then used by the function fl( x )
%
% Step 2: declare fp_t fp_U fp_L fp_method

fp_t = 3;    % three decimal digits of accuracy
fp_U = 4;    % largest exponent allowed
fp_L = -3;   % smallest exponent allowed
fp_method = 'R' % use rounded arithmetic

%
fl( 3.141568 )

fl( 99999.99 )

fl( 999.999 )

fl( 0.0001467 )

fl( 0.00001467 )

% Let's compute the largest root of x^2 + 4.0145 x + 0.01267
% xf will equal the result when the above described floating point number system
% is used.  x will equal the result when matlab precision is used

a = 1
b = 4.0145
c = 0.01267
[ xf x ] = QuadRoot( a, b, c )

% Notice that there is quite a difference in the result.  The RELATIVE error (see text)
% is given by abs( x - xf ) / abs( x )

abs( x - xf ) / abs( x )

% Now, let's see what happens for different values of c.  We will
% let c equal a vector of equally spaced points between 0 and 1 (inclusive)

c = linspace( 0, 1, 20 );

% x(i) will now equal the largest root of x^2 + b x + c(i)

for i=1:20
   [ xf(i) x(i) ] = QuadRoot( a, b, c(i) );
end

% Let's look at a plot of xf and x as a function of c

plot( c, x, '-', c, xf, '-.' )

% Notice that x(1) = 0 (since it is the largest root of x^2 + bx + 0)
% Just so we don't divide by 0, let's set it to 0.00001

input('continue');

x(1) = 0.00001

% Now let's plot the relative error abs( x(i) - xf(i) ) / abs( x(i) )
% Note that the "./" indicates element-wise division of vectors

plot( c, abs( x - xf ) ./ abs( x ) )
xlabel( 'c' )
ylabel( '| x - xf |/| x |' )
title( 'relative error in the largest root of x^2 + 40145 x + c' )

% Notice that the relative error is very large when c is small!

% Exercise: fix the QuadRoot function so it computes it accurately both when 
% b is positive and when it is negative

   
   