function [p, e] = project(A, b)
% project Project a vector b onto the column space of A.
%
% p = project(A, b) returns the orthogonal projection of a
% vector b onto the column space of A.
%
% [p, e] = project(A, b) also returns the vectors e = b - p.
% p is the projection of b onto the column space of A.
% e is the projection of b onto the left nullspace of A.
% Notice that b = p + e and p' * e = 0.
%
% See also projmat.
P = projmat(A);
p = P * b;
e = b - p;