logsumexp <- function(l) {
    n <- length(l);
    L <- sort(l, decreasing = TRUE);
    S <- rep(L[1], n);
    for (k in 1:(n-1)) {
        S[k+1] <- max(L[k+1], S[k]) + log1p(exp(-abs(L[k+1] - S[k]))); }
    S[n]; }

#Generate large number of small log-probabilities
set.seed(1);
n <- 10^6;
l <- rnorm(n, -3000, 1);

#Calculate logsumexp using above function and package
library(matrixStats);
l1 <- logsumexp(l);
l2 <- matrixStats::logSumExp(l);

print(l1, digits = 20)
[1] -2985.6845568559206
print(l2, digits = 20)
[1] -2985.6845568559038

#Test calculation speed of functions (specific to my PC - not reproducible)
TIME1 <- system.time(logsumexp(l));
TIME2 <- system.time(matrixStats::logSumExp(l));

TIME1[3]
elapsed
    0.7

TIME2[3]
elapsed
   0.05