function y = log1p(x)
%LOG1P log(1+x) without underflow
%
% This is for use in Octave. Matlab already has this function, and its version
% is probably better designed and tested.

% Iain Murray, November 2007

% Masks for three different methods:
nan_mask = isnan(x);
small = (abs(x) < 1e-5); % Slightly arbitrary
naive_mask = ~(nan_mask | small);

% x is fairly big. Just do simple thing
y = zeros(size(x));
y(naive_mask) = log(1 + x(naive_mask));

% x is fairly small. Use the trick from:
% http://docs.sun.com/source/806-3568/ncg_goldberg.html
% But I'm just justifying the 'small' range empirically. I haven't bothered to
% work out the exact properties of this and can imagine getting burnt in
% the future. Be careful!
xs = x(small);
p1 = xs + 1;
careful = xs.*log(p1)./(p1-1);
tiny_mask = (p1==1);
careful(tiny_mask) = xs(tiny_mask);
y(small) = careful;

y(nan_mask) = NaN;