% C = LOGSUMEXP(A, B)
%
% Computes LOG(EXP(A)+EXP(B)), with attention to numerical precision.  
%
% A naive computation would overflow if A or B were very positive, or if
% both A and B were very negative, yielding +Inf or -Inf respectively.
% More subtly, if A or B were moderately negative, the naive computation
% would underflow (lose relative precision).  Moderately negative means
% that exp(A) or exp(B) is small but does not itself cause underflow.
%
% Example of positive overflow:
% >> log(exp(-1e5)+exp(1e5))
% ans = Inf
% >> logsumexp(-1e5,1e5)
% ans = 100000
%
% Example of negative overflow:
% >> log(exp(-1e5)+exp(-1e5))
% ans = -Inf
% >> logsumexp(-1e5,1e5)
% ans = -9.9999e+04
%
% Example of underflow:
% >> log(exp(-1e2)+exp(0))
% ans = 0
% >> logsumexp(-1e2,0)
% ans = 3.7201e-44

function c = logsumexp(a, b)

% ensure that A is the larger argument
mask = (a < b);
s = a;
a(mask) = b(mask);
b(mask) = s(mask);

% factor out exp(A), computing A + log(1 + exp(B-A))
c = a + log1p(exp(b-a));