T* histogram;

ASSIGN_OR_RETURN_UNWRAP(&histogram, receiver);

histogram->OnStart(reset ? T::StartFlags::RESET : T::StartFlags::NONE);

T* histogram;

ASSIGN_OR_RETURN_UNWRAP(&histogram, receiver);

histogram->OnStop();

: histogram_(std::move(histogram)) {

// alpha = 1 - 2^(-1/halfLife). With halfLife <= 0, EWMA is disabled.

if (options.half_life > 0) {

ewma_alpha_ = 1.0 - std::exp(-std::log(2.0) / options.half_life);

}

threshold_ = options.threshold;

hdr_histogram* histogram;

if (hdr_init(options.lowest, options.highest, options.figures, &histogram) !=

0) {

return {};

}

return std::make_shared<Histogram>(HistogramPointer(histogram), options);

CHECK_EQ(target->counts_len, source->counts_len);

target->min_value = source->min_value;

target->max_value = source->max_value;

target->normalizing_index_offset = source->normalizing_index_offset;

target->conversion_ratio = source->conversion_ratio;

target->total_count = source->total_count;

std::memcpy(target->counts,

source->counts,

source->counts_len * sizeof(*source->counts));

// The layout is fixed when the histogram is created, so it can be read

// without holding the lock.

hdr_histogram* histogram;

if (hdr_init(histogram_->lowest_discernible_value,

histogram_->highest_trackable_value,

histogram_->significant_figures,

&histogram) != 0) {

return {};

}

return std::make_shared<Histogram>(HistogramPointer(histogram), Options{});

std::shared_ptr<Histogram> clone = CreateWithSameLayout();

if (!clone) return {};

// Every member that holds recorded or statistical state must be copied

// here. The recorded snapshot cache is not copied; the clone builds its own

// on demand.

RwLock::ScopedReadLock lock(mutex_);

CopyRecordedData(clone->histogram_.get(), histogram_.get());

clone->prev_ = prev_;

clone->exceeds_ = exceeds_;

clone->reset_count_ = reset_count_;

clone->ewma_alpha_ = ewma_alpha_;

clone->ewma_mean_ = ewma_mean_;

clone->ewma_variance_ = ewma_variance_;

clone->ewma_initialized_ = ewma_initialized_;

clone->threshold_ = threshold_;

clone->ewma_error_rate_ = ewma_error_rate_;

return clone;

DiffError* error) const {

// Counts are subtracted index by index, so both histograms must map values

// to the same indexes. None of these fields change after creation.

if (!IsCompatible(other) || histogram_->normalizing_index_offset !=

other.histogram_->normalizing_index_offset) {

*error = DiffError::kIncompatible;

return {};

}

std::shared_ptr<Histogram> diff = CreateWithSameLayout();

if (!diff) {

*error = DiffError::kOutOfMemory;

return {};

}

// Only the recorded values and the exceeds count carry over. EWMA and timing

// state cannot be subtracted.

uint64_t reset_count;

{

RwLock::ScopedReadLock lock(mutex_);

CopyRecordedData(diff->histogram_.get(), histogram_.get());

diff->exceeds_ = exceeds_;

reset_count = reset_count_;

}

// `diff` is not shared yet, so only the lock of `other` is needed from here

// on. Never holding both locks at once avoids lock ordering issues.

RwLock::ScopedReadLock lock(other.mutex_);

if (reset_count != other.reset_count_) {

*error = DiffError::kReset;

return {};

}

if (diff->exceeds_ < other.exceeds_) {

*error = DiffError::kNotEarlier;

return {};

}

hdr_histogram* target = diff->histogram_.get();

const hdr_histogram* source = other.histogram_.get();

for (int32_t i = 0; i < target->counts_len; i++) {

if (target->counts[i] < source->counts[i]) {

*error = DiffError::kNotEarlier;

return {};

}

target->counts[i] -= source->counts[i];

}

diff->exceeds_ -= other.exceeds_;

hdr_reset_internal_counters(target);

*error = DiffError::kNone;

return diff;

tracker->TrackFieldWithSize("histogram", GetMemorySize());

tracker->TrackFieldWithSize("qrde_snapshot",

GetCachedRecordedSnapshotMemorySize());

return histogram_->counts_len == other.histogram_->counts_len &&

histogram_->lowest_discernible_value ==

other.histogram_->lowest_discernible_value &&

histogram_->highest_trackable_value ==

other.histogram_->highest_trackable_value &&

histogram_->significant_figures ==

other.histogram_->significant_figures;

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;

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;

}

