diff --git a/quantile.go b/quantile.go
index 8eb5534..96f5f1a 100644
--- a/quantile.go
+++ b/quantile.go
@@ -95,17 +95,16 @@ func (s *Snapshot) ValueAtQuantile(q float64) float64 {
 	}
 
 	// Step 2: Estimate density at each boundary by averaging neighbors.
-	// boundaryDensity has length n+1 (one per boundary).
+	// boundaryDensity has length n+1 (one per boundary). At the outer edges
+	// there is no neighbor on one side, so we use the adjacent bucket's
+	// density directly rather than averaging with zero — otherwise the
+	// rightmost-bucket interpolation gets biased low (which matters a lot
+	// for p99 in long-tailed distributions).
 	boundaryDensity := make([]float64, n+1)
-	for i := range n {
-		switch i {
-		case 0:
-			boundaryDensity[i] = avgDensity[0]
-		case n:
-			boundaryDensity[i] = avgDensity[n-1]
-		default:
-			boundaryDensity[i] = (avgDensity[i-1] + avgDensity[i]) / 2.0
-		}
+	boundaryDensity[0] = avgDensity[0]
+	boundaryDensity[n] = avgDensity[n-1]
+	for i := 1; i < n; i++ {
+		boundaryDensity[i] = (avgDensity[i-1] + avgDensity[i]) / 2.0
 	}
 
 	// Step 3: Walk buckets to find which one contains the target rank,
@@ -220,13 +219,10 @@ func (s *Snapshot) ValuesAtQuantiles(qs []float64) []float64 {
 		}
 	}
 	boundaryDensity := make([]float64, n+1)
