function [u, nx] = householder(x)
% Constructs Householder reflector
%
% [u, nx] = householder(x)
%
% Return a vector u such that H*x = nx*e1, where
% H = I - 2*u*u' and e1 = eye(length(x), 1);
% The vector returned is normalized so that norm(u) = 1, so we can omit
% the division by norm(u)^2 = u'*u later when we use it.

u = x;
nx = -norm(x)*sign(x(1));
u(1) = x(1) - nx;
u = u / norm(u);