v2.3.0 — Clean-room trig, CDF tail fixes & mpmath accuracy characterization

Added

• Accuracy characterization vs mpmath (#46): new
  tools/accuracy_sweep.cpp (deterministic, bit-exact CSV of scalar and
  FORCE_VECTORIZED batch pdf/logpdf/cdf/quantile over
  support-aware grids, 19 distributions × 3 instances, tails to
  p = 1e-300) and tools/accuracy_vs_mpmath.py (mp.dps = 50 oracle,
  42 unconditional self-checks, law-of-F tail normalization), generating
  docs/ACCURACY_CHARACTERIZATION.md. Replaces the issue's pylibstats
  approach, which pins to a released libstats and would characterize the
  wrong code. The oracle carries its own large-parameter machinery where
  mpmath's betainc/lower gammainc hang or raise (clean-room
  continued-fraction incomplete beta, upper-gamma complement, far-tail
  asymptotic guards, safeguarded log-log false-position quantile
  inversion — each anchored by self-checks). Confirms the four pinned
  accuracy gates and surfaces real findings: inconsistent batch NaN
  propagation across 8 distributions, NaN/clamped returns at +-inf
  inputs, quantile NaN/saturation at extreme p, and large-parameter
  CDF accuracy limits (binomial n = 1e6 off 1.3e-2 at the mean) — see
  the doc's Findings section.

Changed

• LogNormal CDF lower tail no longer collapses (#49): every path
  computed 0.5·(1+erf(z/√2)), whose lower tail dies on the 1+erf
  cancellation floor (~1.1e-16 absolute) regardless of erf quality — the
  reason #49's 1-ULP erf swap changed nothing — and returns exactly 0 once
  std::erf saturates at −1 (z ≲ −8.3). True max relative error over the
  issue's benchmark grid was 1.0, not the reported 2.62e-7 (a
  metric-flooring artifact). z < 0 now routes through 0.5·erfc(−z/√2)
  in detail::normal_cdf (whose analysis-helper consumers inherit the
  fix) and both LogNormal batch paths; the SIMD path keeps its vectorized
  erf pipeline with a per-lane erfc recompute below erf-argument −1. New
  mpmath-oracle gate with references down to F ≈ 1.9e-307, budgeted by the
  achievable-accuracy law rel(F) ~ |ln F|·2⁻⁵² (a flat deep-tail budget is
  unachievable in double for this formulation); measured max 0.49 of
  budget. Gaussian's CDF had the same defect independently — fixed in the
  entry below, same change set.

