fix: tdigest Quantile swaps interpolation weights and adds where the reference subtracts

The reference MergingDigest pairs mean[i] with the distance to the right
edge of the bracket and mean[i+1] with the distance to the left edge. We
passed them the other way round, so within a bracket a rising rank pulls
the estimate toward the lower centroid and Quantile is not monotonic.

The right tail branch used max + (...) where the reference uses max - (...),
so it can return a value above MaxValue. Reachable through DecodeDouble,
which accepts a tail centroid of weight above 1.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Author: jaideeppyne

tdigest/double.go

@@ -376,7 +376,7 @@ func (d *Double) Quantile(rank float64) (float64, error) {
 
 	lastWeight := float64(d.centroids[len(d.centroids)-1].weight)
 	if lastWeight > 1 && float64(d.centroidsWeight)-weight <= lastWeight/2.0 {
-		return d.max + (float64(d.centroidsWeight)-weight-1.0)/(lastWeight/2.0-1.0)*(d.max-d.centroids[len(d.centroids)-1].mean), nil
+		return d.max - (float64(d.centroidsWeight)-weight-1.0)/(lastWeight/2.0-1.0)*(d.max-d.centroids[len(d.centroids)-1].mean), nil
 	}
 
 	// interpolate between extremes
@@ -401,7 +401,7 @@ func (d *Double) Quantile(rank float64) (float64, error) {
 			}
 			w1 := weight - weightSoFar - leftWeight
 			w2 := weightSoFar + dw - weight - rightWeight
-			return weightedAverage(d.centroids[i].mean, w1, d.centroids[i+1].mean, w2), nil
+			return weightedAverage(d.centroids[i].mean, w2, d.centroids[i+1].mean, w1), nil
 		}
 		weightSoFar += dw
 	}

tdigest/double_test.go

@@ -19,6 +19,7 @@ package tdigest
 
 import (
 	"math"
+	"math/rand"
 	"testing"
 
 	"github.com/stretchr/testify/assert"
@@ -1023,3 +1024,89 @@ func TestDouble_SerializedSizeBytes(t *testing.T) {
 		assert.Equal(t, 1632, sizeWithoutBuffer)
 	})
 }
+
+func TestDouble_QuantileIsMonotonic(t *testing.T) {
+	distributions := map[string]func(r *rand.Rand) float64{
+		"gaussian":  func(r *rand.Rand) float64 { return r.NormFloat64() },
+		"uniform":   func(r *rand.Rand) float64 { return r.Float64() * 1000 },
+		"lognormal": func(r *rand.Rand) float64 { return math.Exp(r.NormFloat64()) },
+	}
+
+	for name, next := range distributions {
+		t.Run(name, func(t *testing.T) {
+			r := rand.New(rand.NewSource(1))
+			for _, k := range []uint16{10, 20, 50, 100, 200} {
+				for trial := 0; trial < 40; trial++ {
+					sketch, err := NewDouble(k)
+					assert.NoError(t, err)
+					n := 100 + r.Intn(50000)
+					for i := 0; i < n; i++ {
+						assert.NoError(t, sketch.Update(next(r)))
+					}
+
+					previousRank := 0.0
+					previousQuantile, err := sketch.Quantile(0)
+					assert.NoError(t, err)
+					for i := 1; i <= 1000; i++ {
+						rank := float64(i) / 1000.0
+						quantile, err := sketch.Quantile(rank)
+						assert.NoError(t, err)
+						assert.GreaterOrEqualf(t, quantile, previousQuantile,
+							"k=%d n=%d: Quantile(%v)=%v is below Quantile(%v)=%v",
+							k, n, rank, quantile, previousRank, previousQuantile)
+						previousRank = rank
+						previousQuantile = quantile
+					}
+				}
+			}
+		})
+	}
+}
+
+func TestDouble_QuantileIsWithinMinAndMax(t *testing.T) {
+	t.Run("From updates", func(t *testing.T) {
+		r := rand.New(rand.NewSource(2))
+		for trial := 0; trial < 100; trial++ {
+			sketch, err := NewDouble(DefaultK)
+			assert.NoError(t, err)
+			n := 10 + r.Intn(20000)
+			for i := 0; i < n; i++ {
+				assert.NoError(t, sketch.Update(r.NormFloat64()))
+			}
+			minValue, err := sketch.MinValue()
+			assert.NoError(t, err)
+			maxValue, err := sketch.MaxValue()
+			assert.NoError(t, err)
+			for i := 0; i <= 1000; i++ {
+				quantile, err := sketch.Quantile(float64(i) / 1000.0)
+				assert.NoError(t, err)
+				assert.GreaterOrEqual(t, quantile, minValue)
+				assert.LessOrEqual(t, quantile, maxValue)
+			}
+		}
+	})
+
+	// Update and Merge always leave both tail centroids as singletons, so the
+	// tail branches of Quantile are only reachable through deserialization.
+	t.Run("Decoded with heavy tail centroids", func(t *testing.T) {
+		sketch, err := newDoubleFromInternalStates(false, DefaultK, 0, 40,
+			[]doublePrecisionCentroid{{mean: 10, weight: 100}, {mean: 20, weight: 100}, {mean: 30, weight: 100}},
+			300, nil)
+		assert.NoError(t, err)
+		bytes, err := EncodeDouble(sketch, false)
+		assert.NoError(t, err)
+		decoded, err := DecodeDouble(bytes)
+		assert.NoError(t, err)
+
+		quantile, err := decoded.Quantile(0.9)
+		assert.NoError(t, err)
+		assert.InDelta(t, 34.081632653061224, quantile, 1e-12)
+
+		for i := 0; i <= 1000; i++ {
+			quantile, err := decoded.Quantile(float64(i) / 1000.0)
+			assert.NoError(t, err)
+			assert.GreaterOrEqual(t, quantile, 0.0)
+			assert.LessOrEqual(t, quantile, 40.0)
+		}
+	})
+}