perf_hooks: implement qrde analysis support in Histogram

Implements quantile-respectful density estimate calcuation
on Histogram. Helps with tail-focused latency analysis
without retaining raw samples.

Baking this directly into Node.js, based on a 1 million
sample, 1k bin workload, this impl is roughly 150-325x
faster than performing the equivalent in Rscript.

Signed-off-by: James M Snell <jasnell@gmail.com>
Assisted-by: Opencode
PR-URL: https://github.com/nodejs/node/pull/65806
Reviewed-By: Matteo Collina <matteo.collina@gmail.com>

benchmark/perf_hooks/histogram-qrde.js

@@ -0,0 +1,37 @@
+'use strict';
+
+const common = require('../common.js');
+const { createHistogram } = require('perf_hooks');
+
+const bench = common.createBenchmark(main, {
+  n: [5],
+  bins: [100, 1000],
+  samples: [1e6],
+  unique: [100, 1000, 10000],
+  dequantize: ['none', 'hdr', 'all'],
+}, {
+  test: {
+    n: 1,
+    bins: 10,
+    samples: 100,
+    unique: 10,
+  },
+});
+
+async function main({ n, bins, samples, unique, dequantize }) {
+  const histogram = createHistogram();
+  const maximum = 1e12;
+
+  for (let i = 0; i < samples; i++) {
+    const index = i % unique;
+    const rank = unique === 1 ? 0 : index / (unique - 1);
+    histogram.record(Math.max(1, Math.round(maximum ** rank)));
+  }
+
+  await histogram.qrde({ bins, dequantize });
+  bench.start();
+  for (let i = 0; i < n; i++) {
+    await histogram.qrde({ bins, dequantize });
+  }
+  bench.end(n);
+}

doc/api/perf_hooks.md

@@ -2467,6 +2467,82 @@ Returns the values at the specified percentiles, computed in a single
 efficient pass over the histogram data. More efficient than calling
 `histogram.percentile()` multiple times.
 
+### `histogram.qrde([options])`
+
+<!-- YAML
+added: REPLACEME
+-->
+
+* `options` {Object}
+  * `bins` {number} The number of equal-probability density bins to return.
+    Must be between 1 and 1000. Cannot be used with `probabilities`.
+    **Default:** `100`.
+  * `probabilities` {number\[]} Custom probability boundaries. The array must
+    contain between 2 and 1001 strictly increasing values, start with `0`, and
+    end with `1`. Cannot be used with `bins`.
+  * `dequantize` {string} Controls whether repeated bucket values are spread
+    deterministically over their equivalent-value ranges. May be `'none'`,
+    `'hdr'`, or `'all'`. **Default:** `'hdr'`.
+  * `cache` {boolean} When `true`, retains the expanded histogram snapshot for
+    reuse by subsequent calls with `cache: true`. The snapshot is invalidated
+    when the histogram is modified. **Default:** `false`.
+* Returns: {Promise} Fulfills with an {Object} containing:
+  * `probabilities` {Float64Array} The probability boundaries used by the
+    estimate.
+  * `quantiles` {Float64Array} The quantiles at the probability boundaries.
+  * `densities` {Float64Array} The density within each quantile interval.
+  * `count` {bigint} The number of values in the histogram snapshot.
+  * `bucketCount` {number} The number of occupied HDR buckets.
+  * `corrections` {number} The number of non-monotonic floating-point results
+    that were clamped to the preceding quantile.
+  * `dequantize` {string} The selected dequantization mode.
+
+Returns a quantile-respectful density estimate based on the Harrell-Davis
+quantile estimator. By default, `bins` generates equal probability boundaries.
+The `probabilities` option can instead focus the estimate on regions such as
+p90, p99, p99.9, and p99.99. The density for interval `i` contains probability
+mass `probabilities[i + 1] - probabilities[i]`. The histogram is snapshotted
+when the method is called. Snapshot expansion and the estimate are calculated
+in the libuv thread pool. Highly concentrated beta weights use a second-order
+asymptotic approximation to avoid numerical convergence loss at large sample
+counts.
+
+Setting `cache` to `true` avoids repeating snapshot capture and expansion when
+several estimates are requested from an unchanged histogram. The retained
+snapshot uses memory proportional to the number of occupied HDR buckets and is
+released when the histogram is next modified.
+
+QRDE temporarily uses approximately one additional HDR count array plus 32
+bytes per occupied bucket. With `cache: true`, the expanded 32-byte-per-bucket
+snapshot remains allocated. The following estimates use `lowest: 1` and
+`highest: Number.MAX_SAFE_INTEGER` and exclude allocator and JavaScript object
+overhead:
+
+| `figures` | Histogram | Maximum expanded snapshot | Peak cache-miss QRDE |
+| --------- | --------: | ------------------------: | -------------------: |
+| 1         |   6.3 KiB |                    25 KiB |               31 KiB |
+| 2         |    47 KiB |                   188 KiB |              235 KiB |
+| 3         |   352 KiB |                   1.4 MiB |              1.7 MiB |
+| 4         |   5.0 MiB |                    20 MiB |               25 MiB |
+| 5         |    37 MiB |                   148 MiB |              185 MiB |
+
+The maximum snapshot column assumes every representable bucket is occupied.
+Lower `highest` values reduce histogram and temporary copy sizes. Concurrent
+calls that miss the cache each require their own temporary copy and expanded
+snapshot.
+
+HDR histograms aggregate observations into equivalent-value buckets. The
+`'hdr'` dequantization mode models repeated values in buckets wider than one
+unit as a continuous uniform distribution over the bucket resolution. This
+reduces density artifacts introduced by HDR quantization while preserving
+repeated unit-resolution values as point masses. The `'all'` mode also
+dequantizes repeated unit-resolution values. Use `'none'` to calculate the
+grouped Harrell-Davis estimator using bucket midpoints directly.
+
+An empty histogram returns the requested `probabilities` but produces empty
+`quantiles` and `densities` arrays. A non-dequantized interval whose quantile
+boundaries are equal has an infinite density.
+
 ### `histogram.reset()`
 
 <!-- YAML

lib/internal/histogram.js

