From c2ac89dcb6960742e0fa734405d8621adbbac2ce Mon Sep 17 00:00:00 2001
From: sukibaby <163092272+sukibaby@users.noreply.github.com>
Date: Sun, 26 May 2024 15:54:00 -0700
Subject: [PATCH] Improve calc_mean function

A problem with using std::accumulate to calculate the mean is that small numbers get rounded down to zero when dealing with floating point numbers. This is solved by implementing the Kahan summation algorithm (https://en.wikipedia.org/wiki/Kahan_summation_algorithm).
---
 src/RageUtil.cpp | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/src/RageUtil.cpp b/src/RageUtil.cpp
index bf1b883967..04c845db8f 100644
--- a/src/RageUtil.cpp
+++ b/src/RageUtil.cpp
@@ -1190,9 +1190,21 @@ void SortRStringArray( std::vector<RString> &arrayRStrings, const bool bSortAsce
 			bSortAscending?CompareRStringsAsc:CompareRStringsDesc );
 }
 
-float calc_mean( const float *pStart, const float *pEnd )
+float calc_mean(const float* pStart, const float* pEnd)
 {
-	return std::accumulate( pStart, pEnd, 0.f ) / std::distance( pStart, pEnd );
+	/* The Kahan summation algorithm is used here to prevent
+	 * situations where the low order bits may be lost.
+	 * https://en.wikipedia.org/wiki/Kahan_summation_algorithm */
+	float sum = 0.0f;
+	float c = 0.0f;
+	for (const float* p = pStart; p != pEnd; ++p)
+	{
+		float y = *p - c;
+		float t = sum + y;
+		c = (t - sum) - y;
+		sum = t;
+	}
+	return sum / (pEnd - pStart);
 }
 
 float calc_stddev( const float *pStart, const float *pEnd, bool bSample )