function  sz = rtsumsqs(x)
%  sz = rtsumsqs(x) = norm(x)  demonstrates scaling and
%  the correct treatment of infinity and/or NaN in  x :
%    if any(isinf(x(:))),  sz = inf  regardless of  NaN
%    elseif any(isnan(x(:))),  sz = NaN  end
%  Some Matlab versions'  norm  get this wrong.

%                                 W. Kahan,  7 Apr. 2012

%  Constants:  All turn out to be integer powers of  2 .
    rttau = sqrt(realmin) ;  rtomega = sqrt(realmax) ;
    Huge = 1/(eps*rttau) ;   huge = rtomega ;
    Tiny = rttau/eps ;       tiny = sqrt(eps)/rtomega ;

%  Turn input  x  into a real column:
x = x(:) ;  %... turns  x  into a column
if any(imag(x)),  x = [real(x); imag(x)] ;  end

n = length(x) ;  %...  so  x  is  n-by-1 real
if (n < 2),  sz = abs(x) ;  return,  end
big = huge/sqrt(n) ;
ahat = 0 ;  sz = 0 ;  rsz = 0 ;

%  First pass over the data:
for  j = 1:n
    abx = abs(x(j)) ;
    if isinf(abx),  sz = abx ;  return
      elseif isnan(abx),  sz = abx ;
      else  ahat = max(ahat, abx) ;  end
    if (ahat < big)
%...    perform Compensated Summation:
        rsz = abx*abx + rsz ;
        t = sz ;  sz = t + rsz ;
        rsz = (t - sz) + rsz ;  %... compensation
      end,  end  %... of first pass

%  Filter out  NaN :
if isnan(sz),  return,  end

%  Test whether scaling by  sc ~= 1  is necessary:
if (ahat >= big),  sc = tiny ;  %... rteps/rtomega
  elseif (ahat <= Tiny),  sc = Huge ;  %... 1/(eps*rttau)
  else  sc = 1 ;  end  %... scaling will be unnecessary

if (sc ~= 1)  %...  a second pass over the data needed:
    sz = 0 ;  rsz = 0 ;
    for  j = 1:n
        abx = x(j)*sc ;
%...    perform Compensated Summation of scaled values:
        rsz = abx*abx + rsz ;
        t = sz ;  sz = t + rsz ;
        rsz = (t - sz) + rsz ;  %... compensation
      end,  end  %... of necessary second pass

sz = sqrt(sz)/sc ;  %... overunderflows only if it must
%...  End of  rtsumsqs