diff --git a/src/tdigest.c b/src/tdigest.c index 1268c15..f04398a 100644 --- a/src/tdigest.c +++ b/src/tdigest.c @@ -49,55 +49,175 @@ static void inline swap_l(long long *arr, int i, int j) { arr[j] = temp; } -static unsigned int partition(double *means, long long *weights, unsigned int start, - unsigned int end, unsigned int pivot_idx) { - const double pivotMean = means[pivot_idx]; - swap(means, pivot_idx, end); - swap_l(weights, pivot_idx, end); - - int i = start - 1; - - for (unsigned int j = start; j < end; j++) { - // If current element is smaller than the pivot - if (means[j] < pivotMean) { - // increment index of smaller element - i++; - swap(means, i, j); - swap_l(weights, i, j); +// Optional key-comparison counter, used only by the complexity regression test. +// Zero-cost unless TD_INSTRUMENT_SORT is defined at build time. +#ifdef TD_INSTRUMENT_SORT +unsigned long long td_sort_comparisons = 0; +unsigned long long td_sort_heap_fallbacks = 0; +#define TD_SORT_CMP() (++td_sort_comparisons) +#define TD_SORT_FALLBACK() (++td_sort_heap_fallbacks) +#else +#define TD_SORT_CMP() ((void)0) +#define TD_SORT_FALLBACK() ((void)0) +#endif + +// Counted key comparison `a < b`. Keeping the counter inside a dedicated helper (instead of a +// comma expression in the && operands) leaves the sort's boolean conditions side-effect free; +// it inlines to a plain `a < b` when TD_INSTRUMENT_SORT is off. +static inline bool td_key_lt(double a, double b) { + TD_SORT_CMP(); + return a < b; +} + +#define TD_INSORT_THRESHOLD 16 + +// Insertion sort of the inclusive range [lo, hi] (used for small ranges). +static void td_insertion_sort(double *means, long long *weights, int lo, int hi) { + for (int i = lo + 1; i <= hi; i++) { + const double m = means[i]; + const long long w = weights[i]; + int j = i - 1; + while (j >= lo && td_key_lt(m, means[j])) { + means[j + 1] = means[j]; + weights[j + 1] = weights[j]; + j--; } + means[j + 1] = m; + weights[j + 1] = w; } - swap(means, i + 1, end); - swap_l(weights, i + 1, end); - return i + 1; +} + +// Max-heap sift-down over the range starting at `lo` (heap size `n`, root heap +// index `i`), keyed on means with weights moved in lock-step. +static void td_sift_down(double *means, long long *weights, int lo, int i, int n) { + // Internal nodes are [0, n/2); index >= n/2 is a leaf. Testing leaf-ness BEFORE forming + // 2*i+1 keeps the child index from overflowing signed int when n approaches INT_MAX (a + // capacity #41 permits): for i < n/2, 2*i < n so 2*i+1 <= n and cannot overflow. + while (i < n / 2) { + int child = 2 * i + 1; + if (child + 1 < n && td_key_lt(means[lo + child], means[lo + child + 1])) { + child++; + } + if (!td_key_lt(means[lo + i], means[lo + child])) { + break; + } + swap(means, lo + i, lo + child); + swap_l(weights, lo + i, lo + child); + i = child; + } +} + +// Heapsort of the inclusive range [lo, hi]: the O(n log n) worst-case fallback. +static void td_heap_sort(double *means, long long *weights, int lo, int hi) { + const int n = hi - lo + 1; + for (int i = n / 2 - 1; i >= 0; i--) { + td_sift_down(means, weights, lo, i, n); + } + for (int end = n - 1; end > 0; end--) { + swap(means, lo, lo + end); + swap_l(weights, lo, lo + end); + td_sift_down(means, weights, lo, 0, end); + } +} + +// Move the median of means[lo], means[mid], means[hi] into `mid` (weights follow) +// and return it, so an already-sorted / organ-pipe input does not pick a bad pivot. +static double td_median3(double *means, long long *weights, int lo, int mid, int hi) { + if (td_key_lt(means[mid], means[lo])) { + swap(means, lo, mid); + swap_l(weights, lo, mid); + } + if (td_key_lt(means[hi], means[lo])) { + swap(means, lo, hi); + swap_l(weights, lo, hi); + } + if (td_key_lt(means[hi], means[mid])) { + swap(means, mid, hi); + swap_l(weights, mid, hi); + } + return means[mid]; } /** - * Standard quick sort except that sorting rearranges parallel arrays + * Introsort over two parallel arrays keyed on `means` (weights permuted in + * lock-step). Guaranteed O(n log n) time and O(log n) stack: + * - 3-way (Dutch-flag) partitioning collapses runs of equal keys in a single + * pass -- t-digest's common duplicate-heavy input; + * - a median-of-three pivot avoids the trivial already-sorted worst case; + * - a recursion-depth limit falls back to heapsort, so an adversarial + * distinct-value permutation crafted to defeat the pivot (which a + * fixed-position pivot cannot resist) cannot force quadratic time; + * - small ranges finish with insertion sort; + * - recursing only the smaller side bounds the stack to O(log n). * - * @param means Values to sort on - * @param weights The auxillary values to sort. - * @param start The beginning of the values to sort - * @param end The value after the last value to sort + * Note: this is NOT output-identical to a naive single-pivot sort for runs of + * equal keys -- the order within an equal-key run differs, which can change how + * weighted duplicate centroids are subsequently merged. Results stay within + * t-digest's accuracy guarantee; exact centroid-for-centroid reproduction of the + * old sorter is not preserved. */ -static void td_qsort(double *means, long long *weights, unsigned int start, unsigned int end) { - if (start < end) { - // two elements can be directly compared - if ((end - start) == 1) { - if (means[start] > means[end]) { - swap(means, start, end); - swap_l(weights, start, end); - } +static void td_introsort(double *means, long long *weights, int lo, int hi, int depth_limit) { + while (hi - lo > TD_INSORT_THRESHOLD) { + if (depth_limit == 0) { + TD_SORT_FALLBACK(); + td_heap_sort(means, weights, lo, hi); return; } - // generating a random number as a pivot was very expensive vs the array size - // const unsigned int pivot_idx = start + rand()%(end - start + 1); - const unsigned int pivot_idx = (end + start) / 2; // central pivot - const unsigned int new_pivot_idx = partition(means, weights, start, end, pivot_idx); - if (new_pivot_idx > start) { - td_qsort(means, weights, start, new_pivot_idx - 1); + depth_limit--; + const int mid = lo + (hi - lo) / 2; + const double pivot = td_median3(means, weights, lo, mid, hi); + // While scanning: [lo, lt) < pivot, [lt, i) == pivot, (gt, hi] > pivot. + int lt = lo; + int i = lo; + int gt = hi; + while (i <= gt) { + const double v = means[i]; + if (td_key_lt(v, pivot)) { + swap(means, i, lt); + swap_l(weights, i, lt); + lt++; + i++; + } else if (td_key_lt(pivot, v)) { + swap(means, i, gt); + swap_l(weights, i, gt); + gt--; + } else { + i++; + } + } + const int left_size = lt - lo; // count of elements < pivot + const int right_size = hi - gt; // count of elements > pivot + // Recurse the smaller side, loop on the larger (bounds the stack). + if (left_size < right_size) { + if (left_size > 1) { + td_introsort(means, weights, lo, lt - 1, depth_limit); + } + lo = gt + 1; + } else { + if (right_size > 1) { + td_introsort(means, weights, gt + 1, hi, depth_limit); + } + hi = lt - 1; } - td_qsort(means, weights, new_pivot_idx + 1, end); } + td_insertion_sort(means, weights, lo, hi); +} + +static void td_qsort(double *means, long long *weights, unsigned int lo_u, unsigned int hi_u) { + // Indices fit an int: the node arrays are sized by `cap` (an int field), so + // the sorted range is always within [0, node_count). + const int lo = (int)lo_u; + const int hi = (int)hi_u; + if (lo >= hi) { + return; + } + // Depth limit = 2*floor(log2(n)); exceeding it hands the range to heapsort, + // which is what guarantees the O(n log n) worst case. + int depth_limit = 0; + for (int t = hi - lo + 1; t > 1; t >>= 1) { + depth_limit += 2; + } + td_introsort(means, weights, lo, hi, depth_limit); } static inline uint64_t cap_from_compression(uint64_t compression) { @@ -568,6 +688,15 @@ double td_trimmed_mean(td_histogram_t *h, double leftmost_cut, double rightmost_ } int td_add(td_histogram_t *h, double mean, long long weight) { + // Reject non-finite means before any mutation. NaN has no ordering, so it would leave a + // partition unsorted and violate td_compress()'s sorted invariant. +/-Inf sorts fine but is + // not closed under the centroid-merge arithmetic: merging two equal infinities computes + // `delta = Inf - Inf = NaN`, poisoning a centroid mean (which then breaks a later sort). + // Rejecting all non-finite input at ingest keeps every stored mean finite, matching the + // reference t-digest (which rejects NaN in add()). + if (!isfinite(mean)) { + return EINVAL; + } if (should_td_compress(h)) { const int overflow_res = td_compress(h); if (overflow_res != 0) diff --git a/src/tdigest.h b/src/tdigest.h index c07436c..f04eb18 100644 --- a/src/tdigest.h +++ b/src/tdigest.h @@ -100,10 +100,12 @@ void td_reset(td_histogram_t *h); /** * Adds a sample to a histogram. * - * @param val The value to add. + * @param val The value to add. Must be finite: NaN has no ordering (would break the centroid + * sort) and +/-Inf is not closed under the centroid-merge arithmetic (Inf-Inf yields NaN), so + * both are rejected. * @param weight The weight of this point. - * @return 0 on success, EDOM if overflow was detected as a consequence of adding the provided - * weight. + * @return 0 on success, EINVAL if val is not finite (NaN or +/-Inf), EDOM if overflow was + * detected as a consequence of adding the provided weight. * */ int td_add(td_histogram_t *h, double val, long long weight); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 85da305..fe43943 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -23,6 +23,15 @@ if (BUILD_TESTS) target_link_libraries(td_test tdigest m) enable_testing() add_test(td_test td_test) + + # Complexity regression: compiles the library with TD_INSTRUMENT_SORT to count + # key comparisons and asserts one large compress stays O(n log n) (MOD-17228). + add_executable(td_sort_complexity_test unit/td_sort_complexity_test.c) + # CMAKE_CURRENT_LIST_DIR so the src path stays correct under add_subdirectory(). + target_include_directories(td_sort_complexity_test PRIVATE ${CMAKE_CURRENT_LIST_DIR}/../src) + target_compile_definitions(td_sort_complexity_test PRIVATE TD_INSTRUMENT_SORT) + target_link_libraries(td_sort_complexity_test m) + add_test(td_sort_complexity_test td_sort_complexity_test) endif() diff --git a/tests/unit/td_sort_complexity_test.c b/tests/unit/td_sort_complexity_test.c new file mode 100644 index 0000000..1018500 --- /dev/null +++ b/tests/unit/td_sort_complexity_test.c @@ -0,0 +1,321 @@ +/* + * Complexity regression for the centroid sort (MOD-17228). + * + * A correctness-only test cannot tell an O(n log n) sort from an O(n^2) one, and the ordinary + * unit test auto-compresses in small batches so it never sorts one large partition. This test + * compiles the library with TD_INSTRUMENT_SORT (key-comparison and heapsort-fallback counters) + * and asserts the two properties that give td_qsort its worst-case guarantee, each driven + * through the real sort by the input that actually defeats a sort lacking it, across multiple + * sizes so a quadratic implementation cannot slip under a single hand-picked bound: + * + * Part A - duplicate-heavy input (the reported DoS). N identical values is the MOD-17228 + * vector: the pre-PR central-pivot Lomuto sort peels one element per level and does + * exactly n^2/2 comparisons (~5e9 at n=1e5, a ~98 s server freeze). The 3-way + * partition collapses the equal run in a single pass -> exactly 2n comparisons. + * We assert a strict LINEAR bound, which a revert to any sort without equal-run + * handling misses by three-plus orders of magnitude. + * + * Part B - the heapsort fallback, through the FULL td_qsort() path. A median-of-three + * quicksort has an adversarial distinct permutation that forces Theta(n^2); the + * depth-limit -> heapsort fallback is what caps it at O(n log n). We generate that + * permutation with McIlroy's quicksort adversary (see gen_killer) run against a + * faithful mirror of this introsort's pivot/partition sequence, feed it through the + * real td_compress()/td_qsort(), and assert three things across sizes: the heapsort + * fallback is actually reached (td_sort_heap_fallbacks > 0), the comparison count + * stays within an O(n log n) bound (measured ~5.7 n log2 n; without the fallback the + * same input is ~n^2/2, i.e. ~3000 n log2 n at n=1e5), and the normalized ratio + * stays flat. Removing or weakening the fallback fails both the fallback-reached and + * the bound assertion. + * + * Numbers above are measured, not assumed; the asserted constants leave generous margin over + * them while staying far below quadratic. + */ +#ifndef TD_INSTRUMENT_SORT +/* Also set by target_compile_definitions() in tests/CMakeLists.txt; guard so a standalone + * compile still instruments without redefining the macro (which -Werror would reject). */ +#define TD_INSTRUMENT_SORT 1 +#endif + +#include +#include +#include + +#include "tdigest.c" /* brings in the static sort helpers + instrumentation counters */ + +static int failures = 0; + +#define CHECK(cond, ...) \ + do { \ + if (!(cond)) { \ + fprintf(stderr, "FAIL: "); \ + fprintf(stderr, __VA_ARGS__); \ + fprintf(stderr, "\n"); \ + failures++; \ + } \ + } while (0) + +/* ---------- Part A ---------- */ + +/* Sort n identical values in a single compress and return the comparison count. */ +static double compares_all_equal(int n) { + /* cap = 6n + 10 > n, so all n points buffer and one td_compress() sorts them together. */ + td_histogram_t *h = td_new((double)n); + if (h == NULL) { + fprintf(stderr, "allocation failed at n=%d\n", n); + exit(1); + } + for (int i = 0; i < n; ++i) { + if (td_add(h, 42.0, 1) != 0) { + fprintf(stderr, "td_add failed at %d\n", i); + exit(1); + } + } + if (h->unmerged_nodes != n) { + fprintf(stderr, "expected one big compress, but auto-compress ran (unmerged=%d)\n", + h->unmerged_nodes); + exit(1); + } + td_sort_comparisons = 0; + if (td_compress(h) != 0) { + fprintf(stderr, "compress failed at n=%d\n", n); + exit(1); + } + const double c = (double)td_sort_comparisons; + for (int i = 1; i < h->merged_nodes; ++i) { + CHECK(h->nodes_mean[i - 1] <= h->nodes_mean[i], "A: centroids not sorted at %d (n=%d)", i, + n); + } + td_free(h); + return c; +} + +/* ---------- Part B: McIlroy quicksort adversary tailored to this introsort ---------- + * + * The adversary keeps every key "gas" (unassigned) and freezes one only when a comparison + * forces it, always making the pivot just chosen an extreme -> maximally unbalanced partitions. + * It must be driven by the SAME comparison sequence as the target sort, so mirror_qsort below is + * a faithful copy of td_introsort's median-of-three + 3-way partition + recurse-smaller-side + * logic (without the depth limit, so generation elicits the quadratic path). The resulting value + * assignment is the killer permutation for the real td_qsort(). If td_introsort's pivot or + * partition strategy changes, this mirror must change with it. */ +static int *g_val; /* assigned key per identity, or g_gas if still unassigned */ +static int g_gas; /* sentinel meaning "unassigned" */ +static int g_nsolid; /* number of frozen keys */ +static int g_cand; /* current pivot candidate */ + +static int adv_cmp(int x, int y) { /* <0 if key(x)0 if greater, 0 if equal */ + if (g_val[x] == g_gas && g_val[y] == g_gas) { + if (x == g_cand) { + g_val[x] = g_nsolid++; + } else { + g_val[y] = g_nsolid++; + } + } + if (g_val[x] == g_gas) { + g_cand = x; + return 1; + } + if (g_val[y] == g_gas) { + g_cand = y; + return -1; + } + return g_val[x] - g_val[y]; +} + +static void idx_swap(int *a, int i, int j) { + const int t = a[i]; + a[i] = a[j]; + a[j] = t; +} + +/* mirror of td_insertion_sort / td_median3 / td_introsort (see TD_INSORT_THRESHOLD). */ +static void mirror_insertion(int *a, int lo, int hi) { + for (int i = lo + 1; i <= hi; i++) { + const int m = a[i]; + int j = i - 1; + while (j >= lo && adv_cmp(m, a[j]) < 0) { + a[j + 1] = a[j]; + j--; + } + a[j + 1] = m; + } +} + +static int mirror_median3(int *a, int lo, int mid, int hi) { + if (adv_cmp(a[mid], a[lo]) < 0) { + idx_swap(a, lo, mid); + } + if (adv_cmp(a[hi], a[lo]) < 0) { + idx_swap(a, lo, hi); + } + if (adv_cmp(a[hi], a[mid]) < 0) { + idx_swap(a, mid, hi); + } + return a[mid]; +} + +static void mirror_qsort(int *a, int lo, int hi) { + while (hi - lo > TD_INSORT_THRESHOLD) { + const int mid = lo + (hi - lo) / 2; + const int pivot = mirror_median3(a, lo, mid, hi); + int lt = lo; + int i = lo; + int gt = hi; + while (i <= gt) { + const int c = adv_cmp(a[i], pivot); + if (c < 0) { + idx_swap(a, i, lt); + lt++; + i++; + } else if (c > 0) { + idx_swap(a, i, gt); + gt--; + } else { + i++; + } + } + const int left_size = lt - lo; + const int right_size = hi - gt; + if (left_size < right_size) { + if (left_size > 1) { + mirror_qsort(a, lo, lt - 1); + } + lo = gt + 1; + } else { + if (right_size > 1) { + mirror_qsort(a, gt + 1, hi); + } + hi = lt - 1; + } + } + mirror_insertion(a, lo, hi); +} + +/* Fill out[0..n-1] with the killer permutation for the real td_qsort. */ +static void gen_killer(int n, double *out) { + g_val = (int *)calloc((size_t)n, sizeof(int)); + int *a = (int *)calloc((size_t)n, sizeof(int)); + if (g_val == NULL || a == NULL) { + fprintf(stderr, "allocation failed at n=%d\n", n); + exit(1); + } + g_gas = n; + g_nsolid = 0; + g_cand = 0; + for (int i = 0; i < n; i++) { + g_val[i] = g_gas; + a[i] = i; + } + mirror_qsort(a, 0, n - 1); + for (int i = 0; i < n; i++) { + if (g_val[i] == g_gas) { /* never compared: assign any remaining rank */ + g_val[i] = g_nsolid++; + } + out[i] = (double)g_val[i]; + } + free(a); + free(g_val); + g_val = NULL; +} + +/* Sort the killer through the real td_compress()/td_qsort(); report comparisons + fallbacks. */ +static double compares_killer(int n, unsigned long long *fallbacks_out) { + double *killer = (double *)calloc((size_t)n, sizeof(double)); + if (killer == NULL) { + fprintf(stderr, "allocation failed at n=%d\n", n); + exit(1); + } + gen_killer(n, killer); + + td_histogram_t *h = td_new((double)n); + if (h == NULL) { + fprintf(stderr, "allocation failed at n=%d\n", n); + exit(1); + } + for (int i = 0; i < n; ++i) { + if (td_add(h, killer[i], 1) != 0) { + fprintf(stderr, "td_add failed at %d\n", i); + exit(1); + } + } + free(killer); + if (h->unmerged_nodes != n) { + fprintf(stderr, "expected one big compress, but auto-compress ran (unmerged=%d)\n", + h->unmerged_nodes); + exit(1); + } + td_sort_comparisons = 0; + td_sort_heap_fallbacks = 0; + if (td_compress(h) != 0) { + fprintf(stderr, "compress failed at n=%d\n", n); + exit(1); + } + const double c = (double)td_sort_comparisons; + *fallbacks_out = td_sort_heap_fallbacks; + for (int i = 1; i < h->merged_nodes; ++i) { + CHECK(h->nodes_mean[i - 1] <= h->nodes_mean[i], "B: centroids not sorted at %d (n=%d)", i, + n); + } + td_free(h); + return c; +} + +int main(void) { + /* Part A is O(n) so it runs at large sizes. Part B's killer GENERATION runs the adversarial + * mirror sort, which is intentionally O(n^2), so it uses modest sizes; the fallback engages + * and the bounds separate quadratic from O(n log n) well before n gets large. */ + const int a_sizes[] = {25000, 50000, 100000, 200000}; + const int b_sizes[] = {3000, 6000, 12000, 24000}; + const int na = (int)(sizeof(a_sizes) / sizeof(a_sizes[0])); + const int nb = (int)(sizeof(b_sizes) / sizeof(b_sizes[0])); + + /* Part A: duplicate-heavy input must stay LINEAR (3-way partition). New sort = 2n; the + * pre-PR Lomuto sort = n^2/2. Bound 8n leaves 4x margin and is far below quadratic. */ + printf("Part A - duplicate-heavy (MOD-17228 DoS), must be linear:\n"); + for (int i = 0; i < na; ++i) { + const int n = a_sizes[i]; + const double c = compares_all_equal(n); + const double linear_bound = 8.0 * (double)n; + const double quadratic = 0.5 * (double)n * (double)n; + printf(" n=%7d cmps=%12.0f c/n=%5.2f bound(8n)=%.0f (pre-PR ~n^2/2=%.0f)\n", n, c, + c / n, linear_bound, quadratic); + CHECK(c <= linear_bound, "A: n=%d comparisons %.0f exceed linear bound %.0f (quadratic?)", + n, c, linear_bound); + } + + /* Part B: adversarial distinct permutation through the FULL td_qsort() path. The fallback + * must be reached and cap the work at O(n log n) (measured ~5.7 n log2 n); the SAME input + * without the fallback is ~n^2/2 (~3000 n log2 n at n=1e5). Bound 12 n log2 n. */ + printf("Part B - median-of-3 killer through full td_qsort(), fallback must engage:\n"); + double max_ratio = 0.0; + double min_ratio = 1e300; + for (int i = 0; i < nb; ++i) { + const int n = b_sizes[i]; + unsigned long long fallbacks = 0; + const double c = compares_killer(n, &fallbacks); + const double nlogn = (double)n * log2((double)n); + const double ratio = c / nlogn; + if (ratio > max_ratio) { + max_ratio = ratio; + } + if (ratio < min_ratio) { + min_ratio = ratio; + } + printf(" n=%7d cmps=%12.0f c/(n*log2n)=%5.3f fallbacks=%llu bound(12*n*log2n)=%.0f\n", + n, c, ratio, fallbacks, 12.0 * nlogn); + CHECK(fallbacks > 0, "B: n=%d heapsort fallback was never reached (killer ineffective?)", + n); + CHECK(c <= 12.0 * nlogn, "B: n=%d comparisons %.0f exceed O(n log n) bound %.0f", n, c, + 12.0 * nlogn); + } + /* Flatness: for O(n log n) the ratio is ~constant; a quadratic path would grow it ~n/log n. */ + CHECK(max_ratio <= 2.0 * min_ratio, "B: normalized comparison ratio not flat (%.3f..%.3f)", + min_ratio, max_ratio); + + if (failures == 0) { + printf("OK\n"); + return 0; + } + fprintf(stderr, "%d complexity check(s) failed\n", failures); + return 1; +} diff --git a/tests/unit/td_test.c b/tests/unit/td_test.c index 5803e72..66e6366 100644 --- a/tests/unit/td_test.c +++ b/tests/unit/td_test.c @@ -414,6 +414,40 @@ MU_TEST(test_nans) { td_free(t); } +// td_add() rejects every non-finite mean. NaN has no ordering for the centroid sort, and +/-Inf +// is not closed under the centroid-merge arithmetic (merging two equal infinities computes +// Inf - Inf = NaN, poisoning a centroid). Repeated-infinity input previously produced NaN +// centroids; rejecting it at ingest keeps every stored mean finite and the sort invariant intact. +MU_TEST(test_add_nonfinite) { + td_histogram_t *t = td_new(200); + mu_assert(td_add(t, NAN, 1) == EINVAL, "td_add(NaN) must be rejected with EINVAL"); + mu_assert(td_add(t, INFINITY, 1) == EINVAL, "td_add(+Inf) must be rejected with EINVAL"); + mu_assert(td_add(t, -INFINITY, 1) == EINVAL, "td_add(-Inf) must be rejected with EINVAL"); + mu_assert(td_centroid_count(t) == 0, "rejected non-finite input must not be stored"); + + // The old repeated-+Inf reproducer (100 inserts -> NaN centroids) must now add nothing and + // leave a clean, empty digest. + for (int i = 0; i < 100; ++i) { + mu_assert(td_add(t, INFINITY, 1) == EINVAL, "repeated +Inf still rejected"); + } + mu_assert(td_centroid_count(t) == 0, "repeated +Inf must leave the digest empty"); + + // Finite values still work and stay sorted / NaN-free after a compress. + for (int i = 0; i < 50; ++i) { + mu_assert(td_add(t, (double)(i - 25), 1) == 0, "finite insertion"); + } + mu_assert(td_compress(t) == 0, "compress finite values"); + const long long n = td_centroid_count(t); + for (long long i = 0; i < n; ++i) { + mu_assert(isfinite(td_centroids_mean_at(t, (int)i)), "every centroid mean must be finite"); + if (i > 0) { + mu_assert(td_centroids_mean_at(t, (int)(i - 1)) <= td_centroids_mean_at(t, (int)i), + "centroids must stay sorted"); + } + } + td_free(t); +} + MU_TEST(test_two_interp) { td_histogram_t *t = td_new(1000); mu_assert(td_add(t, 1, 1) == 0, "Insertion"); @@ -594,12 +628,105 @@ MU_TEST(test_quantiles_multiple) { td_free(t); } +// Exercises the centroid sort on the inputs that were pathological for the old +// single-pivot quicksort: all-equal, ascending, and descending. The 3-way sort +// must produce a correctly ordered, correct-weight digest for each. (The old +// sort was O(n^2)/O(n)-stack on the all-equal case; this guards correctness of +// the replacement.) +MU_TEST(test_duplicate_heavy_compress) { + const int n = 50000; + + // All-equal: everything collapses onto a single value. + td_histogram_t *eq = td_new(200); + mu_assert(eq != NULL, "created_histogram"); + for (int i = 0; i < n; ++i) { + mu_assert(td_add(eq, 42.0, 1) == 0, "Insertion"); + } + mu_assert(td_compress(eq) == 0, "compress all-equal"); + mu_assert_double_eq((double)n, td_size(eq)); + mu_assert_double_eq(42.0, td_min(eq)); + mu_assert_double_eq(42.0, td_max(eq)); + mu_assert_double_eq(42.0, td_quantile(eq, 0.0)); + mu_assert_double_eq(42.0, td_quantile(eq, 0.5)); + mu_assert_double_eq(42.0, td_quantile(eq, 1.0)); + // Merged centroid means must be non-decreasing. + for (int i = 1; i < eq->merged_nodes; ++i) { + mu_assert(eq->nodes_mean[i - 1] <= eq->nodes_mean[i], "means sorted (all-equal)"); + } + td_free(eq); + + // Ascending and descending must yield the same digest bounds. + td_histogram_t *asc = td_new(200); + td_histogram_t *desc = td_new(200); + mu_assert(asc != NULL && desc != NULL, "created_histograms"); + for (int i = 0; i < n; ++i) { + mu_assert(td_add(asc, (double)i, 1) == 0, "Insertion asc"); + mu_assert(td_add(desc, (double)(n - 1 - i), 1) == 0, "Insertion desc"); + } + mu_assert(td_compress(asc) == 0, "compress asc"); + mu_assert(td_compress(desc) == 0, "compress desc"); + mu_assert_double_eq(0.0, td_min(asc)); + mu_assert_double_eq((double)(n - 1), td_max(asc)); + mu_assert_double_eq(0.0, td_min(desc)); + mu_assert_double_eq((double)(n - 1), td_max(desc)); + for (int i = 1; i < asc->merged_nodes; ++i) { + mu_assert(asc->nodes_mean[i - 1] <= asc->nodes_mean[i], "means sorted (asc)"); + } + for (int i = 1; i < desc->merged_nodes; ++i) { + mu_assert(desc->nodes_mean[i - 1] <= desc->nodes_mean[i], "means sorted (desc)"); + } + // The median of a dense uniform 0..n-1 sits near the middle for both orders. + mu_assert_double_eq_epsilon((double)(n - 1) / 2.0, td_quantile(asc, 0.5), (double)n * 0.02); + mu_assert_double_eq_epsilon((double)(n - 1) / 2.0, td_quantile(desc, 0.5), (double)n * 0.02); + td_free(asc); + td_free(desc); +} + +// Weighted duplicates exercise how runs of equal keys are merged after sorting. +// The introsort changes the intra-run order relative to a naive single-pivot +// sort (so the digest is NOT centroid-for-centroid identical), but the result +// must stay accurate. Uses a distribution whose quantiles are known exactly. +MU_TEST(test_weighted_duplicates_accuracy) { + td_histogram_t *t = td_new(200); + mu_assert(t != NULL, "created_histogram"); + // 5 distinct values, each carrying equal weight, inserted as several weighted + // duplicates so the same mean appears in multiple centroids before merging. + const double vals[5] = {1.0, 2.0, 3.0, 4.0, 5.0}; + long long total = 0; + for (int rep = 0; rep < 4; ++rep) { + for (int i = 0; i < 5; ++i) { + mu_assert(td_add(t, vals[i], 150) == 0, "weighted duplicate insertion"); + total += 150; + } + } + mu_assert(td_compress(t) == 0, "compress"); + mu_assert_double_eq(1.0, td_min(t)); + mu_assert_double_eq(5.0, td_max(t)); + mu_assert_long_eq(total, td_size(t)); // 3000, weight fully accounted + // Each value holds exactly 20% of the mass, so the median is 3 and the + // quantiles land on the values (within interpolation tolerance). + mu_assert_double_eq_epsilon(1.0, td_quantile(t, 0.05), 0.6); + mu_assert_double_eq_epsilon(3.0, td_quantile(t, 0.5), 0.6); + mu_assert_double_eq_epsilon(5.0, td_quantile(t, 0.95), 0.6); + // CDF stays in [0,1] and non-decreasing across the support. + double prev = -1.0; + for (int step = 0; step <= 12; ++step) { + const double x = 0.5 * (double)step; + const double c = td_cdf(t, x); + mu_assert(c >= 0.0 && c <= 1.0, "cdf within [0,1]"); + mu_assert(c >= prev - 1e-9, "cdf non-decreasing"); + prev = c; + } + td_free(t); +} + MU_TEST_SUITE(test_suite) { MU_RUN_TEST(test_basic); MU_RUN_TEST(test_td_init); MU_RUN_TEST(test_compress_small); MU_RUN_TEST(test_compress_large); MU_RUN_TEST(test_nans); + MU_RUN_TEST(test_add_nonfinite); MU_RUN_TEST(test_negative_values); MU_RUN_TEST(test_negative_values_merge); MU_RUN_TEST(test_large_outlier_test); @@ -615,6 +742,8 @@ MU_TEST_SUITE(test_suite) { MU_RUN_TEST(test_trimmed_mean_complex); MU_RUN_TEST(test_overflow); MU_RUN_TEST(test_overflow_merge); + MU_RUN_TEST(test_duplicate_heavy_compress); + MU_RUN_TEST(test_weighted_duplicates_accuracy); } int main(int argc, char *argv[]) {