CopyRecordedData(source.histogram.get(), histogram_.get());

return source;

uint64_t generation, std::shared_ptr<const RecordedSnapshot> snapshot) {

RwLock::ScopedWriteLock lock(mutex_);

if (generation == mutation_generation_) {

recorded_snapshot_cache_ = std::move(snapshot);

}

RwLock::ScopedReadLock lock(mutex_);

if (!recorded_snapshot_cache_) return 0;

return recorded_snapshot_cache_->buckets.capacity() * sizeof(RecordedBucket);

RwLock::ScopedReadLock lock(mutex_);

int64_t total = histogram_->total_count;

if (total == 0) return 0.0;

hdr_iter iter;

hdr_iter_init(&iter, histogram_.get());

while (hdr_iter_next(&iter)) {

if (iter.highest_equivalent_value >= value) {

return static_cast<double>(iter.cumulative_count) /

static_cast<double>(total);

}

// All recorded data accounted for; remaining buckets are empty.

if (iter.cumulative_count >= total) break;

}

return 1.0;

RwLock::ScopedReadLock lock(mutex_);

int64_t total = histogram_->total_count;

if (total < 3) return 0.0;

// Compute mean in one pass, then variance and skewness in a second

// pass. This avoids calling hdr_stddev (which internally recomputes

// hdr_mean), reducing the total from 4 iterations to 2.

double mean = hdr_mean(histogram_.get());

double m2 = 0.0;

double m3 = 0.0;

hdr_iter iter;

hdr_iter_recorded_init(&iter, histogram_.get());

while (hdr_iter_next(&iter)) {

double dev = static_cast<double>(hdr_median_equivalent_value(

histogram_.get(), iter.value)) -

mean;

double d2 = dev * dev;

m2 += static_cast<double>(iter.count) * d2;

m3 += static_cast<double>(iter.count) * d2 * dev;

}

double n = static_cast<double>(total);

double variance = m2 / n;

if (variance == 0.0) return 0.0;

double s3 = variance * std::sqrt(variance); // stddev^3

return (m3 / n) / s3;

RwLock::ScopedReadLock lock(mutex_);

int64_t total = histogram_->total_count;

if (total < 4) return 0.0;

// Same single-pass approach as Skewness: compute mean first, then

// variance and excess kurtosis together in one iteration.

double mean = hdr_mean(histogram_.get());

double m2 = 0.0;

double m4 = 0.0;

hdr_iter iter;

hdr_iter_recorded_init(&iter, histogram_.get());

while (hdr_iter_next(&iter)) {

double dev = static_cast<double>(hdr_median_equivalent_value(

histogram_.get(), iter.value)) -

mean;

double d2 = dev * dev;

m2 += static_cast<double>(iter.count) * d2;

m4 += static_cast<double>(iter.count) * d2 * d2;

}

double n = static_cast<double>(total);

double variance = m2 / n;

if (variance == 0.0) return 0.0;

double s4 = variance * variance; // stddev^4

return (m4 / n) / s4 - 3.0;

auto do_subtract = [&]() -> double {

int64_t dropped = 0;

int32_t len =

std::min(histogram_->counts_len, other.histogram_->counts_len);

for (int32_t i = 0; i < len; i++) {

int64_t count = histogram_->counts[i] - other.histogram_->counts[i];

if (count < 0) {

dropped += -count;

count = 0;

}

histogram_->counts[i] = count;

}

hdr_reset_internal_counters(histogram_.get());

InvalidateRecordedSnapshot();

reset_count_++;

exceeds_ = (exceeds_ > other.exceeds_) ? exceeds_ - other.exceeds_ : 0;

return static_cast<double>(dropped);

};

if (this == &other) {

RwLock::ScopedWriteLock lock(mutex_);

return do_subtract();

}

if (this < &other) {

RwLock::ScopedWriteLock lock1(mutex_);

RwLock::ScopedReadLock lock2(other.mutex_);

return do_subtract();

}

RwLock::ScopedReadLock lock1(other.mutex_);

RwLock::ScopedWriteLock lock2(mutex_);

return do_subtract();

auto do_ks = [&]() -> double {

int64_t n1 = histogram_->total_count;

int64_t n2 = other.histogram_->total_count;

if (n1 == 0 || n2 == 0) return 0.0;

double max_d = 0.0;

int64_t cum1 = 0, cum2 = 0;

int32_t len =

std::max(histogram_->counts_len, other.histogram_->counts_len);

for (int32_t i = 0; i < len; i++) {

if (i < histogram_->counts_len) cum1 += histogram_->counts[i];

if (i < other.histogram_->counts_len) cum2 += other.histogram_->counts[i];

double cdf1 = static_cast<double>(cum1) / static_cast<double>(n1);

double cdf2 = static_cast<double>(cum2) / static_cast<double>(n2);

double d = cdf1 > cdf2 ? cdf1 - cdf2 : cdf2 - cdf1;

if (d > max_d) max_d = d;

}

return max_d;

};

if (this == &other) return 0.0;

if (this < &other) {

RwLock::ScopedReadLock lock1(mutex_);

RwLock::ScopedReadLock lock2(other.mutex_);

return do_ks();

}

RwLock::ScopedReadLock lock1(other.mutex_);

RwLock::ScopedReadLock lock2(mutex_);

return do_ks();

int64_t* values,

size_t length) const {

RwLock::ScopedReadLock lock(mutex_);

hdr_value_at_percentiles(histogram_.get(), percentiles, values, length);

constexpr double FPMIN = 1e-30;

constexpr int MAXIT = 200;

constexpr double EPS = 3e-12;

double qab = a + b;

double qap = a + 1.0;

double qam = a - 1.0;

double c = 1.0;

double d = 1.0 - qab * x / qap;

if (std::fabs(d) < FPMIN) d = FPMIN;

d = 1.0 / d;

double h = d;

for (int m = 1; m <= MAXIT; m++) {

int m2 = 2 * m;

// Even step.

double aa = m * (b - m) * x / ((qam + m2) * (a + m2));

d = 1.0 + aa * d;

if (std::fabs(d) < FPMIN) d = FPMIN;

c = 1.0 + aa / c;

if (std::fabs(c) < FPMIN) c = FPMIN;

d = 1.0 / d;

h *= d * c;

// Odd step.

aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));

