function y = log1p(z)
%LOG1P  Compute log(1+z) accurately.
%   LOG1P(Z) computes log(1+z), compensating for the roundoff in 1+z.
%   For small z, log1p(z) is approximately z, whereas log(1+z) can be zero.
%
%   See also LOG, EXPM1.

%   Algorithm due to W. Kahan, unpublished course notes.
%   Copyright 1984-2003 The MathWorks, Inc. 
%   $Revision: 1.1.6.3 $  $Date: 2004/08/20 19:50:28 $

u = 1+z;
p = (u<=0) | isinf(u) | isnan(u) | (imag(z)~=0);
m = ~(u==1 | p);
y = z;
y(m) = log(u(m)).*(z(m)./(u(m)-1));
y(p) = log(u(p));