% given an n-vector x, this function computes an n-vector v with v(1)
% = 1 such that (I - 2*v*v'/(v'*v))*x is zero in all but the first
% component (from Golub and Van Loan, page 196)

function v = house(x)

n = length(x);
mu = norm(x);
v = x;
if mu ~= 0,
  beta = x(1) + sign(x(1))*mu;
  v(2:n) = v(2:n)/beta;
end
v(1) = 1;