-	for i := range n {
-		switch i {
-		case 0:
-			boundaryDensity[i] = avgDensity[0]
-		default:
-			boundaryDensity[i] = (avgDensity[i-1] + avgDensity[i]) / 2.0
-		}
+	boundaryDensity[0] = avgDensity[0]
+	boundaryDensity[n] = avgDensity[n-1]
+	for i := 1; i < n; i++ {
+		boundaryDensity[i] = (avgDensity[i-1] + avgDensity[i]) / 2.0
 	}
 
 	// Single-pass bucket walk: process all quantiles whose rank falls
diff --git a/quantile_live.go b/quantile_live.go
new file mode 100644
index 0000000..b3bbb03
--- /dev/null
+++ b/quantile_live.go
@@ -0,0 +1,192 @@
+// Copyright 2026 The Cockroach Authors.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+
+package goodhistogram
+
+// ValuesAtQuantilesInto writes estimated values at the given quantiles into dst and
+// returns dst[:len(qs)]. It reads live atomic counters directly without
+// materializing a Snapshot.
+//
+// The result reflects the same eventual consistency Snapshot() already
+// accepts: counters are read independently and may observe a slightly
+// inconsistent total. The inconsistency window here is wider than
+// Snapshot+ValuesAtQuantiles because counters are loaded twice (once to
+// total, once during the walk); for monotonic counters the only effect is
+// that cumulative bucket count may exceed the precomputed total at the
+// tail, which is harmless.
+//
+// qs MUST be sorted in ascending order. dst must have cap >= len(qs); pass a
+// stack-backed slice (e.g. var buf [4]float64; h.ValuesAtQuantilesInto(buf[:0], qs))
+// to make the call fully alloc-free.
+func (h *Histogram) ValuesAtQuantilesInto(dst, qs []float64) []float64 {
+	dst = dst[:len(qs)]
+	if len(qs) == 0 {
+		return dst
+	}
+	cfg := h.cfg
+	n := len(h.counts)
+
+	// Pass 1: load scalars and sum the in-range total.
+	zeroCount := h.ZeroCount.Load()
+	underflow := h.Underflow.Load()
+	overflow := h.Overflow.Load()
+	var inRange uint64
+	for i := 0; i < n; i++ {
+		inRange += h.counts[i].Load()
+	}
+	total := zeroCount + underflow + overflow + inRange
+
+	if total == 0 {
+		for i := range dst {
+			dst[i] = 0
+		}
+		return dst
+	}
+
+	fTotal := float64(total)
+	belowLo := float64(zeroCount + underflow)
+	fInRange := float64(inRange)
+
+	// Classify each quantile. Since qs is sorted ascending and rank = q*total
+	// is monotonic, low-edges come first, walk-eligible middle next, then
+	// high-edges. We resolve edges directly into dst and remember the walk
+	// range as [walkStart, walkEnd).
+	walkStart := len(qs)
+	walkEnd := len(qs)
+	for i, q := range qs {
+		rank := q * fTotal
+		switch {
+		case rank <= 0:
+			if zeroCount+underflow > 0 {
+				dst[i] = cfg.lo
+			} else {
+				dst[i] = cfg.hi
+				for j := 0; j < n; j++ {
+					if h.counts[j].Load() > 0 {
+						dst[i] = cfg.boundaries[j]
+						break
+					}
+				}
+			}
+		case rank >= fTotal:
+			if overflow > 0 {
+				dst[i] = cfg.hi
+			} else {
+				dst[i] = cfg.lo
+				for j := n - 1; j >= 0; j-- {
+					if h.counts[j].Load() > 0 {
+						dst[i] = cfg.boundaries[j+1]
+						break
+					}
+				}
+			}
+		case rank <= belowLo:
+			dst[i] = cfg.lo
+		case rank-belowLo > fInRange:
+			dst[i] = cfg.hi
+		default:
+			if i < walkStart {
+				walkStart = i
+			}
+			walkEnd = i + 1
+		}
+	}
+
+	if walkStart >= walkEnd {
+		return dst
+	}
+
+	// Pass 2: forward walk with a 3-count sliding window. We re-load each
+	// bucket once (peeking ahead by 1) so we have prev/curr/next counts for
+	// computing boundary densities on the fly — no scratch slices.
+	//
+	// boundaryDensity[i]   = (avgDensity[i-1] + avgDensity[i]) / 2
+	// boundaryDensity[i+1] = (avgDensity[i]   + avgDensity[i+1]) / 2
+	// At the outer edges there is no neighbor on one side, so we use the
+	// adjacent bucket's density directly (dL = currD at i==0, dR = currD
+	// at i==n-1) instead of averaging with zero, matching Snapshot's
+	// ValuesAtQuantiles.
+
+	var prevCount, currCount, nextCount uint64
+	var prevW, currW, nextW float64
+
+	currCount = h.counts[0].Load()
+	currW = cfg.boundaries[1] - cfg.boundaries[0]
+	if n > 1 {
+		nextCount = h.counts[1].Load()
+		nextW = cfg.boundaries[2] - cfg.boundaries[1]
+	}
+
+	var cumCount float64
+	wi := walkStart
+
+	for i := 0; i < n && wi < walkEnd; i++ {
+		fc := float64(currCount)
+		nextCum := cumCount + fc
+
+		// Process all walk-eligible quantiles whose adjusted rank falls
+		// in [cumCount, nextCum].
+		for wi < walkEnd {
+			adjRank := qs[wi]*fTotal - belowLo
+			if nextCum < adjRank {
+				break
+			}
+			localRank := adjRank - cumCount
+			lo := cfg.boundaries[i]
+			if currW <= 0 || fc == 0 {
+				dst[wi] = lo
+				wi++
+				continue
+			}
+			currD := fc / currW
+			var dL, dR float64
+			if i == 0 {
+				dL = currD
+			} else {
+				var prevD float64
+				if prevW > 0 && prevCount > 0 {
+					prevD = float64(prevCount) / prevW
+				}
+				dL = (prevD + currD) / 2.0
+			}
+			if i == n-1 {
+				dR = currD
+			} else {
+				var nextD float64
+				if nextW > 0 && nextCount > 0 {
+					nextD = float64(nextCount) / nextW
+				}
+				dR = (currD + nextD) / 2.0
+			}
+			dst[wi] = trapezoidalSolve(lo, currW, fc, dL, dR, localRank)
+			wi++
+		}
+
+		cumCount = nextCum
+
+		// Slide window forward.
+		prevCount, prevW = currCount, currW
+		currCount, currW = nextCount, nextW
+		if i+2 < n {
+			nextCount = h.counts[i+2].Load()
+			nextW = cfg.boundaries[i+3] - cfg.boundaries[i+2]
+		} else {
+			nextCount = 0
+			nextW = 0
+		}
+	}
+
+	// Safety net: any walk-eligible quantiles not yet resolved (shouldn't
+	// happen with monotonic counters, but counters can grow between the two
+	// passes, so the walk may technically fall short).
+	for ; wi < walkEnd; wi++ {
+		dst[wi] = cfg.boundaries[n]
+	}
+
+	return dst
+}
diff --git a/quantile_live_test.go b/quantile_live_test.go
new file mode 100644
index 0000000..2845fc9
--- /dev/null
+++ b/quantile_live_test.go
@@ -0,0 +1,234 @@
+// Copyright 2026 The Cockroach Authors.
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+
+package goodhistogram
+
+import (
+	"fmt"
+	"math/rand"
+	"sync"
+	"sync/atomic"
+	"testing"
+	"time"
+)
+
+// TestValuesAtQuantilesIntoAgreesWithSnapshot checks that
+// ValuesAtQuantilesInto produces exactly the same numbers as
+// Snapshot().ValuesAtQuantiles() across distributions. Equality is
+// bit-for-bit: both paths feed identical arguments to trapezoidalSolve in
+// the same order.
+func TestValuesAtQuantilesIntoAgreesWithSnapshot(t *testing.T) {
+	qs := []float64{0.0, 0.001, 0.01, 0.5, 0.75, 0.9, 0.95, 0.99, 0.999, 1.0}
+
+	for _, dist := range distributions {
+		t.Run(dist.name, func(t *testing.T) {
+			rng := rand.New(rand.NewSource(42))
+			vals := dist.genFn(rng, 100_000)
+
+			h := newGoodHist()
+			for _, v := range vals {
+				h.Record(int64(v))
+			}
+
+			snap := h.Snapshot()
+			want := snap.ValuesAtQuantiles(qs)
+
+			var buf [16]float64
+			got := h.ValuesAtQuantilesInto(buf[:0], qs)
+
+			for i, q := range qs {
+				if got[i] != want[i] {
+					t.Errorf("q=%g: ValuesAtQuantilesInto=%g, ValuesAtQuantiles=%g (diff=%g)",
+						q, got[i], want[i], got[i]-want[i])
+				}
+			}
+		})
+	}
+}
+
+// TestValuesAtQuantilesIntoConcurrentWithRecord runs Record and
+// ValuesAtQuantilesInto concurrently to lock in the lock-free contract
+// under -race: any future change that introduces a data race (e.g.
+// sharing scratch state across callers) will be caught here.
+func TestValuesAtQuantilesIntoConcurrentWithRecord(t *testing.T) {
+	h := newGoodHist()
+	qs := []float64{0.5, 0.9, 0.99}
+
+	// Pre-seed so readers don't observe a transient total==0 (which
+	// returns zeros, not in-range values). The race detector is the
+	// primary signal; the range assertion is just a sanity check.
+	seedRng := rand.New(rand.NewSource(7))
+	for i := 0; i < 1000; i++ {
+		h.Record(int64(benchLo + seedRng.Float64()*benchRange))
+	}
+
+	const writers = 4
+	const readers = 4
+	var stop atomic.Bool
+	var wg sync.WaitGroup
+
+	for w := 0; w < writers; w++ {
+		wg.Add(1)
+		go func(seed int64) {
+			defer wg.Done()
+			rng := rand.New(rand.NewSource(seed))
+			for !stop.Load() {
+				h.Record(int64(benchLo + rng.Float64()*benchRange))
+			}
+		}(int64(w + 1))
+	}
+
+	for r := 0; r < readers; r++ {
+		wg.Add(1)
+		go func() {
+			defer wg.Done()
+			var buf [4]float64
+			for !stop.Load() {
+				got := h.ValuesAtQuantilesInto(buf[:0], qs)
+				for i, v := range got {
+					if v < benchLo || v > benchHi {
+						t.Errorf("q=%g: out-of-range value %g", qs[i], v)
+						return
+					}
+				}
+			}
+		}()
+	}
+
+	time.Sleep(100 * time.Millisecond)
+	stop.Store(true)
+	wg.Wait()
+}
+
+// TestValuesAtQuantilesIntoEdges checks zero-count and edge-only inputs.
+func TestValuesAtQuantilesIntoEdges(t *testing.T) {
+	t.Run("empty histogram", func(t *testing.T) {
+		h := newGoodHist()
+		var buf [4]float64
+		got := h.ValuesAtQuantilesInto(buf[:0], []float64{0.5, 0.99})
+		for i, v := range got {
+			if v != 0 {
+				t.Errorf("empty histogram q[%d]: got %g, want 0", i, v)
+			}
+		}
+	})
+
+	t.Run("only underflow", func(t *testing.T) {
+		// Use values clearly below lo's octave to guarantee underflow.
+		h := newGoodHist()
+		h.Record(1)
+		h.Record(10)
+		var buf [3]float64
+		got := h.ValuesAtQuantilesInto(buf[:0], []float64{0.0, 0.5, 1.0})
+		snap := h.Snapshot()
+		want := snap.ValuesAtQuantiles([]float64{0.0, 0.5, 1.0})
+		for i, v := range got {
+			if v != want[i] {
+				t.Errorf("only-underflow q[%d]: ValuesAtQuantilesInto=%g, ValuesAtQuantiles=%g", i, v, want[i])
+			}
+		}
+	})
+
+	t.Run("only overflow", func(t *testing.T) {
+		// Use values clearly above hi's octave to guarantee overflow.
+		h := newGoodHist()
+		h.Record(int64(benchHi * 4))
+		h.Record(int64(benchHi * 8))
+		var buf [3]float64
+		got := h.ValuesAtQuantilesInto(buf[:0], []float64{0.0, 0.5, 1.0})
+		snap := h.Snapshot()
+		want := snap.ValuesAtQuantiles([]float64{0.0, 0.5, 1.0})
+		for i, v := range got {
+			if v != want[i] {
+				t.Errorf("only-overflow q[%d]: ValuesAtQuantilesInto=%g, ValuesAtQuantiles=%g", i, v, want[i])
+			}
+		}
+	})
+
+	t.Run("empty qs", func(t *testing.T) {
+		h := newGoodHist()
+		h.Record(1000)
+		got := h.ValuesAtQuantilesInto(nil, nil)
+		if len(got) != 0 {
+			t.Errorf("got len=%d, want 0", len(got))
+		}
+	})
+}
+
+// TestValuesAtQuantilesIntoAllocFree verifies zero allocations when dst has cap.
+func TestValuesAtQuantilesIntoAllocFree(t *testing.T) {
+	h := newGoodHist()
+	rng := rand.New(rand.NewSource(42))
+	for i := 0; i < 10_000; i++ {
+		h.Record(int64(benchLo + rng.Float64()*benchRange))
+	}
+	qs := []float64{0.5, 0.99}
+	var buf [4]float64
+
+	allocs := testing.AllocsPerRun(100, func() {
+		_ = h.ValuesAtQuantilesInto(buf[:0], qs)
+	})
+	if allocs != 0 {
+		t.Errorf("ValuesAtQuantilesInto allocated %v times per run, want 0", allocs)
+	}
+}
+
+// BenchmarkQueryPath compares the existing Snapshot+ValuesAtQuantiles path
+// against the new ValuesAtQuantilesInto path. Single-thread is the apples-to-apples
+// comparison since allocation/copy cost is what we're targeting.
+func BenchmarkQueryPath(b *testing.B) {
+	for _, nObs := range []int{1_000, 100_000} {
+		b.Run(fmt.Sprintf("n=%d", nObs), func(b *testing.B) {
+			h := newGoodHist()
+			rng := rand.New(rand.NewSource(42))
+			for i := 0; i < nObs; i++ {
+				h.Record(int64(benchLo + rng.Float64()*benchRange))
+			}
+			qs := []float64{0.5, 0.99}
+
+			b.Run("SnapshotValuesAtQuantiles", func(b *testing.B) {
+				b.ReportAllocs()
+				for i := 0; i < b.N; i++ {
+					snap := h.Snapshot()
+					_ = snap.ValuesAtQuantiles(qs)
+				}
+			})
+			b.Run("ValuesAtQuantilesInto", func(b *testing.B) {
+				b.ReportAllocs()
+				var buf [4]float64
+				for i := 0; i < b.N; i++ {
+					_ = h.ValuesAtQuantilesInto(buf[:0], qs)
+				}
+			})
+		})
+	}
+}
+
+func BenchmarkQueryPathThreeQuantiles(b *testing.B) {
+	h := newGoodHist()
+	rng := rand.New(rand.NewSource(42))
+	for i := 0; i < 10_000; i++ {
+		h.Record(int64(benchLo + rng.Float64()*benchRange))
+	}
+	qs := []float64{0.5, 0.95, 0.99}
+
+	b.Run("SnapshotValuesAtQuantiles", func(b *testing.B) {
+		b.ReportAllocs()
+		for i := 0; i < b.N; i++ {
+			snap := h.Snapshot()
+			_ = snap.ValuesAtQuantiles(qs)
+		}
+	})
+	b.Run("ValuesAtQuantilesInto", func(b *testing.B) {
+		b.ReportAllocs()
+		var buf [4]float64
+		for i := 0; i < b.N; i++ {
+			_ = h.ValuesAtQuantilesInto(buf[:0], qs)
+		}
+	})
+}