function [w,gamma]=house(x)
%
% Computes Householder transformation 
%  H = I - 2ww' where  Hx= gamma e_1, gamma=norm(x)
%     
%  input
%   x- vector to be transformed
%
%  output
%   gamma  - norm(x)
%   w      - Householder vector
%
    n=length(x);
    if n == 1
%
%  Special case of length 1
%
     gamma=abs(x);
     if x < 0
       w=1;
     else
       w=0;
     end;
    else
      alpha=norm(x(2:n));
      gamma=norm([x(1); alpha]);
      sqrt2=sqrt(2);
    if gamma >0  % In case x=0 set w=0
      if x(1) > 0
%
%  Parlett's version, sign unchanged
%
        ratio=x(1)/gamma;
    if alpha > 0
        beta=sqrt2*alpha/sqrt(1+ratio);
        w1=-alpha/(sqrt2*sqrt(gamma)*sqrt(gamma+x(1)));
    else
        w1=0; beta=1;
    end;
      else
        beta=sqrt2*sqrt(gamma)*sqrt(gamma-x(1)); w1=(x(1)-gamma)/beta;
      end;
      w=[w1; x(2:n)/beta];
    else
      w=zeros(n,1);
    end;      
    end;