  function v = house(x)
% x is a column n-vector.
% v is a unit length, column n-vector such that if y = (I - 2vv')x 
% then y(2:n) = 0.

alfa = norm(x);
v = x;
if alfa == 0
    v(1) = 1;
    return
else
    if x(1)>=0
        v(1) = v(1)+alfa;
    else
        v(1) = v(1)-alfa;
    end
    % Since (x +- alfa*e1)'(x +- alfa*e1) = 2alfa^2 + 2|x(1)|alfa
    v = v/sqrt(2*alfa*(alfa + abs(x(1))));
end