perf: recalibrate parallel dispatch thresholds across all four SIMD tiers (issue #50)

docs/SIMD_BENCHMARK_RESULTS.md

@@ -14,6 +14,114 @@ Current release validation targets:
 
 Record new release measurements here after each real-machine validation pass.
 
+## v2.0.2 scipy comparison benchmark
+
+Cross-library throughput and accuracy comparison run via `pylibstats/benchmarks/scipy_comparison.py` (pylibstats v0.3.2, scipy 1.18.0, numpy 2.4.4). Results capture the state after the `perf/dispatch-threshold-recalibration` threshold updates.
+
+### Machines
+
+| Machine | CPU | SIMD | OS | Python |
+|---|---|---|---|---|
+| Asus TUF A16 | AMD Ryzen 7 7445HS (Zen 4) | AVX-512 | Windows 11 | 3.12.10 |
+| MacBook Pro 14,1 (2017) | Intel Core i7-7820HQ (Kaby Lake) | AVX2+FMA | macOS Tahoe | 3.14.6 |
+
+### Peak throughput highlights (pylibstats absolute, at optimal batch size)
+
+| Distribution / Op | Zen4 AVX-512 | Kaby Lake AVX2+FMA | Zen4 peak N |
+|---|---|---|---|
+| Exponential log_PDF | 2.2G/s | 490M/s | 30k |
+| Uniform log_PDF | 1.5G/s | 1.1G/s | 20–30k |
+| Gaussian log_PDF | 2.1G/s | 760M/s | 25k |
+| Exponential PDF | 993M/s | 303M/s | 30k |
+| Uniform PDF | 1.5G/s | 966M/s | 20–30k |
+
+All peak measurements were taken at sizes below the L2→L3 boundary (~50k elements on Zen4; see `SIMD_OPTIMIZATION_REFERENCE.md §Cache hierarchy effects`).
+
+### Throughput ratios vs scipy (selected, at N=100k)
+
+Speedup ratios vary with N. N=100k is representative for sustained throughput above the L3 cache fill threshold on both machines.
+
+| Distribution | Op | Zen4 | Kaby Lake |
+|---|---|---|---|
+| Exponential | log_PDF | 29× | 11× |
+| Gaussian | log_PDF | 6× | 22× |
+| Gamma | log_PDF | 18× | 9× |
+| VonMises | log_PDF | 25× | 24× |
+| StudentT | PDF | 35× | 23× |
+| Uniform | PDF | 19× | 25× |
+| Cauchy | CDF | 0.2× | 1.1× | |
+| VonMises | CDF | 0.1× | 0.1× |
+| Binomial | CDF | 0.9× | 0.4× |
+
+Negative ratios (<1×) indicate distributions with structural algorithm limitations; see `SIMD_OPTIMIZATION_REFERENCE.md §Known structural performance ceilings`.
+
+### Accuracy vs scipy (max relative error, N=50k uniform grid)
+
+All values are bit-identical between Zen4 and Kaby Lake **except VonMises PDF/LogPDF**, which differ due to the Bessel function Tier 1/Tier 2 selection:
+
+| Distribution | pdf | log_pdf | cdf |
+|---|---|---|---|
+| Gaussian | 1.0×10⁻¹⁵ | 4.4×10⁻¹⁶ | 9.7×10⁻¹⁵ |
+| Exponential | 2.2×10⁻¹⁶ | 0 | 1.5×10⁻¹⁴ |
+| Laplace | 0 | 1.4×10⁻¹³ | 0 |
+| Uniform | 0 | 0 | 2.2×10⁻¹⁶ |
+| StudentT | 1.5×10⁻¹⁵ | 5.4×10⁻¹⁶ | 2.7×10⁻⁹ |
+| Gamma | 2.0×10⁻¹⁵ | 3.3×10⁻¹⁶ | 1.8×10⁻⁹ |
+| LogNormal | 7.3×10⁻¹⁵ | 2.7×10⁻¹⁵ | 2.6×10⁻⁷ |
+| VonMises (Zen4/MSVC) | 8.4×10⁻¹⁶ | 4.9×10⁻¹¹ | 4.6×10⁻⁶ |
+| VonMises (Kaby Lake/AppleClang) | 2.3×10⁻⁹ | 3.3×10⁻⁹ | 4.6×10⁻⁶ |
+
+The StudentT/Gamma CDF errors (~10⁻⁹) and LogNormal CDF error (~10⁻⁷) are consistent between machines and reflect approximation limits in the regularised incomplete-beta and erfc paths, not SIMD errors. See issues #49 and related.
+
+### Benchmark command
+
+```bash
+# From the pylibstats repository root:
+python benchmarks/scipy_comparison.py --sizes 1000,10000,100000,1000000
+
+# Fine-grained sweep to detect dispatch threshold issues (issue #50 methodology):
+python benchmarks/scipy_comparison.py --sizes 1000,5000,10000,20000,30000,40000,50000,60000,70000,80000,90000,100000
+```
+
+---
+
+## NEON threshold recalibration sweep
+
+Sweep run via `pylibstats/benchmarks/scipy_comparison.py` (pylibstats v0.3.2, scipy 1.17.1, numpy 2.4.4) to identify and resolve `kNeon` dispatch threshold floor artefacts. Branch: `perf/dispatch-threshold-recalibration`.
+
+### Machine
+
+| Machine | CPU | SIMD | OS | Python |
+|---|---|---|---|---|
+| MacBook Pro (14,1, 2017) | Apple M1 | NEON | macOS | 3.14.6 |
+
+### Key finding: L1 cache boundary at N ≈ 5k
+
+M1 L1 data cache is 64KB per performance core. Two-array footprint (input + output doubles) exceeds L1 at N ≈ 4096 elements, producing a systematic throughput trough at N ≈ 5k across all operations regardless of threshold. Operations with `kNeon` threshold=64 (profiler floor) have parallel always active; their N=5k troughs are entirely cache-driven. See `SIMD_OPTIMIZATION_REFERENCE.md §L1 data cache boundary on NEON/M1`.
+
+### Threshold corrections
+
+Three operations showed genuine dispatch troughs with minimum at N=10k (above the L1 boundary):
+
+| Distribution / Op | Old kNeon | New kNeon | Evidence |
+|---|---|---|---|
+| Laplace PDF | 6144 (floor artefact) | 35000 | Trough at N=10k (268M/s vs 330M/s VECTORIZED); parallel entry first exceeds VECTORIZED at N=35k (+20%) |
+| Rayleigh PDF | 10000 | 20000 | Trough at N=10k (179M/s vs 214M/s VECTORIZED); recovery at N=15k; 20k conservative safe zone |
+| LogNormal PDF | 10000 | 25000 | Trough extends N=10k–15k (180M/s vs 207M/s VECTORIZED); recovery at N=20k; 25k conservative safe zone |
+
+### Benchmark commands
+
+```bash
+# From the pylibstats repository root:
+# Standard sweep to identify candidates:
+python benchmarks/scipy_comparison.py --sizes 1000,5000,10000,20000,30000,40000,50000,60000,70000,80000,90000,100000
+
+# Focused 5k-resolution sweep for trough pinpointing and threshold verification:
+python benchmarks/scipy_comparison.py --sizes 1000,5000,10000,15000,20000,25000,30000,35000,40000,45000,50000
+```
+
+---
+
 ## v1.5.x historical summary
 
 v1.5.0 introduced native transcendentals on AVX2, NEON, and AVX-512 and changed `simd_verification` to report per-operation-type geometric means rather than one composite number.

docs/SIMD_OPTIMIZATION_REFERENCE.md

@@ -21,6 +21,150 @@ Resolved work from v1.4.0 through v2.0.0:
 | NEON | Apple Silicon | M1 validation path |
 | AVX-512 | AMD Zen 4 / supported Intel | Windows validation path |
 
+## Detecting threshold miscalibration via external benchmarks
+
+`scripts/PROFILING_METHOD.md` documents the authoritative `strategy_profile`-based calibration procedure. External language bindings (e.g. `pylibstats`) provide a complementary and faster signal: run a throughput sweep across a dense size grid and look for throughput troughs.
+
+### The trough-at-threshold signature
+
+When a dispatch threshold is too low, the parallel strategy fires before its overhead is amortised. The benchmark exposes this as a V-shaped throughput trough whose minimum occurs at the threshold value:
+
+- Throughput rises with N up to just below the threshold (pure SIMD).
+- At the threshold, throughput drops sharply as threading overhead dominates.
+- Throughput recovers as N grows and parallel work amortises the overhead.
+
+Example from the v2.0.2 AVX-512 recalibration (perf/dispatch-threshold-recalibration):
+
+| Distribution / Op | Threshold | Trough N | Trough throughput | Pre-trough |
+|---|---|---|---|---|
+| Laplace log_PDF | 64 (profiler floor) | 10k | 129M/s | 229M/s at N=5k |
+| Laplace log_PDF | 25k (first pass) | 25k | 170M/s | 433M/s at N=20k |
+| Laplace log_PDF | 50k (final) | none | — | monotone above N=50k |
+| Uniform CDF | 128 | 10k | 109M/s | 833M/s at N=5k |
+| Uniform PDF/LogPDF | 50k | 50k | 438M/s | 1.3G/s at N=40k |
+
+The Laplace log_PDF case illustrates the diagnostic loop: a first-pass threshold of 25k moved the trough from N=10k to N=25k rather than eliminating it. A finer sweep (5k resolution) pinpointed that the threshold only amortises at N=45–50k, motivating the final value of 50k.
+
+### Dispatch effect vs cache boundary effect
+
+A cache hierarchy boundary produces a visually similar throughput drop but cannot be eliminated by adjusting thresholds:
+
+- **Dispatch effect**: trough disappears when threshold is set to NEVER.
+- **Cache boundary**: trough persists even with NEVER threshold.
+
+Example: Uniform CDF on Zen4 (AVX-512) shows a throughput drop from ~880M/s at N=45k to ~550M/s at N=50k with NEVER threshold. This is the L2→L3 boundary: two-array footprint (input + output) of 45k doubles = 720KB (fits in Zen4 L2), 50k doubles = 800KB (does not). No threshold change can raise throughput above the L2 ceiling at sizes that exceed L2 capacity.
+
+A 50k threshold previously placed parallel overhead on top of a cache miss (463M/s instead of 552M/s). Setting NEVER removes the overhead and exposes the true cache-limited floor.
+
+### Profiler floor artefacts in the threshold table
+
+A threshold of 64 in `kAvx512` or any architecture table is a strong signal that the profiler measurement floor was hit rather than a real crossover being measured. See `scripts/PROFILING_METHOD.md §Timer jitter and the sub-64 measurement floor` for the mechanism. When a 64 threshold produces a trough in an external benchmark, replace it with a measured amortisation point from either a fine-grained sweep or a fresh `strategy_profile --large` run.
+
+In the v2.0.2 recalibration, `kAvx512` Laplace PDF/LogPDF thresholds of 64 were both documented as profiler floor artefacts and confirmed by the external benchmark trough at N=10k.
+
+### L1 data cache boundary on NEON/M1
+
+The M1 performance core has a 64KB L1 data cache. The two-array footprint (input + output doubles) crosses the L1→L2 boundary at N ≈ 64KB ÷ 16 bytes = 4096 elements. Throughput sweeps on M1 show a systematic drop at N ≈ 5k regardless of dispatch threshold — this mimics a dispatch trough but is cache-driven and cannot be fixed by raising thresholds.
+
+Two indicators distinguish L1 cache effects from dispatch troughs:
+
+- **Trough position**: L1 boundary troughs are pinned at N ≈ 5k on M1 regardless of threshold. A genuine dispatch trough minimum shifts to match the threshold value.
+- **NEVER test**: set the threshold to NEVER and re-sweep. A dispatch trough disappears; an L1 cache trough persists.
+
+Operations with `kNeon` threshold=64 (profiler floor artefact) have parallel always active; their N=5k throughput drop is entirely the L1→L2 transition and does not warrant recalibration. See `SIMD_BENCHMARK_RESULTS.md §NEON threshold recalibration sweep` for a worked example.
+
+### NEVER thresholds for trivial SIMD operations
+
+Distributions whose per-element SIMD cost is very low (Uniform, Laplace at small N) may show that parallel never sustains an advantage within practical batch sizes. For these, NEVER is correct regardless of what the profiler reports at its measurement ceiling. Indicators:
+
+- Parallel never recovers to pre-threshold SIMD throughput across a 1k–100k sweep.
+- The distribution performs arithmetic or a single transcendental per element with no iterative path.
+- The trough-at-threshold depth exceeds 50% throughput loss.
+
+For AVX-512 v2.0.2, Uniform PDF, Uniform LogPDF, and Uniform CDF are all NEVER.
+
+### Clean-entry criterion for iterative calibration
+
+When iterating threshold values, the target state is a parallel entry that exceeds the last VECTORIZED point in the sweep. If the first benchmark measurement at the new threshold T falls below the VECTORIZED level just below T, the parallel path has not yet amortised its threading overhead — raise the threshold further.
+
+Example from the NEON Laplace PDF recalibration (perf/dispatch-threshold-recalibration, M1):
+
+| Threshold | Last VECTORIZED | Parallel entry (N=T) | Entry delta |
+|---|---|---|---|
+| 25k | 327M/s at N=20k | 316M/s | −3% (parallel losing — raise threshold) |
+| 30k | 385M/s at N=25k | 363M/s | −6% (parallel losing — raise threshold) |
+| 35k | 383M/s at N=30k | 460M/s | +20% (parallel winning — correct threshold) |
+
+At 25k and 30k the parallel entry was below VECTORIZED, meaning the threshold fired before overhead amortised. At 35k parallel immediately exceeded VECTORIZED; 35k is the encoded value.
+
+---
+
+## Cache hierarchy effects on batch throughput
+
+AVX-512's higher per-element throughput means the working set grows faster than on AVX2 machines. The L2→L3 and L3→DRAM transitions are therefore more visible in fine-grained throughput sweeps.
+
+Zen4 (Ryzen 7 7445HS) cache topology:
+- L2: ~1MB per core (effective user-data capacity ~800KB)
+- L3: 16MB shared
+- DRAM: DDR5
+
+Observed thresholds in the v2.0.2 benchmark sweep (two-array footprint = input + output):
+
+| Transition | Two-array size | Approx N (doubles) | Typical throughput drop |
+|---|---|---|---|
+| L2 → L3 | ~800KB | ~50k | 30–60% for compute-bound distributions |
+| L3 → DRAM | ~16MB | ~1M | 30–70% for compute-bound distributions |
+
+The L2 boundary is architecture-specific and must not be assumed to hold on AVX2 (Kaby Lake) or NEON (M1). Throughput comparisons across architectures should be made at sizes that fit in each machine's L2 to avoid confounding compute throughput with memory bandwidth.
+
+---
+
+## Cross-architecture accuracy differences
+
+### Bessel function tier selection
+
+`include/core/bessel.h` provides two tiers for `bessel_i0`, `bessel_i1`, and `log_bessel_i0`:
+
+- **Tier 1** (MSVC/GCC/Clang with `LIBSTATS_HAS_CXX17_BESSEL`): delegates to `std::cyl_bessel_i` (C++17 §29.9.3). Achieves <1 ULP against scipy for κ=2.
+- **Tier 2** (AppleClang/macOS, `LIBSTATS_HAS_CXX17_BESSEL` not defined): A&S §9.8.1–9.8.4 polynomial approximation. Documented precision: <1.6×10⁻⁷.
+
+The tier selection is automatic: CMakeLists.txt probes for `std::cyl_bessel_i` and defines `LIBSTATS_HAS_CXX17_BESSEL` only when available. AppleClang does not implement C++17 special functions as of Xcode 16.
+
+### Implication for VonMises accuracy
+
+VonMises PDF and LogPDF accuracy differs between MSVC and AppleClang builds because the normalisation constant `log(2π·I₀(κ))` is computed via different Bessel tiers. This shows up as:
+
+| Platform | VonMises pdf max rel err | VonMises log_pdf max rel err |
+|---|---|---|
+| Zen4 / MSVC (Tier 1) | ~1 ULP (~8×10⁻¹⁶) | ~5×10⁻¹¹ |
+| Kaby Lake / AppleClang (Tier 2) | ~2×10⁻⁹ | ~3×10⁻⁹ |
+
+The `vector_cos` implementations between AVX2 and AVX-512 are algorithmically identical (same 7-term FMA Horner polynomial), confirmed by cross-machine analysis. The accuracy difference is entirely in the scalar Bessel normalisation path, not in any SIMD kernel. See issue #47 for the proposed Tier 2 upgrade.
+
+### Scipy version independence
+
+All other distributions (Gaussian, Exponential, Laplace, Gamma, etc.) produce bit-identical or near-identical accuracy results across Zen4 and Kaby Lake with the same scipy version. The VonMises accuracy difference was initially suspected to be a scipy version artefact (1.17.1 vs 1.18.0); upgrading to 1.18.0 on Zen4 left the Zen4 accuracy unchanged, ruling out scipy.
+
+---
+
+## Known structural performance ceilings
+
+These distributions have throughput limitations that are inherent to their algorithms, not to dispatch thresholds or SIMD kernel quality:
+
+### VonMises CDF
+
+The CDF has no closed form. The current implementation is a scalar integration loop. Throughput: ~200–900k elements/second vs ~30–50M/s for scipy (which uses Cephes adaptive quadrature, also scalar, but more efficient). No SIMD kernel or threshold change will address this. See issue #51 for a precomputed table design.
+
+### Cauchy CDF
+
+Cauchy is implemented as a delegation wrapper over StudentT(ν=1). The CDF therefore evaluates the regularised incomplete-beta function rather than the trivial closed form `arctan((x-x₀)/γ)/π + 0.5`. This makes Cauchy CDF 2–5× slower than scipy on Zen4 where scipy vectorises the arctan path. See issue #48.
+
+### Binomial CDF and PDF
+
+Binomial PMF uses lgamma log-space evaluation in a scalar loop. There is no SIMD log-gamma primitive (`vector_lgamma` is deferred). Binomial CDF uses PMF summation rather than the regularised incomplete-beta approach used by scipy. Both paths cap throughput at ~5–16M elements/second regardless of batch size or architecture. See issue #52.
+
+---
+
 ## Deferred work
 
 ### vector_floor and vector_blend