• Gaussian CDF lower tail no longer collapses (#49 pattern):
  GaussianDistribution never routed through detail::normal_cdf, so all
  five of its own CDF sites (scalar, the three parallel lambdas, and the
  SIMD batch kernel) independently reproduced the #49 cancellation bug.
  Same fix: x < μ routes through 0.5·erfc(−z/√2); the SIMD batch path
  keeps its vectorized erf pipeline with a per-lane erfc recompute below
  erf-argument −1, deliberately using the scalar path's exact argument
  expression so fixed-up lanes are bit-identical to
  getCumulativeProbability(x) (batch and scalar now differ only in the
  −1 ≤ w < 0 plain-erf band, by ≤ 1 ulp of 1). New mpmath-oracle gate
  (test_gaussian_cdf_accuracy, references to F ≈ 1e-290, a bucket
  covering the dedicated standard-normal code path) under the same
  |ln F|·2⁻⁵² accuracy law; fail-first max relative error 1.0, post-fix
  max 0.287 of budget.

• Von Mises CDF rebuilt on the Bessel series, scalar and batch (#51):
  F(x) = (t+π)/(2π) + Σ b_j·sin(j·t) with b_j = I_j(κ)/(j·π·I₀(κ)),
  t = wrap(x−μ), j_max = ⌈10 + 8.5√κ⌉. Coefficients come from Miller's
  backward recurrence normalized by f_j/f₀ — ratios only, no Bessel
  function evaluated anywhere, so no dependency on the platform Bessel
  tier (#47). Replaces a ≥512-point trapezoid (κ ≤ 50, ~1e-6) and a
  wrapped-normal approximation (κ > 50, ~1e-4). Measured against a
  40-digit mpmath quadrature oracle on Zen 4: scalar max 2.2e-16 and batch
  max 8.9e-16 absolute at every κ bucket through 1000; gates checked in
  (tests/test_vonmises_cdf_accuracy.cpp, budgets pinned 2e-15/4e-15,
  demonstrated failing against the old implementation first). The batch
  path evaluates the series with vector_sin (#95) — one batch sine +
  axpy per term — and the CDF dispatch thresholds are un-NEVERed
  (provisional values matching PDF's, pending benchmark). κ > 1000 keeps
  the wrapped-normal fallback; κ = 0 is the exact uniform closed form.

• Von Mises CDF is now μ-centered — behavior fix for μ ≠ 0: the old
  CDF wrapped x into the absolute interval (−π, π] regardless of μ, which
  contradicted the class's documented F(μ) = 0.5 median invariant and was
  inconsistent with getQuantile (whose grid was always μ-centered). The
  series CDF wraps t = x−μ, restoring F(μ) = 0.5 exactly for every μ and
  making CDF/quantile round-trips consistent. Anyone relying on the old
  absolute-domain CDF values for μ ≠ 0 will see different (now correct)
  results.

• x86 vector_cos rebuilt at every tier; vector_sin added at every tier
  (#95): the four x86 kernels were a 7-term Taylor polynomial at ~1e-10
  absolute (~3e5 ULP; wrong sign near zeros of cos — worst case measured
  ~1e19 ULP on the bit lattice at x ≈ −7.5e6). All four are now the
  clean-room quadrant-reduction kernel already proven in vector_cos_neon
  and in libhmm v4.4.0's x86 retarget of the same derivation: measured on
  Zen 4 against 320-bit mpmath references, max 1 ULP, mean 0.022–0.028
  ULP, all tiers, cos and sin; 0 ULP on the specials set. SSE2 and AVX
  use a plain-arithmetic form that loses nothing in the argument reduction
  (every n·π/2-part product is exact by the 30-bit-split construction);
  AVX2/AVX-512 use explicit FMA. Domain contract matches NEON: |x| ≤ 2²³
  vectorized, per-lane std::cos/std::sin fixup beyond, NaN propagates,
  sin(±0) = ±0 with sign preserved. New per-tier ULP gate
  (tests/test_trig_ulp_gates.cpp, references checked in) demonstrated
  failing against the old kernels before the fix landed; it also runs the
  dispatched entry points and a non-lane-multiple tail span. The per-tier
  kernels are now public statics on VectorOps for direct gate access.
  vector_sin unblocks the #51 von Mises CDF batch path.

• Cauchy CDF uses the closed-form arctan (#48): scalar and batch
  getCumulativeProbability no longer delegate to StudentT(1)'s regularized
  incomplete-beta path — the CDF is 1/2 + atan((x−x₀)/γ)/π, one atan per
  element with no iteration. Accuracy improves from the delegation's ~10⁻⁹
  absolute to ~2 ULP, and the lower tail (x < x₀) is computed
  cancellation-free as atan(−γ/(x−x₀))/π, holding full relative accuracy
  down to F ≈ 10⁻¹⁶ where the naive form caps at ~ulp(1/2)/F. The batch
  overload now runs its own auto-dispatch (scalar loop until a vector_atan
  SIMD primitive exists; parallel/work-stealing above the thresholds) instead
  of allocating a transform buffer and delegating. PDF/LogPDF delegation is
  unchanged. New mpmath-referenced accuracy test (dps=50) that the previous
  implementation fails at up to 2.7×10⁻¹⁰.