%  DEMOFL  demonstrates some phenomena related
%  to floating point arithmetic.
%
%  Matthias Heinkenschloss
%  Department of Computational and Applied Mathematics
%  Rice University
%  Feb 12, 2001
%
clear; clc;
disp('                                                                  ')
disp('  Let the real number  x  be written as                           ')
disp('             0       -1      -2      -    -t       -(t+1)      e  ')
disp('  x = ( d * b + d * b + d * b +..... d * b +  d * b + ... ) * b , ')
disp('         0       1       2            t-1      t                  ')
disp('                                                                  ')
disp('  where   b > 1  is the base and  e  is the exponent. If x is not ')
disp('  equal to zero, then d_0 is normalized to be nonzero             ')
disp('  (The base  b  and the exponent  e  are integer.)                ')
disp('  If the exponent is in [ e_low, e_upp ], then the floating point ')
disp('  representation of x is given by                                 ')
disp('                 0       -1      -2           -t     e            ')
disp('  fl(x) = ( d * b + d * b + d * b +..... d * b  ) * b  ,          ')
disp('             0       1       2            t                       ')
disp('  The t-th digit of  fl(x)  is given by                           ')
disp('         -                                                        ')
disp('         d    = d       if   d    <  0.5*b ,                      ') 
disp('          t-1    t-1          t                                   ')
disp('         -                                                        ')
disp('         d    = d   +1  if   d    >= 0.5 * b.                     ')
disp('          t-1    t-1          t                                   ')
disp('                                                                  ')


disp(' Hit a key to continue '); pause

clear; clc;
disp('                              ')
disp('                              ')
disp('                              ')
disp(' Example 1:  Floating Point Representation of Real Numbers')
disp('                              ')
disp(' Base b = 10                  ')
disp(' Mantissa length t = 4');
disp(' lower boundary of exponent range e_low = 10^(-37)');
disp(' upper boundary of exponent range e_upp = 10^(37)');
b = 10;
t = 4;

x = 1;
while( ~isempty(x) )
  disp('    ')
  x = input(' Enter a number (Enter RETURN to quit)  x = ');
  if( isempty(x) ) break; end
  flx = sfl(x,t,'r');
  if( x == 0 )
      abserr = abs(x - flx);
      fprintf(1, '  x                  = %12.6e  \n', x );
      fprintf(1, '  fl(x)              = %12.6e  \n', flx);
      fprintf(1, '  |x - fl(x) |       = %12.6e  \n', abserr);
  else
      abserr = abs(x - flx);
      relerr = abserr / abs(x);
      fprintf(1, '  x                  = %12.6e  \n', x );
      fprintf(1, '  fl(x)              = %12.6e  \n', flx);
      fprintf(1, '  |x - fl(x) |       = %12.6e  \n', abserr);
      fprintf(1, '  |x - fl(x)| / |x|  = %12.6e  \n', relerr);
  end
end



clear;clc;
disp('                              ')
disp('                              ')
disp('                              ')
disp(' Example 2:  Floating Point Operations with Real Numbers')
disp('                              ')
disp(' Base b = 10                  ')
disp(' Mantissa length t = 4');
disp(' lower boundary of exponent range e_low = 10^(-37)');
disp(' upper boundary of exponent range e_upp = 10^(37)');
b = 10;
t = 4;

x = 1;
while( ~isempty(x) )
  disp('    ') 
  x = input(' Enter a number (hit RETURN to quit)  x = ');
  if( isempty(x) ) break; end
  y = input(' Enter a number (hit RETURN to quit)  y = ');
  if( isempty(x) | isempty(y) ) break; end;
  xpy    =  x + y;
  fxpy   = sflop( x, y, t, 'r', '+' );
  abserr = abs(xpy - fxpy);
  relerr = abserr / abs(xpy);
  fprintf(1, '  x + y              = %12.6e  \n', xpy );
  fprintf(1, '  fl(fl(x)+fl(y))    = %12.6e  \n', fxpy);
  fprintf(1, '  absolute error     = %12.6e  \n', abserr);
  fprintf(1, '  relative error     = %12.6e  \n\n', relerr);
  xmy    =  x - y;
  fxmy   = sflop( x, y, t, 'r', '-' );
  abserr = abs(xmy - fxmy);
  relerr = abserr / abs(xmy);
  fprintf(1, '  x - y              = %12.6e  \n', xmy );
  fprintf(1, '  fl(fl(x)-fl(y))    = %12.6e  \n', fxmy);
  fprintf(1, '  absolute error     = %12.6e  \n', abserr);
  fprintf(1, '  relative error     = %12.6e  \n\n', relerr);
  xty    =  x * y;
  fxty   = sflop( x, y, t, 'r', '*' );
  abserr = abs(xty - fxty);
  relerr = abserr / abs(xty);
  fprintf(1, '  x * y              = %12.6e  \n', xty );
  fprintf(1, '  fl(fl(x)*fl(y))    = %12.6e  \n', fxty);
  fprintf(1, '  absolute error     = %12.6e  \n', abserr);
  fprintf(1, '  relative error     = %12.6e  \n\n', relerr);
  xdy    =  x / y;
  fxdy   = sflop( x, y, t, 'r', '/' );
  abserr = abs(xdy - fxdy);
  relerr = abserr / abs(xdy);
  fprintf(1, '  x / y              = %12.6e  \n', xdy );
  fprintf(1, '  fl(fl(x)/fl(y))    = %12.6e  \n', fxdy);
  fprintf(1, '  absolute error     = %12.6e  \n', abserr);
  fprintf(1, '  relative error     = %12.6e  \n', relerr);
end