A vectorized, numerically-stable implementation of logsumexp on top of Accelerate.framework

logsumexp.c

#include <Accelerate/Accelerate.h>

static double vsum(const double *vector, size_t length) {
	double *const ones = (double *)malloc(sizeof(double) * length);
	const double one = 1;
	vDSP_vfillD(&one, ones, 1, length);

	double sum = 0;
	vDSP_dotprD(ones, 1, vector, 1, &sum, length);

	free(ones);

	return sum;
}

static double vlogsumexp(const double *vector, size_t length) {
	if (vector == NULL || length == 0)
		return 0;

	if (length > INT_MAX) {
		printf("Warning: elements at indices greater than or equal to %ud will not be processed", INT_MAX);
		length = INT_MAX;
	}

	double maxExp = 0;
	vDSP_maxvD(vector, 1, &maxExp, length);

	double *const difference = (double *)malloc(sizeof(double) * length);
	double negativeMaxExp = -maxExp;
	vDSP_vsaddD(vector, 1, &negativeMaxExp, difference, 1, length);

	const int intLength = (int)length;
	vvexp(difference, difference, &intLength);

	const double sum = vsum(difference, length);
	free(difference);

	return log(sum) + maxExp;
}