diff --git a/crates/consistent-choose-k/README.md b/crates/consistent-choose-k/README.md index 45574ab..fd9c105 100644 --- a/crates/consistent-choose-k/README.md +++ b/crates/consistent-choose-k/README.md @@ -37,6 +37,33 @@ Why replication matters - Distributes read/write load across multiple owners, reducing hotspots. - Enables fast recovery and higher tail-latency resilience. +## Permutation APIs and membership semantics + +The existing `ConsistentPermutation` preserves **survivor list order** when +nodes are appended or removed from the end of `0..n`. The experimental +`VirtualPermutation` now supplies the same order-restriction contract using a +different construction: it traverses a single consistent cycle from a permanent +sentinel. Removing the appended node from the larger complete order recovers +the smaller complete order. Both yield distinct nodes and stable `k` prefixes; +replica ranks may shift when membership changes. + +`VirtualPermutation::new(n, seed)` supports `1..=u64::MAX - 1` real nodes, +uses four `u64` fields and no heap state, and offers sequential iteration. +The extra internal label is reserved for the sentinel. `nth(r)` replays +`r + 1` successors; there is no direct rank lookup. Expected `O(k)` prefix +enumeration follows under ideal independent uniform permutations with +constant-cost forward/inverse primitives, **not** as a worst-case guarantee. +The noncryptographic, 64-bit seeded Feistel family approximates that model; +exact uniformity and independence are not claimed. This implementation +replaces the earlier experimental slot-based variants and changes their +output mappings; the existing `ConsistentPermutation` mapping is unchanged. + +See the [algorithm, proof assumptions and API guide](docs/virtual-permutation.md) +and the [reproducible comparison with the existing algorithm](docs/virtual-permutation-performance.md). +The comparison reports fresh-key performance regressions as well as fresh +primary and held-out randomness diagnostics; neither algorithm is replaced +in existing consumers. + ## Applications beyond replication The `ConsistentChooseK` iterator produces a per-key ranking of all `n` nodes in priority order — consistently and with zero memory overhead. This ranking is a strict superset of simple replication and enables drop-in replacements for several well-known algorithms that traditionally require maintaining expensive data structures such as hash rings. diff --git a/crates/consistent-choose-k/benchmarks/Cargo.toml b/crates/consistent-choose-k/benchmarks/Cargo.toml index 5f671fb..a206bce 100644 --- a/crates/consistent-choose-k/benchmarks/Cargo.toml +++ b/crates/consistent-choose-k/benchmarks/Cargo.toml @@ -9,6 +9,12 @@ path = "performance.rs" harness = false test = false +[[bench]] +name = "replica_comparison" +path = "replica_comparison.rs" +harness = false +test = false + [dependencies] consistent-choose-k = { path = "../" } diff --git a/crates/consistent-choose-k/benchmarks/replica_comparison.rs b/crates/consistent-choose-k/benchmarks/replica_comparison.rs new file mode 100644 index 0000000..c39d9df --- /dev/null +++ b/crates/consistent-choose-k/benchmarks/replica_comparison.rs @@ -0,0 +1,238 @@ +//! Paired workloads only in the existing iterator's supported domain. +//! See ../docs/virtual-permutation-performance.md for methodology and results. + +use std::{ + hash::{DefaultHasher, Hash, Hasher}, + hint::black_box, + time::Duration, +}; + +use consistent_choose_k::{ConsistentPermutation, VirtualPermutation}; +use criterion::{ + criterion_group, criterion_main, BatchSize, BenchmarkId, Criterion, SamplingMode, Throughput, +}; +use rand::{rngs::StdRng, RngExt, SeedableRng}; + +const WORKLOAD_SEED: u64 = 0x7065_726d_7574_6531; +const KEY_COUNT: usize = 128; +const NODES: &[u32] = &[ + 1, + 2, + 3, + 4, + 6, + 7, + 8, + 9, + 14, + 15, + 16, + 17, + 30, + 31, + 32, + 33, + 254, + 255, + 256, + 257, + 1000, + 1022, + 1023, + 1024, + 1025, + 65534, + 65535, + 65536, + 65537, + 1_000_000, + (1 << 30) - 2, + (1 << 30) - 1, + 1 << 30, +]; + +fn hash_key(key: u64) -> u64 { + let mut hasher = DefaultHasher::new(); + key.hash(&mut hasher); + hasher.finish() +} + +fn keys() -> Vec { + StdRng::seed_from_u64(WORKLOAD_SEED) + .random_iter() + .take(KEY_COUNT) + .collect() +} + +fn counts(n: u32) -> Vec { + let mut counts = vec![1, 2, 3, 8, 16]; + if n <= 1024 { + counts.extend([n as usize / 4, n as usize]); + } + counts.retain(|&k| k > 0 && k <= n as usize); + counts.sort_unstable(); + counts.dedup(); + counts +} + +fn consume(iter: impl Iterator>, k: usize) { + let sum = iter + .take(k) + .fold(0u64, |sum, node| sum.wrapping_add(node.into())); + black_box(sum); +} + +fn end_to_end(c: &mut Criterion) { + let keys = keys(); + let seeds: Vec<_> = keys.iter().copied().map(hash_key).collect(); + for mode in ["fresh", "seeded"] { + let mut group = c.benchmark_group(format!("sentinel_replicas/{mode}")); + // Both execute exactly one complete query per key, including + // construction, streaming consumption and state destruction. + group.throughput(Throughput::Elements(KEY_COUNT as u64)); + group.sampling_mode(SamplingMode::Flat); + for &n in NODES { + for k in counts(n) { + let input = if mode == "fresh" { &keys } else { &seeds }; + group.bench_function(BenchmarkId::new("layered", format!("n{n}_k{k}")), |b| { + b.iter(|| { + for &key in black_box(input) { + let seed = if mode == "fresh" { hash_key(key) } else { key }; + consume(ConsistentPermutation::new(black_box(n), seed), black_box(k)); + } + }) + }); + group.bench_function(BenchmarkId::new("sentinel", format!("n{n}_k{k}")), |b| { + b.iter(|| { + for &key in black_box(input) { + let seed = if mode == "fresh" { hash_key(key) } else { key }; + consume( + VirtualPermutation::new(u64::from(black_box(n)), seed), + black_box(k), + ); + } + }) + }); + } + } + group.finish(); + } +} + +fn cost_components(c: &mut Criterion) { + let keys = keys(); + let seeds: Vec<_> = keys.iter().copied().map(hash_key).collect(); + let mut setup = c.benchmark_group("sentinel_replicas/setup"); + setup.throughput(Throughput::Elements(KEY_COUNT as u64)); + setup.sampling_mode(SamplingMode::Flat); + setup.bench_function("hash_u64", |b| { + b.iter(|| { + for &key in black_box(&keys) { + black_box(hash_key(key)); + } + }) + }); + for &n in &[17, 257, 1000, 65537, 1 << 30] { + setup.bench_function(BenchmarkId::new("layered", n), |b| { + b.iter(|| { + for &seed in black_box(&seeds) { + black_box(ConsistentPermutation::new(black_box(n), seed)); + } + }) + }); + setup.bench_function(BenchmarkId::new("sentinel", n), |b| { + b.iter(|| { + for &seed in black_box(&seeds) { + black_box(VirtualPermutation::new(u64::from(black_box(n)), seed)); + } + }) + }); + } + setup.finish(); + + for mode in ["stream_only", "collect", "rank_replay"] { + let mut group = c.benchmark_group(format!("sentinel_replicas/{mode}")); + group.throughput(Throughput::Elements(KEY_COUNT as u64)); + group.sampling_mode(SamplingMode::Flat); + for &n in &[17, 257, 1000, 65537, 1 << 30] { + for k in counts(n) { + group.bench_function(BenchmarkId::new("layered", format!("n{n}_k{k}")), |b| { + match mode { + "stream_only" => b.iter_batched_ref( + || { + seeds + .iter() + .map(|&seed| ConsistentPermutation::new(n, seed)) + .collect::>() + }, + |iterators| { + for iter in black_box(iterators) { + consume(iter, black_box(k)); + } + }, + BatchSize::SmallInput, + ), + _ => b.iter(|| { + for &seed in black_box(&seeds) { + let mut iter = ConsistentPermutation::new(black_box(n), seed); + if mode == "collect" { + // Same output width and allocation policy for both. + let mut out = Vec::with_capacity(black_box(k)); + out.extend(iter.take(k).map(u64::from)); + black_box(out); + } else { + // Both methods replay the prefix to answer a rank query. + black_box(iter.nth(black_box(k - 1))); + } + } + }), + } + }); + group.bench_function(BenchmarkId::new("sentinel", format!("n{n}_k{k}")), |b| { + match mode { + "stream_only" => b.iter_batched_ref( + || { + seeds + .iter() + .map(|&seed| VirtualPermutation::new(u64::from(n), seed)) + .collect::>() + }, + |iterators| { + for iter in black_box(iterators) { + consume(iter, black_box(k)); + } + }, + BatchSize::SmallInput, + ), + _ => b.iter(|| { + for &seed in black_box(&seeds) { + let mut iter = + VirtualPermutation::new(u64::from(black_box(n)), seed); + if mode == "collect" { + let mut out = Vec::with_capacity(black_box(k)); + out.extend(iter.take(k)); + black_box(out); + } else { + black_box(iter.nth(black_box(k - 1))); + } + } + }), + } + }); + } + } + group.finish(); + } +} + +criterion_group! { + name = benches; + config = Criterion::default() + .sample_size(20) + .warm_up_time(Duration::from_millis(100)) + .measurement_time(Duration::from_millis(300)) + .nresamples(1000) + .without_plots(); + targets = end_to_end, cost_components +} +criterion_main!(benches); diff --git a/crates/consistent-choose-k/benchmarks/summarize_replica_comparison.py b/crates/consistent-choose-k/benchmarks/summarize_replica_comparison.py new file mode 100644 index 0000000..e43de38 --- /dev/null +++ b/crates/consistent-choose-k/benchmarks/summarize_replica_comparison.py @@ -0,0 +1,58 @@ +#!/usr/bin/env python3 +"""Convert Criterion's batch estimates to ns/query and ns/replica CSV. + +Usage: python3 crates/consistent-choose-k/benchmarks/summarize_replica_comparison.py \ + target/criterion > comparison.csv +""" + +import csv +import json +from pathlib import Path +import re +import sys + + +def summarize(root): + rows = [] + for metadata in root.glob("**/new/benchmark.json"): + benchmark = json.loads(metadata.read_text()) + group = benchmark["group_id"] + if not group.startswith("sentinel_replicas/"): + continue + mode = group.removeprefix("sentinel_replicas/") + algorithm = benchmark["function_id"] + value = benchmark.get("value_str") or "" + match = re.fullmatch(r"n(\d+)_k(\d+)", value) + n, k = map(int, match.groups()) if match else (int(value or 0), 0) + estimates = json.loads((metadata.parent / "estimates.json").read_text()) + mean = estimates["mean"] + interval = mean["confidence_interval"] + divisor = benchmark["throughput"]["Elements"] + rows.append( + ( + mode, + algorithm, + n, + k, + mean["point_estimate"] / divisor, + interval["lower_bound"] / divisor, + interval["upper_bound"] / divisor, + estimates["std_dev"]["point_estimate"] / divisor, + mean["point_estimate"] / divisor / (k if k and mode != "rank_replay" else 1), + ) + ) + if not rows: + raise SystemExit(f"No replica_comparison results found in {root}") + writer = csv.writer(sys.stdout) + writer.writerow( + ["mode", "algorithm", "n", "k", "ns_query", "ci95_low", "ci95_high", + "stddev_ns", "ns_replica"] + ) + for row in sorted(rows): + writer.writerow([*row[:4], *(f"{value:.3f}" for value in row[4:])]) + + +if __name__ == "__main__": + if len(sys.argv) != 2: + raise SystemExit(__doc__) + summarize(Path(sys.argv[1])) diff --git a/crates/consistent-choose-k/docs/permutation-design.md b/crates/consistent-choose-k/docs/permutation-design.md index 6ff55dd..79a9360 100644 --- a/crates/consistent-choose-k/docs/permutation-design.md +++ b/crates/consistent-choose-k/docs/permutation-design.md @@ -4,6 +4,15 @@ This document explains the design of [`ConsistentPermutation`], the per-layer Feistel permutation iterator that this crate uses to drive its `n`-consistent ranking. +This is **survivor-list consistency**. The experimental +[`VirtualPermutation`](virtual-permutation.md) now supplies the same membership +contract through a sentinel-rooted single cycle, with different mappings, +primitive costs and state requirements. See their +[paired performance and randomness comparison](virtual-permutation-performance.md). +Uniformity arguments below model the per-layer bijections as independent +uniform permutations; the actual finite-key, noncryptographic Feistel family +is a practical approximation, not an exact uniform sample from all permutations. + Given a 64-bit `key` and a universe size `n`, the iterator produces a uniformly distributed permutation of `[0, n)` as a streaming iterator satisfying: @@ -450,4 +459,3 @@ Observations matching the theory: pre-build option is `O(k log k)`), which is why the speedup ratio widens with `k`: at `k = 1 000, n = 1 000` the permutation implementation is ~50× faster. - diff --git a/crates/consistent-choose-k/docs/virtual-permutation-performance.md b/crates/consistent-choose-k/docs/virtual-permutation-performance.md new file mode 100644 index 0000000..acb7072 --- /dev/null +++ b/crates/consistent-choose-k/docs/virtual-permutation-performance.md @@ -0,0 +1,314 @@ +# Sentinel-rooted permutation comparison + +The revised `VirtualPermutation` preserves complete survivor-list order, +like the unchanged `ConsistentPermutation`, but is **substantially slower** +in this implementation. It removes heap state, not computational work. +Fresh-key `n=1000,k=3` costs 351.25 versus 36.48 ns/query (9.63x slower); +full enumeration costs 195,818.68 versus 10,761.56 ns (18.20x slower). +No comparably large, repeatable distribution deviations appeared for the +new construction in the primary/held-out diagnostic matrix; that is not +a proof of uniformity or independence. + +All results below were measured anew on the **single-cycle plus sentinel** +construction. They replace the earlier slot-based and matched-network reports. +There is no constant-time output-rank lookup or direct-slot speedup: +both iterators replay a prefix for `nth`. See the +[design, bounds and ideal-model derivation](virtual-permutation.md). +Runner names are `layered` (existing) and `sentinel` (new). + +## Reproduce + +From the repository root, run sequentially, without other CPU-heavy work: + +```sh +cargo test -p consistent-choose-k +cargo test --release -p consistent-choose-k +cargo bench -p consistent-choose-k-benchmarks --bench replica_comparison -- --noplot +python3 crates/consistent-choose-k/benchmarks/summarize_replica_comparison.py \ + target/criterion > comparison.csv +cargo run --release -p consistent-choose-k --example permutation_diagnostics \ + > diagnostics.csv +cargo run --release -p consistent-choose-k --example permutation_diagnostics \ + -- --held-out > held-out.csv +cargo test --release -p consistent-choose-k operation_count_diagnostics \ + -- --ignored --nocapture +``` + +A Criterion filter can select a bounded rerun, for example +`-- 'sentinel_replicas/fresh/.*/n1000_k3$' --noplot`. +The new `sentinel_replicas/` result namespace excludes stale measurements of +the replaced algorithms. Raw samples/estimates remain under `target/criterion`. +The CSV script reports mean ns/query, bootstrap 95% confidence limits, sample +standard deviation and ns/output. Criterion's console times are **batches of +128 queries**: divide by 128 for ns/query, then by `k` for ns/output. +`rank_replay` returns only one output, so its ns/output equals ns/query. + +### Environment and limitations + +Measured September 30, 2026, on Apple M4 Max, native `aarch64-apple-darwin`, +macOS 27.0 build 26A428; `rustc 1.92.0 (ded5c06cf 2025-12-08)`, +LLVM 21.1.3; Apple clang 21.0.0 (`clang-2100.1.1.101`). +The repository bench profile is optimized with debug information and its +configured `-C target-feature=+neon`; no additional LTO, PGO or native-CPU +flags. Criterion 0.8.2 and rand 0.10.3 were resolved locally. + +Each case uses 20 flat samples, 100 ms warmup, 300 ms target measurement, +1,000 bootstrap resamples and no plots. Flat sampling bounds expensive +full-prefix cases; actual durations can exceed the target. All successful +local Cargo commands used +`DEVELOPER_DIR=/Library/Developer/CommandLineTools` to select the independently +installed CLT instead of the default Xcode whose license was unaccepted. +No license was accepted or system setting changed. + +This is a shared host without CPU pinning, frequency control or isolation. +Algorithms run sequentially, existing first in each pair. Intervals concern +repeated timing samples of **one fixed key corpus**, not uncertainty across +all keys, machines or compilers. Some cases are noisy; the wide interval at +`n=257,k=3` is retained rather than discarded. Small differences are not +portable wins. + +### Workload and accounting + +The fixture contains 128 `u64` keys from `StdRng`, seed +`0x7065726d75746531`. Fresh queries hash with `DefaultHasher` inside the timed +region, equally for both algorithms. This is repeated fresh **setup** on a +fixed corpus, not unpredictable new keys every iteration. `StdRng` and +`DefaultHasher` are not cross-version mapping contracts: use the recorded +compiler/dependency versions for the exact corpus. + +| Mode | Timed work | +| --- | --- | +| `fresh` (primary) | Hash key, construct, stream/checksum `k` nodes, destroy | +| `seeded` | Same, but with prehashed seeds; no algorithm-specific cache | +| `setup` | Hash alone, or construct/drop an iterator from a seed | +| `stream_only` | Consume separately prepared iterators with `iter_batched_ref`; construction/destruction excluded for both | +| `collect` | Prehashed construction plus allocate/fill/drop the same-capacity `Vec` | +| `rank_replay` | Prehashed construction and `.nth(k-1)`; both replay `k` successors to return one result | + +All width/key preparation, forward/inverse work, and baseline counter +allocation are charged in the primary comparison. Inputs use `black_box` and +outputs are consumed. Collection uses `u64` elements for both, despite the +existing API's `u32` output. Streaming-only is a component experiment, not +the primary comparison; its prepared-state cache footprints differ. + +The full paired matrix is: + +```text +n = 1,2,3,4,6,7,8,9,14,15,16,17,30,31,32,33, + 254,255,256,257,1000,1022,1023,1024,1025, + 65534,65535,65536,65537,1000000,2^30-2,2^30-1,2^30 +k = valid values from 1,2,3,8,16; also floor(n/4) and n when n<=1024 +``` + +Zeros and duplicates are removed. This covers powers of two and sentinel +boundaries `n+1=2^b`. There are **177 `(n,k)` pairs** in each fresh/seeded +group. Component groups use `n=17,257,1000,65537,2^30`. The completed run +contains **905 estimates**. All paired timings are in the existing iterator's +supported range; larger new domains are not claimed as performance wins. + +## Fresh-query results + +Mean **ns/query [bootstrap 95% confidence interval]**. Ratio is new/existing: +above one is a regression. + +| n | k | Existing layered | New sentinel | Ratio | +| ---: | ---: | ---: | ---: | ---: | +| 1 | 1 | 58.79 [58.58, 59.01] | 14.80 [14.60, 14.95] | 0.25x | +| 7 | 3 | 103.16 [102.76, 103.49] | 864.91 [845.97, 886.02] | 8.38x | +| 8 | 3 | 93.23 [92.57, 94.10] | 2,046.05 [2,036.71, 2,057.79] | 21.95x | +| 8 | 8 | 219.11 [216.88, 221.17] | 6,167.64 [6,071.71, 6,264.64] | 28.15x | +| 17 | 3 | 112.01 [111.87, 112.13] | 1,536.05 [1,533.94, 1,538.30] | 13.71x | +| 33 | 33 | 636.62 [632.07, 641.34] | 23,294.31 [22,161.71, 24,440.97] | 36.59x | +| 254 | 3 | 37.48 [37.27, 37.65] | 468.70 [456.92, 483.16] | 12.51x | +| 255 | 3 | 38.07 [37.72, 38.40] | 494.86 [483.43, 505.74] | 13.00x | +| 256 | 3 | 38.74 [38.62, 38.86] | 956.42 [929.65, 988.09] | 24.69x | +| 257 | 3 | 73.57 [65.12, 88.32] | 903.16 [895.14, 912.31] | 12.28x | +| 1,000 | 1 | 21.30 [21.25, 21.34] | 103.96 [102.76, 105.34] | 4.88x | +| 1,000 | 2 | 28.00 [27.96, 28.04] | 218.57 [215.73, 221.40] | 7.80x | +| 1,000 | 3 | 36.48 [36.11, 36.84] | 351.25 [349.34, 353.43] | 9.63x | +| 1,000 | 8 | 59.65 [59.41, 59.91] | 1,079.05 [1,040.73, 1,151.32] | 18.09x | +| 1,000 | 16 | 101.34 [101.04, 101.63] | 2,393.84 [2,338.73, 2,453.82] | 23.62x | +| 1,000 | 250 | 2,182.76 [2,164.34, 2,198.22] | 44,957.17 [44,920.88, 44,995.31] | 20.60x | +| 1,000 | 1,000 | 10,761.56 [10,554.57, 10,972.06] | 195,818.68 [191,456.60, 200,639.88] | 18.20x | +| 1,022 | 3 | 35.23 [34.97, 35.48] | 343.55 [341.82, 345.57] | 9.75x | +| 1,023 | 3 | 35.90 [35.64, 36.16] | 347.39 [345.30, 349.88] | 9.68x | +| 1,024 | 3 | 36.13 [35.93, 36.30] | 759.81 [758.32, 761.42] | 21.03x | +| 1,024 | 16 | 106.36 [105.76, 107.03] | 5,685.25 [5,593.80, 5,766.20] | 53.45x | +| 1,025 | 3 | 69.08 [68.58, 69.59] | 759.49 [756.87, 762.46] | 10.99x | +| 65,535 | 3 | 38.91 [35.65, 44.28] | 328.64 [318.92, 338.79] | 8.44x | +| 65,536 | 3 | 39.12 [38.49, 39.68] | 747.60 [742.93, 751.59] | 19.11x | +| 65,537 | 3 | 67.87 [66.62, 69.11] | 749.19 [739.39, 760.31] | 11.04x | +| 1,000,000 | 3 | 35.77 [35.24, 36.34] | 328.07 [320.83, 336.15] | 9.17x | +| 2^30 - 2 | 3 | 38.40 [38.06, 38.75] | 305.69 [302.15, 309.49] | 7.96x | +| 2^30 - 1 | 3 | 40.79 [40.24, 41.33] | 300.83 [299.63, 302.02] | 7.37x | +| 2^30 | 3 | 39.51 [39.32, 39.71] | 727.97 [726.74, 729.28] | 18.42x | + +The new implementation is slower in all **176 nontrivial fresh cases**. +The only win is the degenerate `n=1` constant order. Nontrivial ratios range +from 3.51x to 53.45x. At `n=1000,k=3`, the figures are 12.16 versus +117.09 **ns/output**; at `k=1000`, 10.76 versus 195.82 ns/output. +Sample standard deviations are 0.88/4.88 ns for the three-output query and +495.29/10,391.66 ns for full enumeration. + +Every conjugated-cycle step evaluates both Q and its inverse; normalized +lifts also retrace chains. The new Q uses more expensive mixing and more +rounds than the existing iterator's upper layers. This is not an attribution +solely to odd versus even Feistel widths. The sentinel shifts padding +boundaries: `n=1023` has a full internal 1024-label domain, whereas `n=1024` +requires a nearly half-empty 2048-label domain. Operation counts below expose +that cost. Avoiding counter allocation does not offset these extra operations. + +## Setup, streaming, collection and rank replay + +Mean ns/query [95% interval], all prehashed; `n=1000`. +For rank replay, `k` denotes fetching rank `k-1`, not returning `k` outputs. + +| Mode | k | Existing layered | New sentinel | +| --- | ---: | ---: | ---: | +| Constructor + drop | - | 9.59 [9.46, 9.73] | 1.58 [1.57, 1.59] | +| Seeded stream query | 3 | 31.60 [31.27, 31.96] | 336.88 [335.32, 339.24] | +| Seeded stream query | 1,000 | 9,765.81 [9,702.08, 9,834.64] | 186,574.59 [185,570.04, 187,827.35] | +| Stream only | 3 | 14.67 [14.47, 14.87] | 347.11 [338.09, 357.54] | +| Stream only | 1,000 | 9,595.61 [9,574.26, 9,622.33] | 215,975.98 [208,624.79, 224,187.29] | +| Collect | 3 | 42.41 [42.21, 42.60] | 372.48 [371.57, 373.44] | +| Collect | 1,000 | 9,965.63 [9,906.14, 10,036.13] | 206,992.98 [198,395.13, 216,335.19] | +| Rank replay | 3 | 33.34 [32.45, 34.14] | 363.04 [355.80, 371.01] | +| Rank replay | 1,000 | 9,779.07 [9,746.28, 9,820.64] | 187,699.87 [186,407.45, 189,040.15] | + +Hashing a `u64` alone measured 5.41 [4.84, 5.98] ns. Component means are +not additive identities: optimizer behavior, instruction overlap, prepared +state and host noise differ. In particular, component results do not recover +a hidden direct-rank advantage. + +### State and allocation + +On this target the baseline struct is 40 bytes plus one allocation for +`4 * max(1, ceil(log2(n)/2))` bytes of counters: 20 bytes at `n=1000`, +36 at `n=65537`, 60 at `n=2^30`, excluding allocator metadata. +The sentinel iterator is **32 bytes with no heap allocation**. Its temporary +width parameters and scalars are constant-sized, without recursive stack, +prebuilt ring, permutation table or duplicate set. + +Allocation counts follow source inspection, not a custom allocator in the +timed region; the diagnostic prints actual struct sizes. Collection adds one +capacity-`k` `Vec` allocation (8k bytes) to both methods: two total +allocations for the existing iterator and one for the new iterator. + +## Fresh randomness diagnostics + +The example evaluates **3,760,000 complete orders per method per corpus**: +200,000 keys at each `n<=33`, 50,000 at `62,63,64,65`, and 20,000 at +`126,127,128,129,254,255,256,257`. The smaller sizes are +`1,2,3,4,5,6,7,8,9,14,15,16,17,30,31,32,33`. +They include small/odd-width domains and both real-node and sentinel boundaries. + +Primary keys hash `0x7065726d75746531 XOR i`; a disjoint held-out key corpus +hashes `0x686f6c646f757431 XOR i`, using `DefaultHasher`. +The fixed 24/16/8-round Q schedule was not changed or tuned to either corpus +for this construction. Held-out means separate input data, not mathematical +independence supplied by a deterministic hash. + +Metrics cover ranks `0,1,2,n/2,n-2,n-1` where valid, first/middle adjacent +pairs, first/middle and first/last distant pairs, first-three unordered subsets, +all full orders through `n=7`, consecutive application-key primary pairs, +and primaries for `seed` versus `seed XOR (1<<63)`. Duplicate rank choices +are removed. + +Pairs are exact through `n=65`; larger domains use eight contiguous buckets. +For within-key bucket pairs the expected weight is +`size[a]*(size[b] - (a==b))/(n*(n-1))`, not a uniform 64-cell assumption. +Cross-key pairs use `size[a]*size[b]/n^2`. Repeated exact nodes are impossible +within-key and allowed cross-key. Triple histograms are skipped when expected +cell counts would be below ten (all tested sizes above 33); full-order +histograms stop at seven. The actual minimum expected cell count is 11.834. + +The output contains 694 rows across both algorithms per corpus, including +observations, degrees of freedom, minimum expected count, chi-square, maximum +relative cell deviation and empirical total variation (TV). These are +overlapping exploratory diagnostics, **not p-value CI gates**. Consecutive-key +pairs overlap in their keys; rows and nearby sizes are not independent tests. +Bucketed tests can miss correlations inside a bucket. + +### Representative chi-square results + +Values are primary / held-out. All rows use 200,000 observations except +`n=64,65` (50,000) and `n=256` (20,000). + +| n, metric (zero-based ranks) | df | Existing layered | New sentinel | +| --- | ---: | ---: | ---: | +| 7, full order | 5,039 | 5,191.353 / 5,431.710 | 5,188.278 / 5,135.056 | +| 7, first-three subset | 34 | 77.207 / 81.965 | 37.004 / 45.729 | +| 8, first rank | 7 | 4.790 / 4.656 | 4.560 / 9.774 | +| 8, last rank | 7 | 11.501 / 5.734 | 6.189 / 13.977 | +| 8, ordered ranks (0,1) | 55 | 65.968 / 43.725 | 55.089 / 51.515 | +| 8, ordered ranks (4,5) | 55 | 57.991 / 53.807 | 42.932 / 48.067 | +| 8, ordered ranks (0,7) | 55 | 66.922 / 44.737 | 58.735 / 60.360 | +| 8, first-three subset | 55 | 53.474 / 59.828 | 69.037 / 50.819 | +| 8, related-seed primaries | 63 | 86.673 / 51.715 | 83.164 / 76.575 | +| 9, middle rank 4 | 8 | 253.520 / 167.394 | 4.610 / 0.708 | +| 9, ordered ranks (4,5) | 71 | 387.530 / 305.222 | 52.348 / 76.608 | +| 33, middle rank 16 | 32 | 127.629 / 100.534 | 27.809 / 29.831 | +| 33, first-three subset | 5,455 | 5,553.436 / 5,455.392 | 5,429.912 / 5,576.133 | +| 64, ordered ranks (0,1) | 4,031 | 3,967.352 / 3,973.320 | 3,828.813 / 4,140.406 | +| 65, ordered ranks (0,64) | 4,159 | 4,096.640 / 4,155.046 | 4,042.726 / 4,180.339 | +| 256, ordered ranks (0,128), eight buckets | 63 | 500.392 / 488.941 | 74.926 / 77.442 | + +The unchanged baseline has repeatable deviations, especially middle ranks +and distant pairs. At `n=9,rank=4`, its maximum relative cell deviations are +9.91%/7.94%, versus sentinel's 0.90%/0.25%; empirical TV is +1.10%/0.92% versus 0.20%/0.08%. At `n=256`, the distant bucketed pair has +maximum relative deviations 53.96%/53.64% versus 17.76%/14.75%, and TV +5.71%/5.59% versus 2.36%/2.54%. These observations do not justify changing +the user's mapping in this PR. + +Not every new statistic is small. Sentinel's first rank at `n=4` has +chi-square 2.852/13.863 on df=3, maximum relative deviation 0.56%/1.22%. +Its last rank at `n=30` has 55.466/36.187 on df=29, and consecutive-key +primaries at `n=15` have 297.327/203.974 on df=224. Isolated fluctuations +must be interpreted alongside hundreds of correlated checks, not optimized +away by changing the mixer after viewing results. + +TV and maximum cell error include sampling noise and are not corrected +estimates of true family bias: even the new full-order `n=7` histograms have +TV 6.46%/6.39% with only 39.68 expected observations per bin. Diagnostics +cannot prove independence, exact uniformity, unseen-key bounds or +cryptographic security. The new family remains experimental. + +## Sequential operation counts and tails + +The ignored test instruments the actual evaluator **outside timed code**. +It uses 10,000 deterministic mixed seeds per `(n,k)`, 56 cases. Counts are +for whole sequential prefixes, not independent input slots. +One P or P_inverse operation costs two ordinary Q/Q_inverse evaluations; +each Q uses 8, 16 or 24 rounds according to width (width one is XOR). +Percentiles below are nearest-rank percentiles of total **P + P_inverse** +calls per query. `max step` is the largest individual successor evaluation. + +| n | k | Mean P forward | Mean P inverse | Mean Q calls/output | p99 query calls | Max query calls | Max step | +| ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: | +| 1 | 1 | 1.0000 | 0.0000 | 2.0000 | 1 | 1 | 1 | +| 8 | 1 | 3.1700 | 0.6973 | 7.7346 | 12 | 18 | 18 | +| 8 | 8 | 25.2413 | 11.0429 | 9.0710 | 40 | 40 | 20 | +| 255 | 3 | 5.9065 | 0.8688 | 4.5169 | 14 | 20 | 10 | +| 256 | 3 | 11.8445 | 3.8270 | 10.4477 | 32 | 49 | 35 | +| 256 | 256 | 1,012.0098 | 502.0404 | 11.8285 | 1,522 | 1,523 | 49 | +| 1,023 | 3 | 5.9834 | 0.8494 | 4.5552 | 14 | 20 | 12 | +| 1,024 | 3 | 11.9965 | 3.8663 | 10.5752 | 32 | 47 | 31 | +| 65,535 | 3 | 5.9705 | 0.8476 | 4.5454 | 15 | 23 | 15 | +| 65,536 | 3 | 11.9667 | 3.8438 | 10.5403 | 32 | 48 | 36 | +| 2^30 | 3 | 12.0301 | 3.9108 | 10.6273 | 33 | 50 | 36 | +| 2^63 - 1 | 16 | 32.0460 | 11.6256 | 5.4589 | 60 | 75 | 23 | +| u64::MAX - 1 | 16 | 32.1024 | 11.6484 | 5.4688 | 60 | 72 | 25 | + +Observed mean ordinary-PRP calls/output range from 2 to 11.8285. This is +compatible with, but does not prove, the conservative ideal-model expected +bound below 16 described in the design note. The mean visited levels for +three outputs are 5.9834 at `n=1023` versus 8.9764 at `n=1024`; top padding +also adds forward walks and retracing. + +The deterministic invariant test checks that lower-level invocations equal +the selected lower-label subsequence and inverse calls never exceed forward +calls per level along every tested prefix. A constructed full-cycle oracle +still requires more than 4,096 primitive calls for just two outputs at +internal count 2,049. Thus observed tails, and the expected-prefix bound, +are **not a worst-case or adversarial-latency guarantee**. diff --git a/crates/consistent-choose-k/docs/virtual-permutation.md b/crates/consistent-choose-k/docs/virtual-permutation.md new file mode 100644 index 0000000..b700dc8 --- /dev/null +++ b/crates/consistent-choose-k/docs/virtual-permutation.md @@ -0,0 +1,245 @@ +# VirtualPermutation: sentinel-rooted consistent order + +`VirtualPermutation` is an experimental alternative to the unchanged +`ConsistentPermutation`. Both produce distinct nodes, stable `k` prefixes, +and **complete survivor-list restriction**: deleting real node `n` from the +complete order for `n + 1` real nodes recovers the order for `n` real nodes. +Membership is consecutive IDs `0..n`; additions append IDs and removals remove +a suffix. Arbitrary holes, weights and physical-node remapping are out of scope. + +| Property | Existing `ConsistentPermutation` | New `VirtualPermutation` | +| --- | --- | --- | +| Membership changes | Insert/delete entries without reordering survivors | Same order-restriction contract | +| Replica ranks | May shift after insert/delete | May shift after insert/delete | +| Rank query | Replay the iterator | Replay the iterator | +| Supported real-node count | `1..=2^30`, `u32` | `1..=u64::MAX - 1`, `u64` | +| Iterator state | Per-layer heap-allocated counters | Four `u64` fields; no allocation | +| Construction | Interleaved per-layer Feistel streams | Consistent single-cycle successor traversal from a sentinel | + +This revision **replaces the earlier experimental slot-based variants**. +Those variants did not preserve survivor-list order. The experimental output +mapping has changed, the matched-network variant has been removed, and there +is no longer an absolute constant-time rank API. Existing consumers and the +user's `ConsistentPermutation` implementation/mapping are unchanged. + +## API, bounds and key hashing + +```rust +use consistent_choose_k::VirtualPermutation; +use std::hash::{DefaultHasher, Hash, Hasher}; + +let mut hasher = DefaultHasher::new(); +"object-key".hash(&mut hasher); +let seed = hasher.finish(); +let replicas: Vec = VirtualPermutation::new(1000, seed).take(3).collect(); +let third = VirtualPermutation::new(1000, seed).nth(2); // replays three successors +assert_eq!(third, Some(replicas[2])); +let old: Vec<_> = VirtualPermutation::new(1000, seed).collect(); +let restricted: Vec<_> = VirtualPermutation::new(1001, seed) + .filter(|&node| node != 1000).collect(); +assert_eq!(old, restricted); +``` + +The constructor accepts a well-mixed 64-bit seed. Keep it fixed across all +membership and replica counts. `DefaultHasher` matches the examples and +benchmark convention, but Rust does not promise a stable mapping across +versions. Distributed deployments need a specified, versioned key hash and +identical algorithm versions on all participants. Seed collisions give +identical orders. + +Internal label zero is a permanent sentinel; real node `i` has internal label +`i + 1`. Internal count is `n + 1`, so `new(0, seed)` and +`new(u64::MAX, seed)` panic rather than wrap, following the existing +constructor's assertion convention. `n()` returns the original real-node +count. The cloneable, fused iterator reports its remaining size, safely even +when that count exceeds `usize`. `.take(0)` is empty; `.take(n)` is the complete +order; larger requests stop at exhaustion. There is no separate `k` constructor +argument. `nth(r)` uses ordinary iterator replay from the current position; +answering an uncached absolute rank requires replay from a new iterator. + +## Ordinary permutation and guaranteed single cycle + +For each bit width `b`, let `Q(seed,b)` be an ordinary reversible permutation +of `[0,2^b)`, with domain separation depending only on seed and width. Define +the single-cycle primitive by conjugating modular increment: + +```text +P(x) = Q_inverse((Q(x) + 1) mod 2^b) +P_inverse(x) = Q_inverse((Q(x) - 1) mod 2^b) +``` + +Each is two ordinary permutation evaluations: first `Q`, then `Q_inverse`. +Increment is a full cycle, so its conjugate is a full cycle for **every** Q, +not merely with high probability. If Q is uniform on all permutations of a +domain of size `m`, each full cycle has exactly `m` conjugators, so P is uniform +on the `(m-1)!` full cycles. + +The implemented Q is a noncryptographic alternating-XOR Feistel with the +SplitMix64 finalizer as round mixer. It uses 24 rounds at widths 2--4, +16 at widths 5--7, and 8 at widths 8--64; width one is keyed XOR. This fixed +schedule, inherited from the stronger experimental primitive, has been +rediagnosed in the **new construction**, not assumed adequate from old results. +Unequal halves support odd widths. Width seeds are mixed and round keys use +Weyl offsets. The inverse undoes the same updates in reverse order. + +This finite 64-bit family is only a practical pseudorandom approximation, +not exact independent ideal randomness or a proven secure PRP. It cannot +uniformly represent all orders once there are more orders than seeds. +Feistel families also have structural restrictions, such as permutation +parity restrictions on sufficiently large balanced halves. Neither domain +separation, a large round count nor statistical diagnostics proves independence. +Do not use this as encryption or where adversarial keys require cryptographic +security. The existing iterator's different Feistel schedule is not modified. + +All word operations are safe through width 64: modular steps use wrapping +arithmetic followed by masking, half widths never exceed 32, and the largest +dyadic half boundary is `1 << 63`. A conceptual domain cardinality `2^64` +is never stored in a `u64`. + +## Normalized cycle lift + +This is a derived construction/evaluator, not an established production +implementation or an implementation of a published constant-time algorithm. + +Let `A=[0,h)` and let P be a full cycle on `[0,2h)`. Let R be P's cycle +projection onto A: follow P until reaching the next old label. Starting with +`F_1(0)=0`, define: + +```text +F_(2h) = extend(F_h composed with inverse(R), fixing upper labels) composed with P +``` + +For every old label `a`, its P-chain goes through zero or more upper labels +and ends at `R(a)`. The lift changes that final destination to `F_h(a)`. +There are no upper-only cycles in P. Reconnecting all old-node chains using +the lower full cycle therefore produces exactly one full cycle. Its cycle +projection onto A is F_h. For intermediate counts, delete all inactive labels +from the larger cycle. Projection composes, including across dyadic boundaries. + +**Uniform full cycles under the ideal model.** A full cycle P decomposes into +its projected lower cycle R and one ordered upper-node chain attached to each +old label. Every combination of a lower full cycle and such a chain arrangement +corresponds to exactly one full P. Thus uniform P makes R uniform independently +of the arrangement. The lift replaces R with the independent uniform F_h +without changing the arrangement, giving a uniform full F_(2h). Projection of +a uniform full cycle is uniform: each cycle on `m-1` labels has exactly `m-1` +extensions, inserting the new label after any old label. Induction gives +uniform full cycles at every internal count. + +**Rooted list order.** Begin at sentinel zero and repeatedly follow F: + +```text +cursor = 0 +repeat k times, where 0 <= k <= n: + cursor = next_consistent(seed, n + 1, cursor) + emit cursor - 1 +``` + +The sentinel cannot reappear before all `n` real nodes. A full cycle with a +fixed sentinel corresponds bijectively to a real-node order. Cycle deletion +therefore becomes ordinary list deletion: for example `S->A->C->B->S` can +grow into `S->A->D->C->B->S`, preserving the old order. Each ideal real-node +order, ordered prefix, or subset has the appropriate uniform distribution. +Different keys have independent orders **only under independent ideal +primitives across keys**; nodes within a key are sampled without replacement. + +Evaluating successor inputs `0,1,...,k-1` instead of following the cursor would +not implement this contract. Internal input-label successor consistency is +not the public output-rank API. + +## Constant-space successor evaluator + +```text +next_consistent(seed, count, x): + require 0 <= x < count + b = bit_length(count - 1) + while b > 0: + half = 1 << (b - 1) + y = P(seed, b, x) + if y >= count: + x = y + continue + if y >= half: + return y + while x >= half: + x = P_inverse(seed, b, x) + count = half + b -= 1 + return 0 +``` + +Edges whose outputs are upper labels are unchanged by the lift, so walking +inactive upper labels can use P directly. Upon reaching a lower output, +the backward walk recovers the old input at the start of that chain; recursion +then supplies its correct lower-cycle destination. This is the explicit lift +without constructing any tables. The actual implementation uses loops, not a +recursive stack. Walks terminate because P is a full finite cycle intersecting +the retained/lower set. No retry cap or mapping-changing fallback is used. + +## Expected adaptive-prefix work + +The cursor depends on previous outputs. A fixed-input expected-cost bound +would therefore be insufficient. Instead count work over the entire rooted +prefix, under independent uniform ideal Q at each width and constant-cost +forward/inverse ordinary permutation calls. Fix `n,k` before sampling those +primitives. A zero-length prefix does no traversal; the strict bounds below +are for `1<=k<=n`. + +Let `M` be the smallest power of two at least `n+1`. The lifted full cycle +F_M (not the raw primitive P) is uniform. Reaching the first `k` active real outputs visits, in expectation, +`k*M/(n+1)` top-level successors: sample without replacement from the +`M-1` non-sentinel labels until the `k`th of the `n` active labels. This is +less than `2k`. + +At each lower full dyadic domain of size `L=M/2,M/4,...,2`, the successor +invocations form the rooted lower-label subsequence. Their count equals the +number of selected final internal labels in `1..L-1`, so its expectation is +`k*(L-1)/n`. This uses marginal uniformity of the rooted real-node order, not +independence of adaptive calls. Summing these lower-level expectations gives +less than `k*M/n`, which is at most `2k` for `n>=1`. + +At any level, a backward walk retraces an upper-node chain already traversed +by that level's forward prefix. This includes inactive labels just visited +in the top-level walk. Each upper label is retraced at most once; the prefix +does not wrap back through the sentinel. Consequently inverse calls are +bounded pathwise by forward calls across the prefix. Combining the counts +gives a conservative **expected bound below `8k` P/P_inverse calls, or `16k` +ordinary Q/Q_inverse calls**. Width setup is bounded per level invocation and +does not change expected `O(k)` work. This bounds cumulative work from the +sentinel, not work conditioned on an arbitrary already-observed prefix. + +This is an ideal-model derivation, not a published complexity theorem or a +guarantee for every finite-key seed. An unlucky full cycle can force +linear-in-domain work even for a short prefix; there is **no worst-case +`O(k)` or adversarial-latency guarantee**. The iterator uses `O(1)` auxiliary +state excluding returned outputs. There is no prebuilt ring, per-key table, +permutation array, cached prefix or duplicate set. + +## Verification and references + +Tests cover forward/inverse round trips for both Q and P through all 64 +widths, guaranteed single-cycle coverage on small domains, explicit +table-based lift/projection equivalence, rooted traversal and sentinel return, +complete survivor-list restriction, `k` prefixes, default `nth` replay, fused +exhaustion, and small/large dyadic and sentinel boundaries. They explicitly +reject overflowing sentinel counts. A constructed long-walk case checks that +evaluation is not truncated. + +Exhaustive ideal checks cover all 24 ordinary size-four conjugators and all +30,240 independent full-cycle families at sizes 2,4,8. Every real-node order +at sizes 1 through 7 appears equally often and agrees with the explicit +reference and survivor deletion. These are regression checks, not substitutes +for the ideal-model argument or proofs of security. + +The mathematical antecedents concern virtual permutations and cycle projection: + +- Neretin, [Virtual permutations and polymorphisms, section 1.2](https://arxiv.org/html/2202.12978v1#S1): + cycle deletion and equivariance. +- Bourgade, Najnudel and Nikeghbali, + [A unitary extension of virtual permutations, section 1](https://arxiv.org/html/1102.2633v1#S1): + cycle projections, Chinese restaurant construction and uniform coherent families. + +These sources do **not** supply this evaluator, Feistel schedule, single-cycle +sentinel specialization or performance bound. The +[performance and randomness report](virtual-permutation-performance.md) +records actual final-construction measurements and their limitations. diff --git a/crates/consistent-choose-k/examples/permutation_diagnostics.rs b/crates/consistent-choose-k/examples/permutation_diagnostics.rs new file mode 100644 index 0000000..647f430 --- /dev/null +++ b/crates/consistent-choose-k/examples/permutation_diagnostics.rs @@ -0,0 +1,230 @@ +//! Deterministic distribution diagnostics, not statistical CI gates. +//! Run with no arguments for the primary corpus, or --held-out. + +use std::hash::{DefaultHasher, Hash, Hasher}; + +use consistent_choose_k::{ConsistentPermutation, VirtualPermutation}; + +fn key_seed(key: u64, corpus: u64) -> u64 { + let mut hasher = DefaultHasher::new(); + (corpus ^ key).hash(&mut hasher); + hasher.finish() +} + +struct Histogram { + counts: Vec, + probabilities: Vec, +} + +impl Histogram { + fn uniform(cells: usize) -> Self { + Self { + counts: vec![0; cells], + probabilities: vec![1.0 / cells as f64; cells], + } + } + + fn pairs(n: usize, buckets: usize, distinct: bool) -> Self { + let mut sizes = vec![0; buckets]; + for node in 0..n { + sizes[node * buckets / n] += 1; + } + let mut probabilities = vec![]; + for a in 0..buckets { + for b in 0..buckets { + let available = sizes[b] - usize::from(distinct && a == b); + probabilities + .push((sizes[a] * available) as f64 / (n * (n - usize::from(distinct))) as f64); + } + } + Self { + counts: vec![0; buckets * buckets], + probabilities, + } + } + + fn record(&mut self, cell: usize) { + self.counts[cell] += 1; + } + + fn report(&self, algorithm: &str, n: usize, metric: &str) { + let total: u64 = self.counts.iter().sum(); + let mut chi2 = 0.0; + let mut tv = 0.0; + let mut max_relative: f64 = 0.0; + let mut min_expected = f64::INFINITY; + let mut cells = 0; + for (&count, &probability) in self.counts.iter().zip(&self.probabilities) { + if probability == 0.0 { + assert_eq!(count, 0, "impossible pair observed"); + continue; + } + let expected = total as f64 * probability; + assert!(expected >= 9.99, "inadequate expected cell count"); + let delta = (count as f64 - expected).abs(); + chi2 += delta * delta / expected; + tv += delta / total as f64 / 2.0; + max_relative = max_relative.max(delta / expected); + min_expected = min_expected.min(expected); + cells += 1; + } + println!( + "{algorithm},{n},{metric},{total},{},{min_expected:.3},{chi2:.3},{:.3},{tv:.6}", + cells - 1, + max_relative * 100.0 + ); + } +} + +fn order(algorithm: &str, n: usize, seed: u64) -> Vec { + match algorithm { + "layered" => ConsistentPermutation::new(n as u32, seed) + .map(|v| v as usize) + .collect(), + "sentinel" => VirtualPermutation::new(n as u64, seed) + .map(|v| v as usize) + .collect(), + _ => unreachable!(), + } +} + +fn first(algorithm: &str, n: usize, seed: u64) -> usize { + match algorithm { + "layered" => ConsistentPermutation::new(n as u32, seed) + .next() + .expect("nonempty") as usize, + "sentinel" => VirtualPermutation::new(n as u64, seed) + .next() + .expect("nonempty") as usize, + _ => unreachable!(), + } +} + +fn permutation_index(values: &[usize]) -> usize { + let mut index = 0; + for (i, &value) in values.iter().enumerate() { + index = index * (values.len() - i) + values[i + 1..].iter().filter(|&&v| v < value).count(); + } + index +} + +fn main() { + let mut args = std::env::args().skip(1); + let corpus = match args.next().as_deref() { + None => 0x7065_726d_7574_6531, + Some("--held-out") => 0x686f_6c64_6f75_7431, + Some(_) => panic!("usage: permutation_diagnostics [--held-out]"), + }; + assert!( + args.next().is_none(), + "usage: permutation_diagnostics [--held-out]" + ); + println!("# corpus={corpus:#x}; exact bins require >=10 expected observations"); + println!( + "# iterator_bytes: layered={}, sentinel={}; layered also owns heap counters", + std::mem::size_of::(), + std::mem::size_of::() + ); + println!( + "algorithm,n,metric,observations,df,min_expected,chi2,max_relative_percent,total_variation" + ); + for n in [ + 1usize, 2, 3, 4, 5, 6, 7, 8, 9, 14, 15, 16, 17, 30, 31, 32, 33, 62, 63, 64, 65, 126, 127, + 128, 129, 254, 255, 256, 257, + ] { + let samples = if n <= 33 { + 200_000 + } else if n <= 65 { + 50_000 + } else { + 20_000 + }; + let buckets = if samples / (n * n) >= 10 { n } else { 8 }; + let pair_kind = if buckets == n { "exact" } else { "bucket8" }; + let mut ranks = vec![ + 0, + 1.min(n - 1), + 2.min(n - 1), + n / 2, + n.saturating_sub(2), + n - 1, + ]; + ranks.sort_unstable(); + ranks.dedup(); + let mut rank_pairs = vec![ + (0, 1.min(n - 1)), + (n / 2, (n / 2 + 1).min(n - 1)), + (0, n / 2), + (0, n - 1), + ]; + rank_pairs.retain(|(a, b)| a != b); + rank_pairs.sort_unstable(); + rank_pairs.dedup(); + let triple_cells = if n >= 3 { n * (n - 1) * (n - 2) / 6 } else { 0 }; + let triple_ok = triple_cells > 0 && samples / triple_cells >= 10; + if n >= 3 && !triple_ok { + println!( + "# n={n}: skipping exact choose_3; expected={:.3}", + samples as f64 / triple_cells as f64 + ); + } + for algorithm in ["layered", "sentinel"] { + let mut marginals: Vec<_> = ranks.iter().map(|_| Histogram::uniform(n)).collect(); + let mut pairs: Vec<_> = rank_pairs + .iter() + .map(|_| Histogram::pairs(n, buckets, true)) + .collect(); + let mut triples = triple_ok.then(|| Histogram::uniform(triple_cells)); + let mut full = (n <= 7).then(|| Histogram::uniform((1..=n).product())); + let mut adjacent_keys = Histogram::pairs(n, buckets, false); + let mut related_seeds = Histogram::pairs(n, buckets, false); + let mut previous = None; + for key in 0..samples as u64 { + let seed = key_seed(key, corpus); + let values = order(algorithm, n, seed); + assert_eq!(values.len(), n); + for (hist, &rank) in marginals.iter_mut().zip(&ranks) { + hist.record(values[rank]); + } + for (hist, &(a, b)) in pairs.iter_mut().zip(&rank_pairs) { + assert_ne!(values[a], values[b]); + hist.record((values[a] * buckets / n) * buckets + values[b] * buckets / n); + } + if let Some(hist) = &mut triples { + let mut triple = [values[0], values[1], values[2]]; + triple.sort_unstable(); + let [a, b, c] = triple; + hist.record(c * (c - 1) * (c - 2) / 6 + b * (b - 1) / 2 + a); + } + if let Some(hist) = &mut full { + hist.record(permutation_index(&values)); + } + let primary = values[0] * buckets / n; + if let Some(old) = previous { + adjacent_keys.record(old * buckets + primary); + } + previous = Some(primary); + let related = first(algorithm, n, seed ^ (1u64 << 63)) * buckets / n; + related_seeds.record(primary * buckets + related); + } + for (hist, rank) in marginals.iter().zip(&ranks) { + hist.report(algorithm, n, &format!("rank_{rank}")); + } + for (hist, (a, b)) in pairs.iter().zip(&rank_pairs) { + hist.report(algorithm, n, &format!("ordered_{a}_{b}_{pair_kind}")); + } + if let Some(hist) = triples { + hist.report(algorithm, n, "choose_3"); + } + if let Some(hist) = full { + hist.report(algorithm, n, "full_order"); + } + adjacent_keys.report( + algorithm, + n, + &format!("consecutive_key_primaries_{pair_kind}"), + ); + related_seeds.report(algorithm, n, &format!("seed_bitflip_primaries_{pair_kind}")); + } + } +} diff --git a/crates/consistent-choose-k/src/lib.rs b/crates/consistent-choose-k/src/lib.rs index 25b7e0f..4278f83 100644 --- a/crates/consistent-choose-k/src/lib.rs +++ b/crates/consistent-choose-k/src/lib.rs @@ -3,6 +3,7 @@ mod consistent_hash; mod consistent_permutation; mod consistent_reservoir; mod node_map; +mod virtual_permutation; pub use choose_k::ConsistentChooseKHasher; pub use consistent_hash::{ ConsistentHashIterator, ConsistentHashRevIterator, ConsistentHasher, HashSeqBuilder, @@ -11,3 +12,4 @@ pub use consistent_hash::{ pub use consistent_permutation::ConsistentPermutation; pub use consistent_reservoir::ConsistentReservoir; pub use node_map::ConsistentNodeMap; +pub use virtual_permutation::VirtualPermutation; diff --git a/crates/consistent-choose-k/src/virtual_permutation.rs b/crates/consistent-choose-k/src/virtual_permutation.rs new file mode 100644 index 0000000..e12ef26 --- /dev/null +++ b/crates/consistent-choose-k/src/virtual_permutation.rs @@ -0,0 +1,708 @@ +//! Sentinel-rooted, single-cycle consistent permutations. +//! See `docs/virtual-permutation.md` for the construction and assumptions. + +use std::iter::FusedIterator; + +const WEYL: u64 = 0x9E37_79B9_7F4A_7C15; + +#[inline] +fn mix(mut x: u64) -> u64 { + x = (x ^ (x >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); + x = (x ^ (x >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); + x ^ (x >> 31) +} + +trait Permutation { + fn forward(&self, x: u64) -> u64; + fn inverse(&self, x: u64) -> u64; +} + +/// Ordinary noncryptographic Feistel Q, not itself the consistent map or cycle. +struct WordPermutation { + key: u64, + bits: u32, + right_bits: u32, + left_mask: u64, + right_mask: u64, +} + +impl WordPermutation { + #[inline] + fn new(seed: u64, bits: u32) -> Self { + debug_assert!((1..=64).contains(&bits)); + let left_bits = bits / 2; + let right_bits = bits - left_bits; + Self { + key: mix(seed ^ WEYL.wrapping_mul(u64::from(bits))), + bits, + right_bits, + left_mask: (1u64 << left_bits) - 1, + right_mask: (1u64 << right_bits) - 1, + } + } + + #[inline] + fn round(&self, x: u64, round: u64) -> u64 { + mix(x ^ self.key.wrapping_add(WEYL.wrapping_mul(round + 1))) + } + + #[inline] + fn transform_rounds(&self, x: u64) -> u64 { + let mut left = x >> self.right_bits; + let mut right = x & self.right_mask; + for i in 0..PAIRS { + let pair = if INVERSE { PAIRS - 1 - i } else { i }; + if INVERSE { + right ^= self.round(left, 2 * pair + 1) & self.right_mask; + left ^= self.round(right, 2 * pair) & self.left_mask; + } else { + left ^= self.round(right, 2 * pair) & self.left_mask; + right ^= self.round(left, 2 * pair + 1) & self.right_mask; + } + } + (left << self.right_bits) | right + } + + #[inline] + fn transform(&self, x: u64) -> u64 { + match self.bits { + 1 => x ^ (self.key & 1), + 2..=4 => self.transform_rounds::<12, INVERSE>(x), + 5..=7 => self.transform_rounds::<8, INVERSE>(x), + _ => self.transform_rounds::<4, INVERSE>(x), + } + } +} + +impl Permutation for WordPermutation { + #[inline] + fn forward(&self, x: u64) -> u64 { + self.transform::(x) + } + + #[inline] + fn inverse(&self, x: u64) -> u64 { + self.transform::(x) + } +} + +struct SingleCycle { + ordinary: Q, + mask: u64, +} + +impl Permutation for SingleCycle { + #[inline] + fn forward(&self, x: u64) -> u64 { + self.ordinary + .inverse(self.ordinary.forward(x).wrapping_add(1) & self.mask) + } + + #[inline] + fn inverse(&self, x: u64) -> u64 { + self.ordinary + .inverse(self.ordinary.forward(x).wrapping_sub(1) & self.mask) + } +} + +#[inline] +fn cycle(seed: u64, bits: u32) -> SingleCycle { + SingleCycle { + ordinary: WordPermutation::new(seed, bits), + mask: u64::MAX >> (64 - bits), + } +} + +#[inline] +fn evaluate(mut count: u64, mut x: u64, mut layer: impl FnMut(u32) -> P) -> u64 { + debug_assert!(count > 0 && x < count); + let mut bits = u64::BITS - (count - 1).leading_zeros(); + while bits > 0 { + let half = 1u64 << (bits - 1); + let permutation = layer(bits); + loop { + let y = permutation.forward(x); + if y >= count { + x = y; + continue; + } + if y >= half { + return y; + } + // Recover the old input at the start of this chain, rather than + // its old destination y, before evaluating the lower-level cycle. + while x >= half { + x = permutation.inverse(x); + } + break; + } + count = half; + bits -= 1; + } + 0 +} + +/// An allocation-free, per-key ordering of `0..n` preserving survivor order. +/// +/// Every prefix contains distinct nodes. Appending a node inserts it somewhere +/// in the complete order; removing the last node deletes it without reordering +/// survivors. Membership must be consecutive IDs `0..n`. Replica ranks may +/// shift when membership changes. +/// +/// The iterator traverses one consistent cycle from a permanent internal +/// sentinel. It does not evaluate independent rank inputs or cache a prefix. +/// `nth(r)` replays `r + 1` successors from the current position; there is no +/// constant-time absolute-rank API. +/// +/// Under independent uniform ideal permutations per key and bit width, the +/// rooted order is uniform and sequential work is expected `O(k)` for `k` +/// outputs with constant-cost primitives. This is not a worst-case bound. +/// The finite 64-bit seeded, 8/16/24-round Feistel family is noncryptographic; +/// exact uniformity, independence and cryptographic security are not claimed. +/// State is four `u64` fields, with no ring, table, duplicate set or allocation. +/// +/// ``` +/// use consistent_choose_k::VirtualPermutation; +/// +/// let seed = 0x1234_5678_9abc_def0; // normally a well-mixed hash of the key +/// let old: Vec<_> = VirtualPermutation::new(100, seed).collect(); +/// let new: Vec<_> = VirtualPermutation::new(101, seed) +/// .filter(|&node| node != 100).collect(); +/// assert_eq!(old, new); +/// assert_eq!(VirtualPermutation::new(100, seed).take(3).collect::>(), old[..3]); +/// ``` +#[derive(Clone, Debug)] +pub struct VirtualPermutation { + seed: u64, + count: u64, + cursor: u64, + remaining: u64, +} + +impl VirtualPermutation { + /// Construct an ordering of `1..=u64::MAX - 1` real nodes. + /// + /// Supply a well-mixed key hash, held fixed across membership and replica + /// counts. Internal label zero is the sentinel; real node `i` has label + /// `i + 1`. The extra sentinel requires representable internal count `n + 1`. + /// + /// # Panics + /// + /// Panics if `n == 0` or `n == u64::MAX`. + pub fn new(n: u64, seed: u64) -> Self { + assert!(n > 0, "n must be at least 1"); + assert!(n < u64::MAX, "n must be at most u64::MAX - 1"); + Self { + seed, + count: n + 1, + cursor: 0, + remaining: n, + } + } + + /// Number of real nodes, independent of iterator position. + pub fn n(&self) -> u64 { + self.count - 1 + } +} + +impl Iterator for VirtualPermutation { + type Item = u64; + + fn next(&mut self) -> Option { + if self.remaining == 0 { + return None; + } + self.cursor = evaluate(self.count, self.cursor, |bits| cycle(self.seed, bits)); + debug_assert!(self.cursor > 0 && self.cursor < self.count); + self.remaining -= 1; + Some(self.cursor - 1) + } + + fn size_hint(&self) -> (usize, Option) { + match usize::try_from(self.remaining) { + Ok(remaining) => (remaining, Some(remaining)), + Err(_) => (usize::MAX, None), + } + } +} + +impl FusedIterator for VirtualPermutation {} + +#[cfg(test)] +mod tests { + use super::*; + use std::cell::Cell; + + struct CountedCycle<'a> { + cycle: SingleCycle, + calls: &'a Cell<[u64; 3]>, + } + + impl Permutation for CountedCycle<'_> { + fn forward(&self, x: u64) -> u64 { + let mut v = self.calls.get(); + v[0] += 1; + self.calls.set(v); + self.cycle.forward(x) + } + fn inverse(&self, x: u64) -> u64 { + let mut v = self.calls.get(); + v[1] += 1; + self.calls.set(v); + self.cycle.inverse(x) + } + } + + fn seed(i: u64) -> u64 { + mix(i.wrapping_add(WEYL)) + } + + fn successor(count: u64, key: u64, x: u64) -> u64 { + evaluate(count, x, |bits| cycle(key, bits)) + } + + #[test] + fn ordinary_and_cycle_round_trips_all_widths() { + for bits in 1..=64 { + let mask = u64::MAX >> (64 - bits); + for key in [0, 1, u64::MAX, seed(42)] { + let q = WordPermutation::new(key, bits); + let p = cycle(key, bits); + for x in [0, 1, mask / 2, mask / 2 + 1, mask] + .into_iter() + .chain((0..64).map(|i| seed(i) & mask)) + { + assert!(q.forward(x) <= mask && p.forward(x) <= mask); + assert_eq!(q.inverse(q.forward(x)), x); + assert_eq!(q.forward(q.inverse(x)), x); + assert_eq!(p.inverse(p.forward(x)), x); + assert_eq!(p.forward(p.inverse(x)), x); + } + } + } + } + + #[test] + fn conjugates_are_single_cycles() { + for bits in 1..=9 { + for key in 0..16 { + let p = cycle(seed(key), bits); + let mut seen = vec![false; 1 << bits]; + let mut x = 0; + for _ in 0..1 << bits { + assert!(!seen[x as usize]); + seen[x as usize] = true; + assert_eq!(p.inverse(p.forward(x)), x); + x = p.forward(x); + } + assert_eq!(x, 0); + assert!(seen.into_iter().all(|v| v)); + } + } + } + + fn project(p: &[u64], count: usize) -> Vec { + (0..count) + .map(|x| { + let mut y = p[x]; + while y >= count as u64 { + y = p[y as usize]; + } + y + }) + .collect() + } + + fn lift(lower: &[u64], p: &[u64]) -> Vec { + let h = lower.len(); + let r = project(p, h); + let mut inverse = vec![0; h]; + for (x, &y) in r.iter().enumerate() { + inverse[y as usize] = x; + } + p.iter() + .map(|&y| { + if y < h as u64 { + lower[inverse[y as usize]] + } else { + y + } + }) + .collect() + } + + fn rooted_order(p: &[u64]) -> Vec { + let mut x = 0; + let mut order = vec![]; + for _ in 1..p.len() { + x = p[x as usize]; + assert_ne!(x, 0, "sentinel reached early"); + order.push(x - 1); + } + assert_eq!(p[x as usize], 0, "cycle does not close at sentinel"); + order + } + + #[test] + fn complete_orders_prefixes_and_survivor_restriction() { + for key in 0..16 { + let key = seed(key); + let mut previous_order = vec![]; + let mut previous_map = vec![0]; + for n in 1..=257 { + let iter = VirtualPermutation::new(n, key); + let order: Vec<_> = iter.clone().collect(); + let mut sorted = order.clone(); + sorted.sort_unstable(); + assert_eq!(sorted, (0..n).collect::>()); + assert_eq!( + order + .iter() + .copied() + .filter(|&node| node < n - 1) + .collect::>(), + previous_order + ); + for k in [0, 1, 2, 3, 8, n / 2, n] { + if k <= n { + assert_eq!( + iter.clone().take(k as usize).collect::>(), + order[..k as usize] + ); + } + } + for rank in [0, n / 2, n - 1] { + let mut replay = iter.clone(); + assert_eq!(replay.nth(rank as usize), Some(order[rank as usize])); + assert_eq!(replay.next(), order.get(rank as usize + 1).copied()); + } + let map: Vec<_> = (0..=n).map(|x| successor(n + 1, key, x)).collect(); + assert_eq!(rooted_order(&map), order); + assert_eq!(project(&map, n as usize), previous_map); + previous_map = map; + previous_order = order; + } + } + } + + #[test] + fn matches_explicit_single_cycle_lifts() { + for key in 0..8 { + let key = seed(key); + let mut full = vec![0]; + for bits in 1..=8 { + let p = cycle(key, bits); + let table: Vec<_> = (0..1 << bits).map(|x| p.forward(x)).collect(); + full = lift(&full, &table); + for count in full.len() / 2 + 1..=full.len() { + let projected = project(&full, count); + assert_eq!( + (0..count as u64) + .map(|x| successor(count as u64, key, x)) + .collect::>(), + projected + ); + assert_eq!( + VirtualPermutation::new(count as u64 - 1, key).collect::>(), + rooted_order(&projected) + ); + } + } + } + } + + #[test] + fn large_boundaries_and_sentinel_overflow() { + let mut sizes = vec![u64::MAX - 2, u64::MAX - 1]; + for bits in 1..64 { + let power = 1u64 << bits; + sizes.extend([power.saturating_sub(2).max(1), power - 1, power, power + 1]); + } + for n in sizes { + for key in [0, 1, u64::MAX, seed(42)] { + let k = n.min(16) as usize; + let values: Vec<_> = VirtualPermutation::new(n, key).take(k).collect(); + assert!(values.iter().all(|&v| v < n)); + let mut distinct = values.clone(); + distinct.sort_unstable(); + distinct.dedup(); + assert_eq!(distinct.len(), k); + if n < u64::MAX - 1 { + let restricted: Vec<_> = VirtualPermutation::new(n + 1, key) + .take(k + 1) + .filter(|&node| node < n) + .take(k) + .collect(); + assert_eq!(values, restricted); + } + } + } + } + + #[test] + fn sequential_level_subsequences_and_inverse_charging() { + for key in 0..8 { + for n in 1u64..=129 { + let bits = u64::BITS - n.leading_zeros(); + let calls: Vec<_> = (0..=bits).map(|_| Cell::new([0; 3])).collect(); + let mut selected = vec![0; bits as usize + 1]; + let mut cursor = 0; + for k in 1..=n { + cursor = evaluate(n + 1, cursor, |width| CountedCycle { + cycle: cycle(seed(key), width), + calls: &calls[width as usize], + }); + for width in 1..=bits { + let [forward, inverse, _] = calls[width as usize].get(); + assert!( + inverse <= forward, + "each upper chain is retraced at most once" + ); + if width < bits { + selected[width as usize] += u64::from(cursor < 1 << width); + assert_eq!(forward, selected[width as usize]); + } else { + assert!(forward >= k && forward < 1 << bits); + } + } + } + } + } + } + + #[test] + fn iterator_boundaries() { + let mut p = VirtualPermutation::new(1, 0); + assert_eq!(p.n(), 1); + assert_eq!(p.size_hint(), (1, Some(1))); + assert_eq!(p.clone().take(0).count(), 0); + assert_eq!(p.next(), Some(0)); + assert_eq!(successor(p.count, p.seed, p.cursor), 0); + assert_eq!(p.size_hint(), (0, Some(0))); + assert_eq!(p.next(), None); + assert_eq!(p.next(), None); + assert_eq!(p.nth(usize::MAX), None); + let p = VirtualPermutation::new(u64::MAX - 1, 0); + let expected = usize::try_from(u64::MAX - 1).ok(); + assert_eq!(p.size_hint(), (expected.unwrap_or(usize::MAX), expected)); + } + + #[test] + #[should_panic(expected = "n must be at least 1")] + fn invalid_empty_domain() { + VirtualPermutation::new(0, 0); + } + + #[test] + #[should_panic(expected = "n must be at most u64::MAX - 1")] + fn invalid_sentinel_overflow() { + VirtualPermutation::new(u64::MAX, 0); + } + + fn permutations(mut values: Vec) -> Vec> { + fn visit(values: &mut [u64], start: usize, out: &mut Vec>) { + if start == values.len() { + out.push(values.to_vec()); + } else { + for i in start..values.len() { + values.swap(start, i); + visit(values, start + 1, out); + values.swap(start, i); + } + } + } + let mut out = vec![]; + visit(&mut values, 0, &mut out); + out + } + + fn all_cycles(count: usize) -> Vec> { + permutations((1..count as u64).collect()) + .into_iter() + .map(|order| { + let mut map = vec![0; count]; + let mut x = 0; + for y in order { + map[x] = y; + x = y as usize; + } + map + }) + .collect() + } + + struct Table<'a>(&'a [u64]); + + impl Permutation for Table<'_> { + fn forward(&self, x: u64) -> u64 { + self.0[x as usize] + } + fn inverse(&self, x: u64) -> u64 { + self.0.iter().position(|&y| y == x).expect("bijection") as u64 + } + } + + #[test] + fn exhaustive_ideal_conjugators_and_rooted_orders() { + use std::collections::BTreeMap; + + let mut conjugates = BTreeMap::new(); + for q in permutations((0..4).collect()) { + let p = SingleCycle { + ordinary: Table(&q), + mask: 3, + }; + let table: Vec<_> = (0..4).map(|x| p.forward(x)).collect(); + *conjugates.entry(table).or_insert(0) += 1; + } + assert_eq!(conjugates.len(), 6); + assert!(conjugates.values().all(|&v| v == 4)); + + let p2 = [1, 0]; + let cycles8 = all_cycles(8); + let mut counts: Vec, usize>> = (1..=7).map(|_| BTreeMap::new()).collect(); + for p4 in all_cycles(4) { + let f4 = lift(&p2, &p4); + for p8 in &cycles8 { + let f8 = lift(&f4, p8); + let mut previous = vec![]; + for n in 1..=7 { + let order = rooted_order(&project(&f8, n + 1)); + let mut actual = vec![]; + let mut cursor = 0; + for _ in 0..n { + cursor = evaluate(n as u64 + 1, cursor, |bits| { + Table(match bits { + 1 => &p2, + 2 => &p4, + 3 => p8, + _ => unreachable!(), + }) + }); + assert_ne!(cursor, 0); + actual.push(cursor - 1); + } + assert_eq!(actual, order); + assert_eq!( + order + .iter() + .copied() + .filter(|&v| v < n as u64 - 1) + .collect::>(), + previous + ); + previous = order.clone(); + *counts[n - 1].entry(order).or_insert(0) += 1; + } + } + } + for (i, counts) in counts.iter().enumerate() { + let factorial: usize = (1..=i + 1).product(); + assert_eq!(counts.len(), factorial); + assert!(counts.values().all(|&v| v == 30240 / factorial)); + } + } + + #[test] + fn unlucky_sequential_walks_are_not_truncated() { + struct ReverseCycle<'a> { + mask: u64, + calls: &'a Cell, + } + impl Permutation for ReverseCycle<'_> { + fn forward(&self, x: u64) -> u64 { + self.calls.set(self.calls.get() + 1); + x.wrapping_sub(1) & self.mask + } + fn inverse(&self, x: u64) -> u64 { + self.calls.set(self.calls.get() + 1); + x.wrapping_add(1) & self.mask + } + } + let calls = Cell::new(0); + let factory = |bits| ReverseCycle { + mask: (1u64 << bits) - 1, + calls: &calls, + }; + let first = evaluate(2049, 0, factory); + assert_eq!(first, 2048); + assert_eq!(evaluate(2049, first, factory), 2047); + assert!(calls.get() > 4096); + } + + #[test] + #[ignore = "deterministic sequential operation-count report, not a timing CI gate"] + fn operation_count_diagnostics() { + println!("n,k,mean_forward,mean_inverse,mean_q_per_output,mean_levels,p50_calls,p99_calls,max_calls,max_step_calls"); + for n in [ + 1, + 2, + 3, + 7, + 8, + 9, + 255, + 256, + 257, + 1023, + 1024, + 1025, + 65535, + 65536, + 65537, + 1 << 30, + (1 << 63) - 1, + u64::MAX - 1, + ] { + let mut ks = vec![1, 3, 16]; + if n <= 256 { + ks.extend([n / 4, n]); + } + ks.retain(|&k| k > 0 && k <= n); + ks.sort_unstable(); + ks.dedup(); + for k in ks { + let calls = Cell::new([0; 3]); + let mut totals = [0u64; 3]; + let mut samples = vec![]; + let mut max_step = 0; + for key in 0..10_000 { + calls.set([0; 3]); + let mut cursor = 0; + for _ in 0..k { + let before = calls.get()[0] + calls.get()[1]; + cursor = evaluate(n + 1, cursor, |bits| { + let mut v = calls.get(); + v[2] += 1; + calls.set(v); + CountedCycle { + cycle: cycle(seed(key), bits), + calls: &calls, + } + }); + assert!(cursor > 0 && cursor <= n); + max_step = max_step.max(calls.get()[0] + calls.get()[1] - before); + } + let v = calls.get(); + for i in 0..3 { + totals[i] += v[i]; + } + samples.push(v[0] + v[1]); + } + samples.sort_unstable(); + println!( + "{n},{k},{:.4},{:.4},{:.4},{:.4},{},{},{},{}", + totals[0] as f64 / 10_000.0, + totals[1] as f64 / 10_000.0, + 2.0 * (totals[0] + totals[1]) as f64 / (10_000 * k) as f64, + totals[2] as f64 / 10_000.0, + samples[4999], + samples[9899], + samples[9999], + max_step + ); + } + } + } +}