function y = expm1(z)
%EXPM1  Compute exp(z)-1 accurately.
%    EXPM1(Z) computes exp(z)-1, compensating for the roundoff in exp(z).
%    For small z, expm1(z) is approximately z, whereas exp(z)-1 can be zero.
%
%    (The matrix exponential demo function EXPM1 is now EXPMDEMO1.)
%
%    See also EXP, LOG1P, EXPMDEMO1.

%   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:27 $

u = exp(z);
p = (u==0) | isinf(u) | isnan(u) | abs(imag(z))>pi/2;
m = ~(u==1 | p);
y = z;
y(m) = (u(m)-1).*(z(m)./log(u(m)));
y(p) = u(p)-1;