d = 1.0 + aa * d;

if (std::fabs(d) < FPMIN) d = FPMIN;

c = 1.0 + aa / c;

if (std::fabs(c) < FPMIN) c = FPMIN;

d = 1.0 / d;

double del = d * c;

h *= del;

if (std::fabs(del - 1.0) <= EPS) break;

}

return h;

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;

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))));

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,

};

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,

};

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);

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);

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;

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 < params.symmetry_point) {

return front * BetaContinuedFraction(params.a, params.b, x) / params.a;

}

return 1.0 -

front * BetaContinuedFraction(params.b, params.a, 1.0 - x) / params.b;

double squared = t * t;

double denominator = df + squared;

if (squared < df) {

double x = squared / denominator;

double ibeta = RegularizedIncompleteBeta(0.5, df / 2.0, x);

return t >= 0.0 ? 0.5 + 0.5 * ibeta : 0.5 - 0.5 * ibeta;

}

double x = df / denominator;

double ibeta = RegularizedIncompleteBeta(df / 2.0, 0.5, x);

return t >= 0.0 ? 1.0 - 0.5 * ibeta : 0.5 * ibeta;

double lo = 0.0;

double hi = 1.0;

while (StudentTCdf(-hi, df) > p) hi *= 2.0;

for (int i = 0; i < 100; i++) {

double mid = (lo + hi) / 2.0;

if (StudentTCdf(-mid, df) > p) {

lo = mid;

} else {

hi = mid;

}

}

return (lo + hi) / 2.0;

if (k < 0) return 0.0;

if (k >= n) return 1.0;

return RegularizedIncompleteBeta(

static_cast<double>(n - k), static_cast<double>(k + 1), 1.0 - p);

std::vector<double> quantiles;

std::vector<double> densities;

size_t corrections = 0;

kNone,

kHdr,

kAll,

QrdeDequantization dequantization) {

return bucket.count > 1 && (dequantization == QrdeDequantization::kAll ||

(dequantization == QrdeDequantization::kHdr &&

bucket.resolution > 1.0));

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};

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;

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;

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;

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());

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()) {