%
%  PURPOSE:
%  SFL  computes the floating point represenation of x.
%  This is just a simulation of what the floating point 
%  represenation of x would be. (sfl = Simulate FLoating 
%  point)
% 
%  It assumed that the base is 10 and that there
%  is no restriction on the exponent range.
%  The mantissa length is m, i.e. if
%
%                                 -               e
%          x = d  . d  d  d ..... d    d   ... * 10 ,
%               0    1  2  3       m-1  m
%
%  then
%                                         e
%     sfl(x) = d  . d  d  d ..... d   * 10 .
%               0    1  2  3       m-1   
%
%  The m-th digit of  sfl(x)  may be different from the
%  m-th digit of. If FLOATING POINT ARITHMETIC is used, 
%  then
%                -
%         d    = d       if   d  <= 4 ,  
%          m-1    m            m
%                -
%         d    = d +1    if   d  => 5 .
%          m-1    m            m
%
%  If CHOPPED ARITHMETIC is used, then
%              -
%         d  = d  .  
%          m    m         
%   
% 
%  CALLING SEQUENCE:
%
%  function [flx] = sfl( x, m, arith )
% 
%  INPUT PARAMETERS:
% 
%  x      real
%         floating point operation
% 
%  m      integer
%         mantissa length
% 
%  arith  string
%         arith = 'c'   chopped arithmetic
%         arith = 'r'   rounded arithmetic
% 
%
%
%
%
%  Matthias Heinkenschloss
%  Department of Computational and Applied Mathematics
%  Rice University
%  Feb 12, 2001
%
% 


function [flx] = sfl( x, m, arith )

%  The following internal variables are used:
%
%  iexp   is the exponent e,  
%  xm     is the mantissa (with leading 0.)
%  xrd    is the rounded value of x  (rounded to m digits)
%  xchp   is the chopped value of x  (chopped to m digits)

%  if  x = 0. we don't have to do anything

if( x == 0. ) 
    flx = 0.;
    return
end


%  Determine the exponent
tmp = log10(abs(x));
if ( tmp > 0. ) 
     ixexp = 1 + fix( tmp );
else
     ixexp = fix( tmp );
end
   
 
%  determine the mantissa
xm = x * 10^(-ixexp);

if( xm > 0 ) 
    xrd  = fix( xm * 10^m + 0.5 ) * 10.^(ixexp-m);
    xchp = fix( xm * 10^m ) * 10^(ixexp-m);
end
if( xm < 0 ) 
    xrd  = fix( xm * 10^m - 0.5 ) * 10.^(ixexp-m);
    xchp = fix( xm * 10^m ) * 10^(ixexp-m);
end

if( arith == 'c' )  
    flx = xchp;
end
if( arith == 'r' )  
    flx = xrd;
end