@@ -33,9 +33,11 @@ const {
 
 const {
   validateArray,
+  validateBoolean,
   validateInteger,
   validateNumber,
   validateObject,
+  validateOneOf,
 } = require('internal/validators');
 
 const {
@@ -45,6 +47,7 @@ const {
 const kDestroy = Symbol('kDestroy');
 const kHandle = Symbol('kHandle');
 const kRecordable = Symbol('kRecordable');
+const kQrdeDequantizationModes = ['none', 'hdr', 'all'];
 
 const {
   kClone,
@@ -594,6 +597,73 @@ class Histogram {
     return map;
   }
 
+  /**
+   * Builds a quantile-respectful density estimate using Harrell-Davis
+   * quantiles. Values can be spread deterministically within their HDR
+   * equivalent-value ranges to reduce quantization artifacts.
+   * @param {{ bins?: number, probabilities?: number[],
+   *   dequantize?: 'none'|'hdr'|'all', cache?: boolean }} [options]
+   * @returns {Promise<{
+   *   probabilities: Float64Array,
+   *   quantiles: Float64Array,
+   *   densities: Float64Array,
+   *   count: bigint,
+   *   bucketCount: number,
+   *   corrections: number,
+   *   dequantize: 'none'|'hdr'|'all',
+   * }>}
+   */
+  qrde(options = kEmptyObject) {
+    if (!isHistogram(this))
+      throw new ERR_INVALID_THIS('Histogram');
+    validateObject(options, 'options');
+    const { bins, probabilities, dequantize = 'hdr', cache = false } = options;
+    if (bins !== undefined && probabilities !== undefined) {
+      throw new ERR_INVALID_ARG_VALUE(
+        'options', options, '"bins" and "probabilities" are mutually exclusive');
+    }
+
+    let boundaries;
+    if (probabilities === undefined) {
+      const binCount = bins ?? 100;
+      validateInteger(binCount, 'options.bins', 1, 1000);
+      boundaries = new Float64Array(binCount + 1);
+      for (let i = 0; i <= binCount; i++) boundaries[i] = i / binCount;
+    } else {
+      validateArray(probabilities, 'options.probabilities', 2);
+      const length = probabilities.length;
+      if (length > 1001) {
+        throw new ERR_OUT_OF_RANGE(
+          'options.probabilities.length', '>= 2 && <= 1001', length);
+      }
+
+      boundaries = new Float64Array(length);
+      let previous = -1;
+      for (let i = 0; i < length; i++) {
+        const probability = probabilities[i];
+        validateNumber(probability, `options.probabilities[${i}]`, 0, 1);
+        if (probability <= previous) {
+          throw new ERR_INVALID_ARG_VALUE(
+            'options.probabilities', probabilities, 'must be strictly increasing');
+        }
+        boundaries[i] = probability;
+        previous = probability;
+      }
+      if (boundaries[0] !== 0 || boundaries[length - 1] !== 1) {
+        throw new ERR_INVALID_ARG_VALUE(
+          'options.probabilities', probabilities, 'must start with 0 and end with 1');
+      }
+    }
+
+    validateOneOf(dequantize,
+                  'options.dequantize', kQrdeDequantizationModes);
+    validateBoolean(cache, 'options.cache');
+    let mode = 0;
+    if (dequantize === 'hdr') mode = 1;
+    else if (dequantize === 'all') mode = 2;
+    return this[kHandle]?.qrde(boundaries, mode, cache);
+  }
+
   /**
    * @returns {void}
    */

src/histogram-inl.h

@@ -33,9 +33,15 @@ void Histogram::UpdateEwma(double value) {
   }
 }
 
+void Histogram::InvalidateRecordedSnapshot() {
+  mutation_generation_++;
+  recorded_snapshot_cache_.reset();
+}
+
 void Histogram::Reset() {
   RwLock::ScopedWriteLock lock(mutex_);
   hdr_reset(histogram_.get());
+  InvalidateRecordedSnapshot();
   exceeds_ = 0;
   prev_ = 0;
   ewma_mean_ = 0;
@@ -49,8 +55,10 @@ double Histogram::Add(const Histogram& other) {
     exceeds_ += other.exceeds_;
     if (other.prev_ > prev_) prev_ = other.prev_;
     // hdr_add merges all bucket counts and total_count internally.
-    return static_cast<double>(
-        hdr_add(histogram_.get(), other.histogram_.get()));
+    const double dropped =
+        static_cast<double>(hdr_add(histogram_.get(), other.histogram_.get()));
+    InvalidateRecordedSnapshot();
+    return dropped;
   };
 
   // When adding a histogram to itself, a single write lock suffices.
@@ -146,8 +154,10 @@ bool Histogram::RecordCorrected(int64_t value, int64_t expected_interval) {
       hdr_record_corrected_value(histogram_.get(), value, expected_interval);
   if (!recorded)
     exceeds_++;
-  else
+  else {
+    InvalidateRecordedSnapshot();
     UpdateEwma(static_cast<double>(value));
+  }
   return recorded;
 }
 
@@ -156,8 +166,10 @@ bool Histogram::Record(int64_t value) {
   bool recorded = hdr_record_value(histogram_.get(), value);
   if (!recorded)
     exceeds_++;
-  else
+  else {
+    InvalidateRecordedSnapshot();
     UpdateEwma(static_cast<double>(value));
+  }
   return recorded;
 }
 
@@ -170,8 +182,10 @@ uint64_t Histogram::RecordDelta() {
     delta = time - prev_;
     if (!hdr_record_value(histogram_.get(), delta))
       exceeds_++;
-    else
+    else {
+      InvalidateRecordedSnapshot();
       UpdateEwma(static_cast<double>(delta));
+    }
   }
   prev_ = time;
   return delta;

src/histogram.cc

@@ -5,31 +5,42 @@
 #include "node_debug.h"
 #include "node_errors.h"
 #include "node_external_reference.h"
+#include "threadpoolwork-inl.h"
 #include "util.h"
 #include "v8-typed-array.h"
 
+#include <algorithm>
 #include <cmath>
+#include <cstring>
 #include <numbers>
 #include <vector>
 
 namespace node {
 
 using v8::Array;
+using v8::ArrayBuffer;
 using v8::BigInt;
 using v8::CFunction;
 using v8::Context;
+using v8::Exception;
 using v8::Float64Array;
 using v8::FunctionCallbackInfo;
 using v8::FunctionTemplate;
+using v8::Global;
+using v8::HandleScope;
 using v8::Integer;
 using v8::Isolate;
 using v8::Local;
 using v8::Map;
+using v8::Name;
+using v8::Null;
 using v8::Number;
 using v8::Object;
 using v8::ObjectTemplate;
+using v8::Promise;
 using v8::String;
 using v8::Uint32;
+using v8::Uint8Array;
 using v8::Value;
 
 template <typename T>
@@ -66,6 +77,8 @@ std::shared_ptr<Histogram> Histogram::Create(const Options& options) {
 
 void Histogram::MemoryInfo(MemoryTracker* tracker) const {
   tracker->TrackFieldWithSize("histogram", GetMemorySize());
+  tracker->TrackFieldWithSize("qrde_snapshot",
+                              GetCachedRecordedSnapshotMemorySize());
 }
 
 bool Histogram::IsCompatible(const Histogram& other) const {
@@ -78,6 +91,99 @@ bool Histogram::IsCompatible(const Histogram& other) const {
              other.histogram_->significant_figures;
 }
 
+Histogram::RecordedSnapshot Histogram::BuildRecordedSnapshot(
+    const hdr_histogram* histogram) {
+  RecordedSnapshot snapshot;
+  snapshot.total_count = histogram->total_count;
+  if (snapshot.total_count == 0) return snapshot;
+
+  const size_t occupied =
+      std::count_if(histogram->counts,
+                    histogram->counts + histogram->counts_len,
+                    [](int64_t count) { return count != 0; });
+  snapshot.buckets.reserve(occupied);
+  hdr_iter iter;
+  hdr_iter_recorded_init(&iter, histogram);
+  while (hdr_iter_next(&iter)) {
+    snapshot.buckets.push_back({
+        static_cast<double>(iter.median_equivalent_value),
+        static_cast<double>(
+            hdr_size_of_equivalent_value_range(histogram, iter.value)),
+        iter.count,
+        iter.cumulative_count,
+    });
+  }
+  return snapshot;
+}
+
+Histogram::RecordedSnapshotSource Histogram::GetRecordedSnapshotSource(
+    bool use_cache) const {
+  RecordedSnapshotSource source;
+  int64_t lowest;
+  int64_t highest;
+  int figures;
+  {
+    RwLock::ScopedReadLock lock(mutex_);
+    source.generation = mutation_generation_;
+    if (use_cache && recorded_snapshot_cache_) {
+      source.snapshot = recorded_snapshot_cache_;
+      source.cache_hit = true;
+      return source;
+    }
+    if (histogram_->total_count == 0) {
+      source.snapshot = std::make_shared<const RecordedSnapshot>();
+      return source;
+    }
+    lowest = histogram_->lowest_discernible_value;
+    highest = histogram_->highest_trackable_value;
+    figures = histogram_->significant_figures;
+  }
+
+  hdr_histogram* copy;
+  if (hdr_init(lowest, highest, figures, &copy) != 0) return source;
+  source.histogram.reset(copy);
+
+  RwLock::ScopedReadLock lock(mutex_);
+  source.generation = mutation_generation_;
+  if (use_cache && recorded_snapshot_cache_) {
+    source.histogram.reset();
+    source.snapshot = recorded_snapshot_cache_;
+    source.cache_hit = true;
+    return source;
+  }
+  if (histogram_->total_count == 0) {
+    source.histogram.reset();
+    source.snapshot = std::make_shared<const RecordedSnapshot>();
+    return source;
+  }
+
+  CHECK_EQ(source.histogram->counts_len, histogram_->counts_len);
+  source.histogram->min_value = histogram_->min_value;
+  source.histogram->max_value = histogram_->max_value;
+  source.histogram->normalizing_index_offset =
+      histogram_->normalizing_index_offset;
+  source.histogram->conversion_ratio = histogram_->conversion_ratio;
+  source.histogram->total_count = histogram_->total_count;
+  std::memcpy(source.histogram->counts,
+              histogram_->counts,
+              histogram_->counts_len * sizeof(*histogram_->counts));
+  return source;
+}
+
+void Histogram::CacheRecordedSnapshot(
+    uint64_t generation, std::shared_ptr<const RecordedSnapshot> snapshot) {
+  RwLock::ScopedWriteLock lock(mutex_);
+  if (generation == mutation_generation_) {
+    recorded_snapshot_cache_ = std::move(snapshot);
+  }
+}
+
+size_t Histogram::GetCachedRecordedSnapshotMemorySize() const {
+  RwLock::ScopedReadLock lock(mutex_);
+  if (!recorded_snapshot_cache_) return 0;
+  return recorded_snapshot_cache_->buckets.capacity() * sizeof(RecordedBucket);
+}
+
 double Histogram::Cdf(int64_t value) const {
   RwLock::ScopedReadLock lock(mutex_);
   int64_t total = histogram_->total_count;
@@ -169,6 +275,7 @@ double Histogram::Subtract(const Histogram& other) {
       histogram_->counts[i] = count;
     }
     hdr_reset_internal_counters(histogram_.get());
+    InvalidateRecordedSnapshot();
     exceeds_ = (exceeds_ > other.exceeds_) ? exceeds_ - other.exceeds_ : 0;
     return static_cast<double>(dropped);
   };
@@ -273,22 +380,157 @@ static double BetaContinuedFraction(double a, double b, double x) {
   return h;
 }
 
+struct BetaParameters {
+  double a;
+  double b;
+  double log_normalization;
+  double mean;
+  double standard_deviation;
+  double skewness;
+  double excess_kurtosis;
+  double log_scale;
+  double symmetry_point;
+  bool use_asymptotic_front;
+  bool use_asymptotic_cdf;
+};
+
+static double StirlingCorrection(double x) {
+  const double inverse = 1.0 / x;
+  const double inverse_squared = inverse * inverse;
+  return inverse * (1.0 / 12.0 +
+                    inverse_squared *
+                        (-1.0 / 360.0 +
+                         inverse_squared *
+                             (1.0 / 1260.0 +
+                              inverse_squared *
+                                  (-1.0 / 1680.0 + inverse_squared / 1188.0))));
+}
+
+static BetaParameters MakeBetaParameters(double a,
+                                         double b,
+                                         bool approximate_large_shapes) {
+  const double sum = a + b;
+  const double mean = a / sum;
+  const bool use_asymptotic_front =
+      approximate_large_shapes && a >= 8.0 && b >= 8.0;
+  // The continued fraction needs progressively more iterations as both
+  // beta shapes grow. At this concentration the Edgeworth error is smaller
+  // than the continued-fraction truncation error.
+  const bool use_asymptotic_cdf =
+      approximate_large_shapes && a * b / sum >= 25'000.0;
+  const double variance = mean * (1.0 - mean) / (sum + 1.0);
+  const double skewness =
+      2.0 * (b - a) * std::sqrt(sum + 1.0) / ((sum + 2.0) * std::sqrt(a * b));
+  const double excess_kurtosis =
+      6.0 * ((a - b) * (a - b) * (sum + 1.0) - a * b * (sum + 2.0)) /
+      (a * b * (sum + 2.0) * (sum + 3.0));
+  return {
+      a,
+      b,
+      use_asymptotic_front ? 0.0
+                           : std::lgamma(sum) - std::lgamma(a) - std::lgamma(b),
+      mean,
+      std::sqrt(variance),
+      skewness,
+      excess_kurtosis,
+      use_asymptotic_front ? 0.5 * std::log(sum * mean * (1.0 - mean) /
+                                            (2.0 * std::numbers::pi)) +
+                                 StirlingCorrection(sum) -
+                                 StirlingCorrection(a) - StirlingCorrection(b)
+                           : 0.0,
+      (a + 1.0) / (sum + 2.0),
+      use_asymptotic_front,
+      use_asymptotic_cdf,
+  };
+}
+
+static BetaParameters ReflectBetaParameters(const BetaParameters& params) {
+  return {
+      params.b,
+      params.a,
+      params.log_normalization,
+      1.0 - params.mean,
+      params.standard_deviation,
+      -params.skewness,
+      params.excess_kurtosis,
+      params.log_scale,
+      1.0 - params.symmetry_point,
+      params.use_asymptotic_front,
+      params.use_asymptotic_cdf,
+  };
+}
+
+static double BetaFront(const BetaParameters& params, double x) {
+  if (x <= 0.0 || x >= 1.0) return 0.0;
+  double log_front;
+  if (params.use_asymptotic_front) {
+    log_front = params.a * std::log1p((x - params.mean) / params.mean) +
+                params.b * std::log1p((params.mean - x) / (1.0 - params.mean)) +
+                params.log_scale;
+  } else {
+    log_front = params.log_normalization + params.a * std::log(x) +
+                params.b * std::log(1.0 - x);
+  }
+  return std::exp(log_front);
+}
+
+static double AsymptoticBetaCdf(const BetaParameters& params,
+                                double x,
+                                double* front_out) {
+  const double z = (x - params.mean) / params.standard_deviation;
+  const double normal_cdf = 0.5 * std::erfc(-z / std::numbers::sqrt2);
+  if (std::fabs(z) >= 8.0) {
+    if (front_out != nullptr) *front_out = 0.0;
+    return normal_cdf;
+  }
+
+  const double z2 = z * z;
+  const double z3 = z2 * z;
+  const double normal_pdf =
+      std::exp(-0.5 * z2) / std::sqrt(2.0 * std::numbers::pi);
+  if (front_out != nullptr) {
+    const double z4 = z2 * z2;
+    const double z6 = z3 * z3;
+    const double density_correction =
+        1.0 + params.skewness / 6.0 * (z3 - 3.0 * z) +
+        params.excess_kurtosis / 24.0 * (z4 - 6.0 * z2 + 3.0) +
+        params.skewness * params.skewness / 72.0 *
+            (z6 - 15.0 * z4 + 45.0 * z2 - 15.0);
+    *front_out = x * (1.0 - x) * normal_pdf / params.standard_deviation *
+                 std::max(0.0, density_correction);
+  }
+  const double correction = params.skewness / 6.0 * (1.0 - z2) -
+                            params.excess_kurtosis / 24.0 * (z3 - 3.0 * z) -
+                            params.skewness * params.skewness / 72.0 *
+                                (z3 * z2 - 10.0 * z3 + 15.0 * z);
+  return std::clamp(normal_cdf + normal_pdf * correction, 0.0, 1.0);
+}
+
 // Regularized incomplete beta function I_x(a, b).
 // Returns the probability that a Beta(a,b) random variable is <= x.
-static double RegularizedIncompleteBeta(double a, double b, double x) {
+static double RegularizedIncompleteBeta(const BetaParameters& params,
+                                        double x,
+                                        double* front_out = nullptr) {
+  if (front_out != nullptr) *front_out = 0.0;
   if (x <= 0.0) return 0.0;
   if (x >= 1.0) return 1.0;
-
-  double ln_front = std::lgamma(a + b) - std::lgamma(a) - std::lgamma(b) +
-                    a * std::log(x) + b * std::log(1.0 - x);
-  double bt = std::exp(ln_front);
+  if (params.use_asymptotic_cdf) {
+    return AsymptoticBetaCdf(params, x, front_out);
+  }
+  const double front = BetaFront(params, x);
+  if (front_out != nullptr) *front_out = front;
 
   // Use the symmetry relation to ensure the continued fraction
   // converges in the region where it is most accurate.
-  if (x < (a + 1.0) / (a + b + 2.0)) {
-    return bt * BetaContinuedFraction(a, b, x) / a;
+  if (x < params.symmetry_point) {
+    return front * BetaContinuedFraction(params.a, params.b, x) / params.a;
   }
-  return 1.0 - bt * BetaContinuedFraction(b, a, 1.0 - x) / b;
+  return 1.0 -
+         front * BetaContinuedFraction(params.b, params.a, 1.0 - x) / params.b;
+}
+
+static double RegularizedIncompleteBeta(double a, double b, double x) {
+  return RegularizedIncompleteBeta(MakeBetaParameters(a, b, false), x);
 }
 
 // Standard normal CDF: Phi(x) = P(Z <= x).
@@ -338,6 +580,381 @@ static double BinomialCdf(int64_t k, int64_t n, double p) {
       static_cast<double>(n - k), static_cast<double>(k + 1), 1.0 - p);
 }
 
+struct QrdeResult {
+  std::vector<double> quantiles;
+  std::vector<double> densities;
+  size_t corrections = 0;
+};
+
+enum class QrdeDequantization : uint32_t {
+  kNone,
+  kHdr,
+  kAll,
+};
+
+static bool QrdeShouldDequantize(const Histogram::RecordedBucket& bucket,
+                                 QrdeDequantization dequantization) {
+  return bucket.count > 1 && (dequantization == QrdeDequantization::kAll ||
+                              (dequantization == QrdeDequantization::kHdr &&
+                               bucket.resolution > 1.0));
+}
+
+static std::pair<double, double> QrdeBucketRange(
+    const Histogram::RecordedSnapshot& snapshot,
+    size_t index,
+    QrdeDequantization dequantization) {
+  const auto& bucket = snapshot.buckets[index];
+  if (!QrdeShouldDequantize(bucket, dequantization)) {
+    return {bucket.value, bucket.value};
+  }
+
+  const double half_resolution = bucket.resolution / 2.0;
+  if (snapshot.buckets.size() == 1) {
+    return {bucket.value - half_resolution, bucket.value + half_resolution};
+  }
+  if (index == 0) {
+    return {bucket.value, bucket.value + half_resolution};
+  }
+  if (index + 1 == snapshot.buckets.size()) {
+    return {bucket.value - half_resolution, bucket.value};
+  }
+  return {bucket.value - half_resolution, bucket.value + half_resolution};
+}
+
+static bool QrdeExactSupport(const Histogram::RecordedSnapshot& snapshot,
+                             const BetaParameters& params,
+                             const BetaParameters& reflected,
+                             double count,
+                             size_t* begin_index,
+                             size_t* end_index) {
+  // Collapsing both tails changes a quantile by at most 2^-79 times the
+  // histogram range, which is less than 2^-16 over the int64 value domain.
+  constexpr double tail_mass = 0x1p-80;
+  size_t lo = 0;
+  size_t hi = snapshot.buckets.size();
+  while (lo < hi) {
+    const size_t mid = lo + (hi - lo) / 2;
+    const double rank =
+        static_cast<double>(snapshot.buckets[mid].cumulative_count) / count;
+    const double cdf = RegularizedIncompleteBeta(params, rank);
+    if (!std::isfinite(cdf)) return false;
+    if (cdf < tail_mass) {
+      lo = mid + 1;
+    } else {
+      hi = mid;
+    }
+  }
+  if (lo == snapshot.buckets.size()) return false;
+  const size_t first = lo;
+
+  lo = first;
+  hi = snapshot.buckets.size();
+  while (lo < hi) {
+    const size_t mid = lo + (hi - lo) / 2;
+    const double remaining_rank =
+        static_cast<double>(snapshot.total_count -
+                            snapshot.buckets[mid].cumulative_count) /
+        count;
+    const double survival =
+        RegularizedIncompleteBeta(reflected, remaining_rank);
+    if (!std::isfinite(survival)) return false;
+    if (survival > tail_mass) {
+      lo = mid + 1;
+    } else {
+      hi = mid;
+    }
+  }
+  if (lo == snapshot.buckets.size() || lo < first) return false;
+  *begin_index = first;
+  *end_index = lo + 1;
+  return true;
+}
+
+static double QrdeQuantile(const Histogram::RecordedSnapshot& snapshot,
+                           double p,
+                           QrdeDequantization dequantization) {
+  const double count = static_cast<double>(snapshot.total_count);
+  const double a = (count + 1.0) * p;
+  const double b = (count + 1.0) * (1.0 - p);
+  const BetaParameters mass_params = MakeBetaParameters(a, b, true);
+  const BetaParameters reflected_params = ReflectBetaParameters(mass_params);
+
+  size_t begin_index = 0;
+  size_t end_index = snapshot.buckets.size();
+  if (mass_params.use_asymptotic_cdf) {
+    // Outside this interval the normal-tail mass is below 2e-33.
+    constexpr double support_deviations = 12.0;
+    const double lower_rank = std::max(
+        0.0,
+        mass_params.mean - support_deviations * mass_params.standard_deviation);
+    const double upper_rank = std::min(
+        1.0,
+        mass_params.mean + support_deviations * mass_params.standard_deviation);
+    const auto first = std::lower_bound(
+        snapshot.buckets.begin(),
+        snapshot.buckets.end(),
+        lower_rank * count,
+        [](const Histogram::RecordedBucket& bucket, double cumulative) {
+          return static_cast<double>(bucket.cumulative_count) < cumulative;
+        });
+    const auto last = std::lower_bound(
+        first,
+        snapshot.buckets.end(),
+        upper_rank * count,
+        [](const Histogram::RecordedBucket& bucket, double cumulative) {
+          return static_cast<double>(bucket.cumulative_count) < cumulative;
+        });
+    begin_index = first - snapshot.buckets.begin();
+    if (last != snapshot.buckets.end()) {
+      end_index = last - snapshot.buckets.begin() + 1;
+    }
+  } else if (snapshot.buckets.size() >= 512) {
+    QrdeExactSupport(snapshot,
+                     mass_params,
+                     reflected_params,
+                     count,
+                     &begin_index,
+                     &end_index);
+  }
+
+  const bool dequantize = dequantization != QrdeDequantization::kNone;
+  const auto first_range = QrdeBucketRange(snapshot, 0, dequantization);
+  const auto last_range =
+      QrdeBucketRange(snapshot, snapshot.buckets.size() - 1, dequantization);
+  const double lower_endpoint = first_range.first;
+  const double upper_endpoint = last_range.second;
+
+  double previous_cdf = 0.0;
+  double previous_survival = 0.0;
+  double previous_front = 0.0;
+  int64_t previous_count =
+      begin_index == 0 ? 0 : snapshot.buckets[begin_index - 1].cumulative_count;
+  if (previous_count != 0) {
+    const double previous_rank = static_cast<double>(previous_count) / count;
+    previous_cdf = RegularizedIncompleteBeta(
+        mass_params, previous_rank, dequantize ? &previous_front : nullptr);
+    if (!std::isfinite(previous_cdf) ||
+        (dequantize && !std::isfinite(previous_front))) {
+      begin_index = 0;
+      end_index = snapshot.buckets.size();
+      previous_count = 0;
+      previous_cdf = 0.0;
+      previous_front = 0.0;
+    } else {
+      previous_cdf = std::clamp(previous_cdf, 0.0, 1.0);
+    }
+  }
+
+  double quantile = lower_endpoint * previous_cdf;
+  bool using_survival = false;
+
+  for (size_t i = begin_index; i < end_index; i++) {
+    const auto& bucket = snapshot.buckets[i];
+    const double u0 = static_cast<double>(previous_count) / count;
+    const double u1 = static_cast<double>(bucket.cumulative_count) / count;
+    double front = 0.0;
+    double mass;
+    // Use the reflected CDF above the symmetry point so small upper-tail
+    // interval weights are not lost to subtraction from one.
+    if (u1 <= mass_params.symmetry_point) {
+      const double cdf =
+          std::clamp(RegularizedIncompleteBeta(
+                         mass_params, u1, dequantize ? &front : nullptr),
+                     previous_cdf,
+                     1.0);
+      mass = cdf - previous_cdf;
+      previous_cdf = cdf;
+    } else {
+      const double remaining_rank =
+          static_cast<double>(snapshot.total_count - bucket.cumulative_count) /
+          count;
+      double survival = RegularizedIncompleteBeta(
+          reflected_params, remaining_rank, dequantize ? &front : nullptr);
+      if (using_survival) {
+        survival = std::clamp(survival, 0.0, previous_survival);
+        mass = previous_survival - survival;
+      } else {
+        survival = std::clamp(survival, 0.0, 1.0 - previous_cdf);
+        mass = 1.0 - previous_cdf - survival;
+        using_survival = true;
+      }
+      previous_survival = survival;
+    }
+
+    if (!QrdeShouldDequantize(bucket, dequantization)) {
+      quantile += bucket.value * mass;
+    } else {
+      // I_x(a + 1, b) = I_x(a, b) - front / a, so the first moment
+      // within this rank interval does not require a second beta CDF.
+      const double sum = a + b;
+      const double rank_width = static_cast<double>(bucket.count) / count;
+      const double local_moment =
+          std::clamp((a / sum - u0) * mass - (front - previous_front) / sum,
+                     0.0,
+                     rank_width * mass);
+      const auto [lower, upper] = QrdeBucketRange(snapshot, i, dequantization);
+      quantile += lower * mass + (upper - lower) * local_moment / rank_width;
+    }
+
+    previous_front = front;
+    previous_count = bucket.cumulative_count;
+  }
+
+  const double remaining_mass =
+      using_survival ? previous_survival : 1.0 - previous_cdf;
+  quantile += upper_endpoint * std::clamp(remaining_mass, 0.0, 1.0);
+  return quantile;
+}
+
+static QrdeResult CalculateQrde(const Histogram::RecordedSnapshot& snapshot,
+                                const std::vector<double>& probabilities,
+                                QrdeDequantization dequantization) {
+  QrdeResult result;
+  if (snapshot.buckets.empty()) return result;
+
+  result.quantiles.resize(probabilities.size());
+  result.densities.resize(probabilities.size() - 1);
+
+  const auto first_range = QrdeBucketRange(snapshot, 0, dequantization);
+  const auto last_range =
+      QrdeBucketRange(snapshot, snapshot.buckets.size() - 1, dequantization);
+  result.quantiles.front() = first_range.first;
+  result.quantiles.back() = last_range.second;
+
+  for (size_t i = 1; i + 1 < probabilities.size(); i++) {
+    double value = QrdeQuantile(snapshot, probabilities[i], dequantization);
+    value =
+        std::clamp(value, result.quantiles.front(), result.quantiles.back());
+    if (value < result.quantiles[i - 1]) {
+      value = result.quantiles[i - 1];
+      result.corrections++;
+    }
+    result.quantiles[i] = value;
+  }
+
+  for (size_t i = 0; i < result.densities.size(); i++) {
+    const double width = result.quantiles[i + 1] - result.quantiles[i];
+    const double probability_mass = probabilities[i + 1] - probabilities[i];
+    result.densities[i] = width == 0.0 ? std::numeric_limits<double>::infinity()
+                                       : probability_mass / width;
+  }
+
+  return result;
+}
+
+static Local<Float64Array> ToFloat64Array(Isolate* isolate,
+                                          const std::vector<double>& values) {
+  auto store =
+      ArrayBuffer::NewBackingStore(isolate, values.size() * sizeof(double));
+  if (!values.empty()) {
+    memcpy(store->Data(), values.data(), values.size() * sizeof(double));
+  }
+  Local<ArrayBuffer> buffer = ArrayBuffer::New(isolate, std::move(store));
+  return Float64Array::New(buffer, 0, values.size());
+}
+
+static const char* QrdeDequantizationName(QrdeDequantization dequantization) {
+  switch (dequantization) {
+    case QrdeDequantization::kNone:
+      return "none";
+    case QrdeDequantization::kHdr:
+      return "hdr";
+    case QrdeDequantization::kAll:
+      return "all";
+  }
+  UNREACHABLE();
+}
+
+class QrdeJob final : public ThreadPoolWork {
+ public:
+  QrdeJob(Environment* env,
+          Local<Promise::Resolver> resolver,
+          std::shared_ptr<Histogram> histogram,
+          Histogram::RecordedSnapshotSource snapshot_source,
+          std::vector<double> probabilities,
+          QrdeDequantization dequantization)
+      : ThreadPoolWork(env, "histogram.qrde"),
+        histogram_(std::move(histogram)),
+        snapshot_source_(std::move(snapshot_source)),
+        probabilities_(std::move(probabilities)),
+        dequantization_(dequantization) {
+    resolver_.Reset(env->isolate(), resolver);
+  }
+
+  void DoThreadPoolWork() override {
+    if (snapshot_source_.snapshot) {
+      snapshot_ = std::move(snapshot_source_.snapshot);
+    } else {
+      snapshot_ = std::make_shared<const Histogram::RecordedSnapshot>(
+          Histogram::BuildRecordedSnapshot(snapshot_source_.histogram.get()));
+      snapshot_source_.histogram.reset();
+    }
+    if (histogram_) {
+      histogram_->CacheRecordedSnapshot(snapshot_source_.generation, snapshot_);
+      histogram_.reset();
+    }
+    result_ = CalculateQrde(*snapshot_, probabilities_, dequantization_);
+  }
+
+  void AfterThreadPoolWork(int status) override {
+    std::unique_ptr<QrdeJob> self(this);
+    Environment* env = ThreadPoolWork::env();
+    CHECK(status == 0 || status == UV_ECANCELED);
+    if (!env->can_call_into_js()) {
+      resolver_.Reset();
+      return;
+    }
+
+    Isolate* isolate = env->isolate();
+    HandleScope handle_scope(isolate);
+    Context::Scope context_scope(env->context());
+    InternalCallbackScope callback_scope(
+        env, Object::New(isolate), {0, 0}, InternalCallbackScope::kNoFlags);
+    Local<Promise::Resolver> resolver =
+        Local<Promise::Resolver>::New(isolate, resolver_);
+
+    if (status == UV_ECANCELED) {
+      USE(resolver->Reject(
+          env->context(),
+          Exception::Error(OneByteString(isolate, "QRDE job was canceled"))));
+      resolver_.Reset();
+      return;
+    }
+
+    Local<Name> names[] = {
+        FIXED_ONE_BYTE_STRING(isolate, "probabilities"),
+        FIXED_ONE_BYTE_STRING(isolate, "quantiles"),
+        FIXED_ONE_BYTE_STRING(isolate, "densities"),
+        FIXED_ONE_BYTE_STRING(isolate, "count"),
+        FIXED_ONE_BYTE_STRING(isolate, "bucketCount"),
+        FIXED_ONE_BYTE_STRING(isolate, "corrections"),
+        FIXED_ONE_BYTE_STRING(isolate, "dequantize"),
+    };
+    Local<Value> values[] = {
+        ToFloat64Array(isolate, probabilities_),
+        ToFloat64Array(isolate, result_.quantiles),
+        ToFloat64Array(isolate, result_.densities),
+        BigInt::New(isolate, snapshot_->total_count),
+        Number::New(isolate, static_cast<double>(snapshot_->buckets.size())),
+        Number::New(isolate, static_cast<double>(result_.corrections)),
+        OneByteString(isolate, QrdeDequantizationName(dequantization_)),
+    };
+    Local<Object> value =
+        Object::New(isolate, Null(isolate), names, values, arraysize(names));
+    USE(resolver->Resolve(env->context(), value));
+    resolver_.Reset();
+  }
+
+ private:
+  std::shared_ptr<Histogram> histogram_;
+  Histogram::RecordedSnapshotSource snapshot_source_;
+  std::shared_ptr<const Histogram::RecordedSnapshot> snapshot_;
+  const std::vector<double> probabilities_;
+  const QrdeDequantization dequantization_;
+  QrdeResult result_;
+  Global<Promise::Resolver> resolver_;
+};
+
 // -----------------------------------------------------------------------
 // Minimal CBOR encoder/decoder (RFC 8949) -- just enough types for
 // histogram export/import: unsigned int, float64, array, and map.
@@ -1176,6 +1793,7 @@ void HistogramImpl::AddMethods(Isolate* isolate, Local<FunctionTemplate> tmpl) {
   SetProtoMethodNoSideEffect(isolate, tmpl, "cohensD", GetCohensD);
   SetProtoMethodNoSideEffect(isolate, tmpl, "cliffsD", GetCliffsD);
   SetProtoMethodNoSideEffect(isolate, tmpl, "percentileCI", GetPercentileCI);
+  SetProtoMethod(isolate, tmpl, "qrde", GetQrde);
   SetFastMethodNoSideEffect(
       isolate, instance, "ewmaMean", GetEwmaMean, &fast_get_ewma_mean_);
   SetFastMethodNoSideEffect(
@@ -1230,6 +1848,7 @@ void HistogramImpl::RegisterExternalReferences(
   registry->Register(GetCohensD);
   registry->Register(GetCliffsD);
   registry->Register(GetPercentileCI);
+  registry->Register(GetQrde);
   registry->Register(GetEwmaMean);
   registry->Register(GetEwmaStddev);
   registry->Register(GetEwmaErrorRate);
@@ -1244,12 +1863,10 @@ void HistogramImpl::RegisterExternalReferences(
   is_registered = true;
 }
 
-HistogramBase::HistogramBase(
-    Environment* env,
-    Local<Object> wrap,
-    const Histogram::Options& options)
-    : BaseObject(env, wrap),
-      HistogramImpl(options) {
+HistogramBase::HistogramBase(Environment* env,
+                             Local<Object> wrap,
+                             const Histogram::Options& options)
+    : BaseObject(env, wrap), HistogramImpl(options) {
   MakeWeak();
   wrap->SetAlignedPointerInInternalField(
       HistogramImpl::InternalFields::kImplField,
@@ -1257,12 +1874,10 @@ HistogramBase::HistogramBase(
       EmbedderDataTag::kDefault);
 }
 
-HistogramBase::HistogramBase(
-    Environment* env,
-    Local<Object> wrap,
-    std::shared_ptr<Histogram> histogram)
-    : BaseObject(env, wrap),
-      HistogramImpl(std::move(histogram)) {
+HistogramBase::HistogramBase(Environment* env,
+                             Local<Object> wrap,
+                             std::shared_ptr<Histogram> histogram)
+    : BaseObject(env, wrap), HistogramImpl(std::move(histogram)) {
   MakeWeak();
   wrap->SetAlignedPointerInInternalField(
       HistogramImpl::InternalFields::kImplField,
@@ -1292,8 +1907,8 @@ void HistogramBase::Record(const FunctionCallbackInfo<Value>& args) {
   CHECK_IMPLIES(!args[0]->IsNumber(), args[0]->IsBigInt());
   bool lossless = true;
   int64_t value = args[0]->IsBigInt()
-      ? args[0].As<BigInt>()->Int64Value(&lossless)
-      : static_cast<int64_t>(args[0].As<Number>()->Value());
+                      ? args[0].As<BigInt>()->Int64Value(&lossless)
+                      : static_cast<int64_t>(args[0].As<Number>()->Value());
   if (!lossless || value < 1)
     return THROW_ERR_OUT_OF_RANGE(env, "value is out of range");
   HistogramBase* histogram;
@@ -1356,8 +1971,7 @@ void HistogramBase::RecordCorrected(const FunctionCallbackInfo<Value>& args) {
 }
 
 BaseObjectPtr<HistogramBase> HistogramBase::Create(
-    Environment* env,
-    const Histogram::Options& options) {
+    Environment* env, const Histogram::Options& options) {
   Local<Object> obj;
   if (!GetConstructorTemplate(env->isolate_data())
            ->InstanceTemplate()
@@ -1370,8 +1984,7 @@ BaseObjectPtr<HistogramBase> HistogramBase::Create(
 }
 
 BaseObjectPtr<HistogramBase> HistogramBase::Create(
-    Environment* env,
-    std::shared_ptr<Histogram> histogram) {
+    Environment* env, std::shared_ptr<Histogram> histogram) {
   Local<Object> obj;
   if (!GetConstructorTemplate(env->isolate_data())
            ->InstanceTemplate()
@@ -1538,8 +2151,9 @@ BaseObjectPtr<IntervalHistogram> IntervalHistogram::Create(
     AsyncWrap::ProviderType type) {
   Local<Object> obj;
   if (!GetConstructorTemplate(env)
-          ->InstanceTemplate()
-          ->NewInstance(env->context()).ToLocal(&obj)) {
+           ->InstanceTemplate()
+           ->NewInstance(env->context())
+           .ToLocal(&obj)) {
     return nullptr;
   }
 
@@ -1563,8 +2177,7 @@ void IntervalHistogram::MemoryInfo(MemoryTracker* tracker) const {
 void IntervalHistogram::OnStart(StartFlags flags) {
   if (enabled_ || IsHandleClosing()) return;
   enabled_ = true;
-  if (flags == StartFlags::RESET)
-    histogram()->Reset();
+  if (flags == StartFlags::RESET) histogram()->Reset();
   uv_timer_start(&timer_, TimerCB, interval_, interval_);
   uv_unref(reinterpret_cast<uv_handle_t*>(&timer_));
 }
@@ -2011,6 +2624,47 @@ void HistogramImpl::GetPercentileCI(const FunctionCallbackInfo<Value>& args) {
   args.GetReturnValue().Set(arr);
 }
 
+void HistogramImpl::GetQrde(const FunctionCallbackInfo<Value>& args) {
+  Environment* env = Environment::GetCurrent(args);
+  HistogramImpl* histogram = HistogramImpl::FromJSObject(args.This());
+  CHECK(args[0]->IsFloat64Array());
+  CHECK(args[1]->IsUint32());
+  CHECK(args[2]->IsBoolean());
+  Local<Float64Array> input = args[0].As<Float64Array>();
+  CHECK_GE(input->Length(), 2);
+  CHECK_LE(input->Length(), 1001);
+  auto backing = input->Buffer()->GetBackingStore();
+  const double* data = reinterpret_cast<const double*>(
+      static_cast<const char*>(backing->Data()) + input->ByteOffset());
+  std::vector<double> probabilities(data, data + input->Length());
+  const uint32_t dequantization = args[1].As<Uint32>()->Value();
+  CHECK_LE(dequantization, static_cast<uint32_t>(QrdeDequantization::kAll));
+  const bool cache_snapshot = args[2]->IsTrue();
+
+  Local<Promise::Resolver> resolver;
+  if (!Promise::Resolver::New(env->context()).ToLocal(&resolver)) return;
+
+  auto histogram_ptr = histogram->histogram();
+  auto snapshot_source =
+      histogram_ptr->GetRecordedSnapshotSource(cache_snapshot);
+  if (!snapshot_source.histogram && !snapshot_source.snapshot) {
+    THROW_ERR_MEMORY_ALLOCATION_FAILED(env);
+    return;
+  }
+  std::shared_ptr<Histogram> cache_target;
+  if (cache_snapshot && !snapshot_source.cache_hit) {
+    cache_target = histogram_ptr;
+  }
+  auto* job = new QrdeJob(env,
+                          resolver,
+                          std::move(cache_target),
+                          std::move(snapshot_source),
+                          std::move(probabilities),
+                          static_cast<QrdeDequantization>(dequantization));
+  args.GetReturnValue().Set(resolver->GetPromise());
+  job->ScheduleWork();
+}
+
 void HistogramImpl::GetEwmaMean(const FunctionCallbackInfo<Value>& args) {
   HistogramImpl* histogram = HistogramImpl::FromJSObject(args.This());
   args.GetReturnValue().Set((*histogram)->EwmaMean());
@@ -2049,10 +2703,10 @@ void HistogramImpl::DoExport(const FunctionCallbackInfo<Value>& args) {
   HistogramImpl* histogram = HistogramImpl::FromJSObject(args.This());
   std::vector<uint8_t> data = (*histogram)->Export();
 
-  auto store = v8::ArrayBuffer::NewBackingStore(env->isolate(), data.size());
+  auto store = ArrayBuffer::NewBackingStore(env->isolate(), data.size());
   memcpy(store->Data(), data.data(), data.size());
-  auto buf = v8::ArrayBuffer::New(env->isolate(), std::move(store));
-  auto arr = v8::Uint8Array::New(buf, 0, data.size());
+  auto buf = ArrayBuffer::New(env->isolate(), std::move(store));
+  auto arr = Uint8Array::New(buf, 0, data.size());
   args.GetReturnValue().Set(arr);
 }
 
@@ -2062,7 +2716,7 @@ void HistogramImpl::DoImport(const FunctionCallbackInfo<Value>& args) {
     THROW_ERR_INVALID_ARG_TYPE(env, "data must be a Uint8Array");
     return;
   }
-  Local<v8::Uint8Array> input = args[0].As<v8::Uint8Array>();
+  Local<Uint8Array> input = args[0].As<Uint8Array>();
   auto backing = input->Buffer()->GetBackingStore();
   const uint8_t* data =
       static_cast<const uint8_t*>(backing->Data()) + input->ByteOffset();
@@ -2163,8 +2817,8 @@ std::unique_ptr<worker::TransferData> IterationHistogram::CloneForMessaging()
   return std::make_unique<HistogramBase::HistogramTransferData>(histogram());
 }
 
-std::unique_ptr<worker::TransferData>
-IntervalHistogram::CloneForMessaging() const {
+std::unique_ptr<worker::TransferData> IntervalHistogram::CloneForMessaging()
+    const {
   return std::make_unique<HistogramBase::HistogramTransferData>(histogram());
 }
 

src/histogram.h

@@ -37,8 +37,27 @@ class Histogram : public MemoryRetainer {
                             // exceeding this threshold.
   };
 
+  struct RecordedBucket {
+    double value;
+    double resolution;
+    int64_t count;
+    int64_t cumulative_count;
+  };
+
+  struct RecordedSnapshot {
+    std::vector<RecordedBucket> buckets;
+    int64_t total_count = 0;
+  };
+
   using HistogramPointer = DeleteFnPtr<hdr_histogram, hdr_close>;
 
+  struct RecordedSnapshotSource {
+    HistogramPointer histogram;
+    std::shared_ptr<const RecordedSnapshot> snapshot;
+    uint64_t generation = 0;
+    bool cache_hit = false;
+  };
+
   // Factory method that returns nullptr on hdr_init failure.
   static std::shared_ptr<Histogram> Create(const Options& options);
 
@@ -129,6 +148,10 @@ class Histogram : public MemoryRetainer {
   void LogBuckets(int64_t first_bucket, double log_base, Iterator&& fn) const;
 
   bool IsCompatible(const Histogram& other) const;
+  RecordedSnapshotSource GetRecordedSnapshotSource(bool use_cache) const;
+  static RecordedSnapshot BuildRecordedSnapshot(const hdr_histogram* histogram);
+  void CacheRecordedSnapshot(uint64_t generation,
+                             std::shared_ptr<const RecordedSnapshot> snapshot);
 
   void MemoryInfo(MemoryTracker* tracker) const override;
   SET_MEMORY_INFO_NAME(Histogram)
@@ -136,6 +159,8 @@ class Histogram : public MemoryRetainer {
 
  private:
   inline void UpdateEwma(double value);
+  inline void InvalidateRecordedSnapshot();
+  size_t GetCachedRecordedSnapshotMemorySize() const;
 
   HistogramPointer histogram_;
   uint64_t prev_ = 0;
@@ -151,6 +176,9 @@ class Histogram : public MemoryRetainer {
   int64_t threshold_ = 0;
   double ewma_error_rate_ = 0;
 
+  uint64_t mutation_generation_ = 0;
+  std::shared_ptr<const RecordedSnapshot> recorded_snapshot_cache_;
+
   RwLock mutex_;
 };
 
@@ -203,6 +231,7 @@ class HistogramImpl {
   static void GetCohensD(const v8::FunctionCallbackInfo<v8::Value>& args);
   static void GetCliffsD(const v8::FunctionCallbackInfo<v8::Value>& args);
   static void GetPercentileCI(const v8::FunctionCallbackInfo<v8::Value>& args);
+  static void GetQrde(const v8::FunctionCallbackInfo<v8::Value>& args);
   static void GetEwmaMean(const v8::FunctionCallbackInfo<v8::Value>& args);
   static void GetEwmaStddev(const v8::FunctionCallbackInfo<v8::Value>& args);
   static void GetEwmaErrorRate(const v8::FunctionCallbackInfo<v8::Value>& args);

test/fixtures/qrde-r-oracle.R

@@ -0,0 +1,133 @@
+# Regenerate qrde-r-oracle.json with:
+# Rscript test/fixtures/qrde-r-oracle.R > test/fixtures/qrde-r-oracle.json
+
+beta_front <- function(x, a, b) {
+  result <- numeric(length(x))
+  interior <- x > 0 & x < 1
+  result[interior] <- x[interior] * (1 - x[interior]) *
+    dbeta(x[interior], a, b)
+  result
+}
+
+bucket_ranges <- function(values, resolutions, counts, mode) {
+  lower <- upper <- values
+  spread <- counts > 1 &
+    (mode == "all" | (mode == "hdr" & resolutions > 1))
+  half <- resolutions / 2
+  lower[spread] <- values[spread] - half[spread]
+  upper[spread] <- values[spread] + half[spread]
+  if (length(values) > 1) {
+    if (spread[1]) lower[1] <- values[1]
+    if (spread[length(values)]) upper[length(values)] <- values[length(values)]
+  }
+  list(lower = lower, upper = upper, spread = spread)
+}
+
+hd_quantile <- function(values, resolutions, counts, p, mode) {
+  n <- sum(counts)
+  a <- (n + 1) * p
+  b <- (n + 1) * (1 - p)
+  u1 <- cumsum(counts) / n
+  u0 <- c(0, head(u1, -1))
+  threshold <- (a + 1) / (a + b + 2)
+  lower <- u1 <= threshold
+  upper <- u0 >= threshold
+  crossing <- !(lower | upper)
+  mass <- numeric(length(values))
+  mass[lower] <- pbeta(u1[lower], a, b) - pbeta(u0[lower], a, b)
+  mass[upper] <- pbeta(1 - u0[upper], b, a) -
+    pbeta(1 - u1[upper], b, a)
+  mass[crossing] <- 1 - pbeta(u0[crossing], a, b) -
+    pbeta(1 - u1[crossing], b, a)
+
+  ranges <- bucket_ranges(values, resolutions, counts, mode)
+  if (!any(ranges$spread)) return(sum(values * mass))
+
+  front0 <- beta_front(u0, a, b)
+  front1 <- beta_front(u1, a, b)
+  local <- (p - u0) * mass - (front1 - front0) / (a + b)
+  rank_width <- counts / n
+  local <- pmax(0, pmin(rank_width * mass, local))
+  sum(ifelse(ranges$spread,
+             ranges$lower * mass +
+               (ranges$upper - ranges$lower) * local / rank_width,
+             values * mass))
+}
+
+make_case <- function(name, values, resolutions, counts, bins, mode, indices) {
+  ranges <- bucket_ranges(values, resolutions, counts, mode)
+  expected <- vapply(indices, function(index) {
+    if (index == 0) return(ranges$lower[1])
+    if (index == bins) return(tail(ranges$upper, 1))
+    hd_quantile(values, resolutions, counts, index / bins, mode)
+  }, numeric(1))
+  list(name = name, indices = indices, expected = expected)
+}
+
+make_probability_case <- function(name, values, resolutions, counts,
+                                  probabilities, mode) {
+  ranges <- bucket_ranges(values, resolutions, counts, mode)
+  expected <- vapply(probabilities, function(probability) {
+    if (probability == 0) return(ranges$lower[1])
+    if (probability == 1) return(tail(ranges$upper, 1))
+    hd_quantile(values, resolutions, counts, probability, mode)
+  }, numeric(1))
+  list(name = name, probabilities = probabilities, expected = expected)
+}
+
+cases <- list(
+  make_case("small-none", c(1, 2, 4, 16, 100), rep(1, 5),
+            c(1, 2, 5, 3, 1), 10, "none", 0:10),
+  make_case("small-all", c(1, 2, 4, 16, 100), rep(1, 5),
+            c(1, 2, 5, 3, 1), 10, "all", 0:10),
+  make_case("exact-support", 1:1024, rep(1, 1024), rep(1, 1024),
+            100, "none", c(1, 2, 10, 25, 50, 75, 90, 98, 99)),
+  make_probability_case("tail-focused", 1:1024, rep(1, 1024), rep(1, 1024),
+                        c(0, 0.5, 0.9, 0.99, 0.999, 0.9999, 1), "none"),
+  make_case("approximation-threshold", c(1, 131071), c(1, 1),
+            c(26239, 973761), 1000, "none",
+            c(1, 10, 24, 25, 26, 27, 50, 500, 950, 973, 974, 975, 976,
+              990, 999)),
+  make_case("multimodal-none", c(1, 10, 100, 1000), rep(1, 4),
+            c(900000, 90000, 9000, 1000), 100, "none",
+            c(1, 10, 25, 50, 75, 89, 90, 91, 98, 99)),
+  make_case("multimodal-all", c(1, 10, 100, 1000), rep(1, 4),
+            c(900000, 90000, 9000, 1000), 100, "all",
+            c(1, 10, 25, 50, 75, 89, 90, 91, 98, 99)),
+  make_case("wide-none",
+            c(1049088, 1100048498688, 4505798650626048, 9005000231485440),
+            c(1024, 1073741824, 4398046511104, 4398046511104),
+            c(700000, 200000, 99999, 1), 1000, "none",
+            c(1, 100, 500, 699, 700, 701, 899, 900, 901, 990, 999)),
+  make_case("wide-hdr",
+            c(1049088, 1100048498688, 4505798650626048, 9005000231485440),
+            c(1024, 1073741824, 4398046511104, 4398046511104),
+            c(700000, 200000, 99999, 1), 1000, "hdr",
+            c(1, 100, 500, 699, 700, 701, 899, 900, 901, 990, 999)),
+  make_case("huge-count", c(1, 1000), c(1, 1),
+            c(2^50, 2 * 2^50), 3, "none", c(1, 2))
+)
+
+number_array <- function(values) {
+  paste0("[", paste(sprintf("%.17g", values), collapse = ", "), "]")
+}
+
+cat("{\n")
+cat(sprintf("  \"rVersion\": \"%s\",\n", getRversion()))
+cat("  \"implementation\": \"stats::pbeta/stats::dbeta\",\n")
+cat("  \"generator\": \"test/fixtures/qrde-r-oracle.R\",\n")
+cat("  \"cases\": {\n")
+for (i in seq_along(cases)) {
+  case <- cases[[i]]
+  cat(sprintf("    \"%s\": {\n", case$name))
+  if (is.null(case$indices)) {
+    cat(sprintf("      \"probabilities\": %s,\n",
+                number_array(case$probabilities)))
+  } else {
+    cat(sprintf("      \"indices\": %s,\n", number_array(case$indices)))
+  }
+  cat(sprintf("      \"quantiles\": %s\n", number_array(case$expected)))
+  cat(sprintf("    }%s\n", if (i == length(cases)) "" else ","))
+}
+cat("  }\n")
+cat("}\n")

test/fixtures/qrde-r-oracle.json

@@ -0,0 +1,47 @@
+{
+  "rVersion": "4.3.3",
+  "implementation": "stats::pbeta/stats::dbeta",
+  "generator": "test/fixtures/qrde-r-oracle.R",
+  "cases": {
+    "small-none": {
+      "indices": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
+      "quantiles": [1, 1.5965008927561191, 2.4535177038411491, 3.2717431130126258, 4.0354421521907682, 5.2952901873257829, 7.993414878390019, 12.95756351062192, 25.2075453528239, 59.34721654185303, 100]
+    },
+    "small-all": {
+      "indices": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
+      "quantiles": [1, 1.5164508130154255, 2.3591698112327482, 3.1776385959166684, 3.9790848992171997, 5.3180579814394306, 8.051827450668231, 13.007946863612533, 25.296537882899184, 59.450958793096312, 100]
+    },
+    "exact-support": {
+      "indices": [1, 2, 10, 25, 50, 75, 90, 98, 99],
+      "quantiles": [10.740000001653616, 20.979999999999997, 102.90000000000001, 256.5, 512.5, 768.5, 922.09999999999991, 1004.02, 1014.2599999983464]
+    },
+    "tail-focused": {
+      "probabilities": [0, 0.5, 0.90000000000000002, 0.98999999999999999, 0.999, 0.99990000000000001, 1],
+      "quantiles": [1, 512.5, 922.09999999999991, 1014.2599999983464, 1023.3985498813823, 1023.9671500705408, 1024]
+    },
+    "approximation-threshold": {
+      "indices": [1, 10, 24, 25, 26, 27, 50, 500, 950, 973, 974, 975, 976, 990, 999],
+      "quantiles": [1, 1, 1, 1.0000000003635539, 8767.6505315696159, 131070.85769690933, 131071, 131071, 131071, 131071, 131071, 131071, 131071, 131071, 131071]
+    },
+    "multimodal-none": {
+      "indices": [1, 10, 25, 50, 75, 89, 90, 91, 98, 99],
+      "quantiles": [1, 1, 1, 1, 1, 1, 5.5031915372522544, 10, 10, 55.117879933163145]
+    },
+    "multimodal-all": {
+      "indices": [1, 10, 25, 50, 75, 89, 90, 91, 98, 99],
+      "quantiles": [1.0055555555555555, 1.0555555555555556, 1.1388888888888888, 1.2777777777777777, 1.4166666666666665, 1.4944444444444445, 5.5041002375271573, 9.6111111111111107, 10.388888888888889, 55.120539546632884]
+    },
+    "wide-none": {
+      "indices": [1, 100, 500, 699, 700, 701, 899, 900, 901, 990, 999],
+      "quantiles": [1049088, 1049088, 1049088, 16044533331.731945, 550152461987.60876, 1084084847693.1234, 3073520372530.2637, 2255046784384464.5, 4503906146402988, 4505798650626048, 4505798650626048]
+    },
+    "wide-hdr": {
+      "indices": [1, 100, 500, 699, 700, 701, 899, 900, 901, 990, 999],
+      "quantiles": [1049088.7314285715, 1049161.142857143, 1049453.7142857143, 16036716594.329391, 549884945967.66187, 1083561149070.3787, 3073089774188.3379, 2253952024297723, 4501752029575391.5, 4507557908813304, 4507953736957585]
+    },
+    "huge-count": {
+      "indices": [1, 2],
+      "quantiles": [500.49999838367461, 1000]
+    }
+  }
+}

test/parallel/test-perf-hooks-histogram-qrde-oracle.js

@@ -0,0 +1,157 @@
+'use strict';
+
+const common = require('../common');
+const fixtures = require('../common/fixtures');
+const assert = require('assert');
+const { createHistogram } = require('perf_hooks');
+
+const oracle = JSON.parse(fixtures.readSync('qrde-r-oracle.json', 'utf8'));
+
+function buildHistogram(entries, options = {}) {
+  const histogram = createHistogram(options);
+  for (const [value, count] of entries) {
+    if (count === 1) {
+      histogram.record(value);
+      continue;
+    }
+
+    const block = createHistogram(options);
+    block.record(value);
+    let remaining = count;
+    while (remaining > 0) {
+      if (remaining % 2 === 1) histogram.add(block);
+      remaining = Math.floor(remaining / 2);
+      if (remaining > 0) block.add(block);
+    }
+  }
+  return histogram;
+}
+
+function assertMatchesOracle(definition, result) {
+  const expected = oracle.cases[definition.name];
+  const span = result.quantiles.at(-1) - result.quantiles[0];
+  const tolerance = Math.max(1e-10, span * definition.relativeTolerance);
+  assert.strictEqual(result.bucketCount, definition.bucketCount);
+  const indices = expected.indices ?? expected.probabilities.map((_, i) => i);
+  if (expected.probabilities !== undefined) {
+    assert.deepStrictEqual(result.probabilities,
+                           new Float64Array(expected.probabilities));
+  }
+  for (let i = 0; i < indices.length; i++) {
+    const index = indices[i];
+    const difference = Math.abs(result.quantiles[index] -
+                                expected.quantiles[i]);
+    assert.ok(difference <= tolerance,
+              `${definition.name} p${result.probabilities[index]}: ` +
+              `${result.quantiles[index]} != ${expected.quantiles[i]} ` +
+              `(difference ${difference}, tolerance ${tolerance})`);
+  }
+}
+
+const exactTolerance = 2e-12;
+const asymptoticTolerance = 2e-8;
+const cases = [
+  {
+    name: 'small-none',
+    bins: 10,
+    dequantize: 'none',
+    entries: [[1, 1], [2, 2], [4, 5], [16, 3], [100, 1]],
+    bucketCount: 5,
+    relativeTolerance: exactTolerance,
+  },
+  {
+    name: 'small-all',
+    bins: 10,
+    dequantize: 'all',
+    entries: [[1, 1], [2, 2], [4, 5], [16, 3], [100, 1]],
+    bucketCount: 5,
+    relativeTolerance: exactTolerance,
+  },
+  {
+    name: 'exact-support',
+    bins: 100,
+    dequantize: 'none',
+    entries: Array.from({ length: 1024 }, (_, index) => [index + 1, 1]),
+    bucketCount: 1024,
+    relativeTolerance: exactTolerance,
+  },
+  {
+    name: 'tail-focused',
+    probabilities: [0, 0.5, 0.9, 0.99, 0.999, 0.9999, 1],
+    dequantize: 'none',
+    entries: Array.from({ length: 1024 }, (_, index) => [index + 1, 1]),
+    bucketCount: 1024,
+    relativeTolerance: exactTolerance,
+  },
+  {
+    name: 'approximation-threshold',
+    bins: 1000,
+    dequantize: 'none',
+    entries: [[1, 26239], [131071, 973761]],
+    options: { highest: 131071, figures: 5 },
+    bucketCount: 2,
+    relativeTolerance: asymptoticTolerance,
+  },
+  {
+    name: 'multimodal-none',
+    bins: 100,
+    dequantize: 'none',
+    entries: [[1, 900000], [10, 90000], [100, 9000], [1000, 1000]],
+    bucketCount: 4,
+    relativeTolerance: asymptoticTolerance,
+  },
+  {
+    name: 'multimodal-all',
+    bins: 100,
+    dequantize: 'all',
+    entries: [[1, 900000], [10, 90000], [100, 9000], [1000, 1000]],
+    bucketCount: 4,
+    relativeTolerance: asymptoticTolerance,
+  },
+  {
+    name: 'wide-none',
+    bins: 1000,
+    dequantize: 'none',
+    entries: [
+      [2 ** 20, 700000],
+      [2 ** 40, 200000],
+      [2 ** 52, 99999],
+      [Number.MAX_SAFE_INTEGER, 1],
+    ],
+    bucketCount: 4,
+    relativeTolerance: asymptoticTolerance,
+  },
+  {
+    name: 'wide-hdr',
+    bins: 1000,
+    dequantize: 'hdr',
+    entries: [
+      [2 ** 20, 700000],
+      [2 ** 40, 200000],
+      [2 ** 52, 99999],
+      [Number.MAX_SAFE_INTEGER, 1],
+    ],
+    bucketCount: 4,
+    relativeTolerance: asymptoticTolerance,
+  },
+  {
+    name: 'huge-count',
+    bins: 3,
+    dequantize: 'none',
+    entries: [[1, 2 ** 50], [1000, 2 ** 51]],
+    bucketCount: 2,
+    relativeTolerance: asymptoticTolerance,
+  },
+];
+
+(async () => {
+  assert.strictEqual(oracle.implementation, 'stats::pbeta/stats::dbeta');
+  for (const definition of cases) {
+    const histogram = buildHistogram(definition.entries, definition.options);
+    const options = { dequantize: definition.dequantize };
+    if (definition.probabilities === undefined) options.bins = definition.bins;
+    else options.probabilities = definition.probabilities;
+    const result = await histogram.qrde(options);
+    assertMatchesOracle(definition, result);
+  }
+})().then(common.mustCall());

test/parallel/test-perf-hooks-histogram-qrde.js

@@ -0,0 +1,222 @@
+'use strict';
+
+const common = require('../common');
+const assert = require('assert');
+const { createHistogram } = require('perf_hooks');
+
+function assertClose(actual, expected, tolerance = 1e-12) {
+  assert.ok(Math.abs(actual - expected) <= tolerance,
+            `${actual} != ${expected}`);
+}
+
+(async () => {
+  const empty = createHistogram();
+  const emptyResult = await empty.qrde();
+  assert.strictEqual(Object.getPrototypeOf(emptyResult), null);
+  const defaultProbabilities = new Float64Array(101);
+  for (let i = 0; i <= 100; i++) defaultProbabilities[i] = i / 100;
+  assert.deepStrictEqual(emptyResult.probabilities, defaultProbabilities);
+  assert.deepStrictEqual(emptyResult.quantiles, new Float64Array());
+  assert.deepStrictEqual(emptyResult.densities, new Float64Array());
+  assert.strictEqual(emptyResult.count, 0n);
+  assert.strictEqual(emptyResult.bucketCount, 0);
+  assert.strictEqual(emptyResult.corrections, 0);
+  assert.strictEqual(emptyResult.dequantize, 'hdr');
+
+  assert.throws(() => empty.qrde(null), {
+    code: 'ERR_INVALID_ARG_TYPE',
+  });
+  assert.throws(() => empty.qrde({ bins: 0 }), {
+    code: 'ERR_OUT_OF_RANGE',
+  });
+  assert.throws(() => empty.qrde({ bins: 1001 }), {
+    code: 'ERR_OUT_OF_RANGE',
+  });
+  assert.throws(() => empty.qrde({
+    bins: 10,
+    probabilities: [0, 1],
+  }), {
+    code: 'ERR_INVALID_ARG_VALUE',
+  });
+  assert.throws(() => empty.qrde({ probabilities: [0] }), {
+    code: 'ERR_INVALID_ARG_VALUE',
+  });
+  assert.throws(() => empty.qrde({ probabilities: [0.1, 1] }), {
+    code: 'ERR_INVALID_ARG_VALUE',
+  });
+  assert.throws(() => empty.qrde({ probabilities: [0, 0.9] }), {
+    code: 'ERR_INVALID_ARG_VALUE',
+  });
+  assert.throws(() => empty.qrde({ probabilities: [0, 0.5, 0.5, 1] }), {
+    code: 'ERR_INVALID_ARG_VALUE',
+  });
+  assert.throws(() => empty.qrde({ probabilities: [0, 1.1, 1] }), {
+    code: 'ERR_OUT_OF_RANGE',
+  });
+  assert.throws(() => empty.qrde({ probabilities: new Array(1002) }), {
+    code: 'ERR_OUT_OF_RANGE',
+  });
+  assert.throws(() => empty.qrde({ dequantize: true }), {
+    code: 'ERR_INVALID_ARG_VALUE',
+  });
+  assert.throws(() => empty.qrde({ cache: 1 }), {
+    code: 'ERR_INVALID_ARG_TYPE',
+  });
+
+  const sample = createHistogram();
+  sample.record(3);
+  sample.record(4);
+  sample.record(7);
+  const pending = sample.qrde({ bins: 10, dequantize: 'none' });
+  assert.ok(pending instanceof Promise);
+  const result = await pending;
+  assert.ok(result.quantiles instanceof Float64Array);
+  assert.ok(result.densities instanceof Float64Array);
+  const tenBinProbabilities = new Float64Array(11);
+  for (let i = 0; i <= 10; i++) tenBinProbabilities[i] = i / 10;
+  assert.deepStrictEqual(result.probabilities, tenBinProbabilities);
+  assert.strictEqual(result.quantiles.length, 11);
+  assert.strictEqual(result.densities.length, 10);
+  assert.strictEqual(result.count, 3n);
+  assert.strictEqual(result.bucketCount, 3);
+  assert.strictEqual(result.corrections, 0);
+  assert.strictEqual(result.dequantize, 'none');
+  assertClose(result.quantiles[5], 122 / 27);
+  for (let i = 0; i < result.densities.length; i++) {
+    const width = result.quantiles[i + 1] - result.quantiles[i];
+    assertClose(result.densities[i] * width, 0.1);
+  }
+
+  const customInput = [0, 0.2, 0.5, 0.9, 1];
+  const customPending = sample.qrde({
+    probabilities: customInput,
+    dequantize: 'none',
+  });
+  customInput[1] = 0.4;
+  const custom = await customPending;
+  assert.deepStrictEqual(custom.probabilities,
+                         new Float64Array([0, 0.2, 0.5, 0.9, 1]));
+  assert.strictEqual(custom.quantiles.length, 5);
+  assert.strictEqual(custom.densities.length, 4);
+  assertClose(custom.quantiles[1], result.quantiles[2]);
+  assertClose(custom.quantiles[2], result.quantiles[5]);
+  assertClose(custom.quantiles[3], result.quantiles[9]);
+  for (let i = 0; i < custom.densities.length; i++) {
+    const width = custom.quantiles[i + 1] - custom.quantiles[i];
+    const mass = custom.probabilities[i + 1] - custom.probabilities[i];
+    assertClose(custom.densities[i] * width, mass);
+  }
+
+  const pointMass = createHistogram();
+  for (let i = 0; i < 100; i++) pointMass.record(100);
+
+  const pointMassNone =
+    await pointMass.qrde({ bins: 10, dequantize: 'none' });
+  const pointMassHdr = await pointMass.qrde({ bins: 10 });
+  assert.deepStrictEqual(pointMassNone.quantiles,
+                         new Float64Array(11).fill(100));
+  assert.deepStrictEqual(pointMassHdr.quantiles, pointMassNone.quantiles);
+  assert.ok(pointMassNone.densities.every((value) => value === Infinity));
+  assert.ok(pointMassHdr.densities.every((value) => value === Infinity));
+
+  const pointMassAll =
+    await pointMass.qrde({ bins: 10, dequantize: 'all' });
+  for (let i = 0; i <= 10; i++) {
+    assertClose(pointMassAll.quantiles[i], 99.5 + i / 10);
+  }
+  for (const density of pointMassAll.densities) assertClose(density, 1);
+  assert.strictEqual(pointMassAll.corrections, 0);
+
+  const wideBucket = createHistogram();
+  for (let i = 0; i < 100; i++) wideBucket.record(100_000);
+  const [wideNone, wideHdr, wideAll] = await Promise.all([
+    wideBucket.qrde({ bins: 10, dequantize: 'none' }),
+    wideBucket.qrde({ bins: 10, dequantize: 'hdr' }),
+    wideBucket.qrde({ bins: 10, dequantize: 'all' }),
+  ]);
+  assert.ok(wideNone.densities.every((value) => value === Infinity));
+  assert.ok(wideHdr.densities.every(Number.isFinite));
+  assert.deepStrictEqual(wideHdr.quantiles, wideAll.quantiles);
+
+  const snapshot = createHistogram();
+  snapshot.record(3);
+  snapshot.record(4);
+  snapshot.record(7);
+  const snapshotPending =
+    snapshot.qrde({ bins: 1000, dequantize: 'none' });
+  snapshot.record(1000);
+  const snapshotResult = await snapshotPending;
+  assert.strictEqual(snapshotResult.count, 3n);
+  assert.strictEqual(snapshotResult.quantiles.at(-1), 7);
+
+  const cachedSnapshot = createHistogram();
+  cachedSnapshot.record(3);
+  cachedSnapshot.record(4);
+  cachedSnapshot.record(7);
+  const cachedPending = cachedSnapshot.qrde({ bins: 1000, cache: true });
+  cachedSnapshot.record(1000);
+  const cachedResult = await cachedPending;
+  assert.strictEqual(cachedResult.count, 3n);
+  assert.strictEqual(cachedResult.quantiles.at(-1), 7);
+  const refreshedCache = await cachedSnapshot.qrde({
+    probabilities: [0, 0.5, 1],
+    cache: true,
+  });
+  assert.strictEqual(refreshedCache.count, 4n);
+  assert.strictEqual(refreshedCache.quantiles.at(-1), 1000);
+  const reusedCache = await cachedSnapshot.qrde({ bins: 4, cache: true });
+  assert.strictEqual(reusedCache.count, 4n);
+  assert.strictEqual(reusedCache.quantiles.at(-1), 1000);
+  cachedSnapshot.reset();
+  const invalidatedCache = await cachedSnapshot.qrde({ cache: true });
+  assert.strictEqual(invalidatedCache.count, 0n);
+
+  const arithmeticCache = createHistogram();
+  arithmeticCache.record(1);
+  await arithmeticCache.qrde({ cache: true });
+  const operand = createHistogram();
+  operand.record(2);
+  arithmeticCache.add(operand);
+  assert.strictEqual(
+    (await arithmeticCache.qrde({ cache: true })).count, 2n);
+  arithmeticCache.subtract(operand);
+  assert.strictEqual(
+    (await arithmeticCache.qrde({ cache: true })).count, 1n);
+  arithmeticCache.recordCorrected(10, 5);
+  assert.strictEqual(
+    (await arithmeticCache.qrde({ cache: true })).count, 3n);
+
+  const exactTails = createHistogram();
+  for (let i = 1; i <= 1024; i++) exactTails.record(i);
+  const [exactNone, exactAll] = await Promise.all([
+    exactTails.qrde({ bins: 100, dequantize: 'none' }),
+    exactTails.qrde({ bins: 100, dequantize: 'all' }),
+  ]);
+  assert.strictEqual(exactNone.bucketCount, 1024);
+  assertClose(exactNone.quantiles[1], 10.740000001653616, 1e-9);
+  assertClose(exactNone.quantiles[99], 1014.2599999983464, 1e-9);
+  assertClose(exactAll.quantiles[1], exactNone.quantiles[1]);
+  assertClose(exactAll.quantiles[99], exactNone.quantiles[99]);
+
+  const tinyUpperTail = createHistogram({ highest: 32768, figures: 4 });
+  for (let i = 1; i <= 31; i++) tinyUpperTail.record(i);
+  tinyUpperTail.record(16384);
+  const tinyUpperTailResult =
+    await tinyUpperTail.qrde({ bins: 2, dequantize: 'none' });
+  assertClose(tinyUpperTailResult.quantiles[1],
+              16.500000000000938, 2e-13);
+
+  const largeCount = createHistogram();
+  largeCount.record(1);
+  largeCount.record(3);
+  for (let i = 0; i < 52; i++) {
+    const copy = createHistogram();
+    copy.add(largeCount);
+    largeCount.add(copy);
+  }
+  largeCount.record(1);
+  const largeCountResult =
+    await largeCount.qrde({ bins: 2, dequantize: 'none' });
+  assert.strictEqual(largeCountResult.count, (1n << 53n) + 1n);
+  assertClose(largeCountResult.quantiles[1], 2);
+})().then(common.mustCall());