Add histogram_quantile() PromQL query for percentiles

The PromQL layer only had a fixed-width histogram() (per-bucket
counts), with no way to ask for p50/p95/p99. Adds
histogram_quantile(p, metric[range]), interpolating over each
series's retained raw sample values (linear interpolation between
closest ranks) rather than building a separate pre-aggregated
bucketed-histogram representation — exact for the data actually kept.
An optional bracketed range narrows the computation to a trailing
window, mirroring rate()'s range selector.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01TeyhRfMYjMditek6geDRoQ

src/engine/snapshot/mod.rs

@@ -1287,10 +1287,76 @@ impl EngineSnapshot {
                     })
                     .collect()
             }
+            Some(AggOp::Quantile) => {
+                let quantile = parsed
+                    .quantile
+                    .expect("histogram_quantile() always carries a quantile");
+                // The bracketed range, if any, narrows each series to its trailing window
+                // (mirroring rate()'s range selector) rather than the full panel window —
+                // interpolating over whatever raw values are retained for that series, since
+                // this engine doesn't keep a separate pre-aggregated bucketed-histogram
+                // representation to compute quantiles from.
+                let range_ms = parsed.range.map(|d| d.as_millis() as i64);
+                let last_ts = series_points
+                    .values()
+                    .flat_map(|points| points.iter().map(|(ts, _)| *ts))
+                    .max()
+                    .unwrap_or(0);
+
+                let mut groups: BTreeMap<BTreeMap<String, String>, Vec<f64>> = BTreeMap::new();
+                for (labels, points) in &series_points {
+                    let group_key = project_labels(labels, &parsed.by_labels);
+                    let window_start = match range_ms {
+                        Some(range_ms) => {
+                            let window_end = points.last().map(|(ts, _)| *ts).unwrap_or(0);
+                            Some(window_end - range_ms)
+                        }
+                        None => None,
+                    };
+                    groups.entry(group_key).or_default().extend(
+                        points
+                            .iter()
+                            .filter(|(ts, _)| window_start.is_none_or(|start| *ts >= start))
+                            .map(|(_, value)| *value),
+                    );
+                }
+
+                groups
+                    .into_iter()
+                    .filter_map(|(labels, mut values)| {
+                        let quantile_value = interpolated_quantile(&mut values, quantile)?;
+                        Some(PromqlSeriesDto {
+                            labels,
+                            points: vec![(last_ts, quantile_value)],
+                        })
+                    })
+                    .collect()
+            }
         })
     }
 }
 
+/// Linear interpolation between closest ranks over the raw retained sample values (the "R-7"
+/// method, matching numpy/Excel's default `PERCENTILE` behavior) — exact for the data actually
+/// kept, unlike a bucketed-histogram approximation.
+fn interpolated_quantile(values: &mut [f64], quantile: f64) -> Option<f64> {
+    if values.is_empty() {
+        return None;
+    }
+    values.sort_by(|a, b| a.partial_cmp(b).unwrap());
+    if values.len() == 1 {
+        return Some(values[0]);
+    }
+    let rank = quantile * (values.len() - 1) as f64;
+    let lower = rank.floor() as usize;
+    let upper = rank.ceil() as usize;
+    if lower == upper {
+        return Some(values[lower]);
+    }
+    let fraction = rank - lower as f64;
+    Some(values[lower] + (values[upper] - values[lower]) * fraction)
+}
+
 // --- sort: sort-helper index maintenance ---
 impl EngineSnapshot {
     pub(super) fn build_sort_helper_index(