wolffd@0: function result = logsum(logv)
wolffd@0: 
wolffd@0: len = length(logv);
wolffd@0: if (len<2);
wolffd@0:   error('Subroutine logsum cannot sum less than 2 terms.');
wolffd@0: end;
wolffd@0: 
wolffd@0: % First two terms
wolffd@0: if (logv(2)<logv(1)),
wolffd@0:   result = logv(1) + log( 1 + exp( logv(2)-logv(1) ) );
wolffd@0: else,
wolffd@0:   result = logv(2) + log( 1 + exp( logv(1)-logv(2) ) );
wolffd@0: end;
wolffd@0: 
wolffd@0: % Remaining terms
wolffd@0: for (i=3:len),
wolffd@0:   term = logv(i);
wolffd@0:   if (result<term),
wolffd@0:     result = term   + log( 1 + exp( result-term ) );
wolffd@0:   else,
wolffd@0:     result = result + log( 1 + exp( term-result ) );
wolffd@0:   end;    
wolffd@0: end