use kahan summation for numerical stability

Signed-off-by: darshanime <deathbullet@gmail.com>

promql/functions.go

@@ -439,11 +439,14 @@ func funcMinOverTime(vals []parser.Value, args parser.Expressions, enh *EvalNode
 // === sum_over_time(Matrix parser.ValueTypeMatrix) Vector ===
 func funcSumOverTime(vals []parser.Value, args parser.Expressions, enh *EvalNodeHelper) Vector {
 	return aggrOverTime(vals, enh, func(values []Point) float64 {
-		var sum float64
+		var sum, c float64
 		for _, v := range values {
-			sum += v.V
+			sum, c = kahanSummationIter(v.V, sum, c)
 		}
-		return sum
+		if math.IsInf(sum, 0) {
+			return sum
+		}
+		return sum + c
 	})
 }
 
@@ -675,23 +678,52 @@ func funcTimestamp(vals []parser.Value, args parser.Expressions, enh *EvalNodeHe
 	return enh.Out
 }
 
+func kahanSummation(samples []float64) float64 {
+	sum, c := 0.0, 0.0
+
+	for _, v := range samples {
+		sum, c = kahanSummationIter(v, sum, c)
+	}
+	return sum + c
+}
+
+func kahanSummationIter(v, sum, c float64) (float64, float64) {
+	t := sum + v
+	// using Neumaier improvement, swap if next term larger than sum
+	if math.Abs(sum) >= math.Abs(v) {
+		c += (sum - t) + v
+	} else {
+		c += (v - t) + sum
+	}
+	sum = t
+	return sum, c
+}
+
 // linearRegression performs a least-square linear regression analysis on the
 // provided SamplePairs. It returns the slope, and the intercept value at the
 // provided time.
 func linearRegression(samples []Point, interceptTime int64) (slope, intercept float64) {
 	var (
-		n            float64
-		sumX, sumY   float64
-		sumXY, sumX2 float64
+		n          float64
+		sumX, cX   float64
+		sumY, cY   float64
+		sumXY, cXY float64
+		sumX2, cX2 float64
 	)
 	for _, sample := range samples {
-		x := float64(sample.T-interceptTime) / 1e3
 		n += 1.0
-		sumY += sample.V
-		sumX += x
-		sumXY += x * sample.V
-		sumX2 += x * x
+		x := float64(sample.T-interceptTime) / 1e3
+		sumX, cX = kahanSummationIter(x, sumX, cX)
+		sumY, cY = kahanSummationIter(sample.V, sumY, cY)
+		sumXY, cXY = kahanSummationIter(x*sample.V, sumXY, cXY)
+		sumX2, cX2 = kahanSummationIter(x*x, sumX2, cX2)
 	}
+
+	sumX = sumX + cX
+	sumY = sumY + cY
+	sumXY = sumXY + cXY
+	sumX2 = sumX2 + cX2
+
 	covXY := sumXY - sumX*sumY/n
 	varX := sumX2 - sumX*sumX/n
 

promql/functions_test.go

@@ -15,6 +15,7 @@ package promql
 
 import (
 	"context"
+	"math"
 	"testing"
 	"time"
 
@@ -71,3 +72,9 @@ func TestFunctionList(t *testing.T) {
 		require.True(t, ok, "function %s exists in parser package, but not in promql package", i)
 	}
 }
+
+func TestKahanSummation(t *testing.T) {
+	vals := []float64{1.0, math.Pow(10, 100), 1.0, -1 * math.Pow(10, 100)}
+	expected := 2.0
+	require.Equal(t, expected, kahanSummation(vals))
+}