bessel.h Tier 2 fallback (A&S polynomial) limits VonMises accuracy to ~10⁻⁷ on macOS/AppleClang

[OldCrow]

Summary

On macOS/AppleClang builds, LIBSTATS_HAS_CXX17_BESSEL is not defined (AppleClang does not implement std::cyl_bessel_i), so bessel_i0, bessel_i1, and log_bessel_i0 fall back to the Abramowitz & Stegun polynomial approximations in include/core/bessel.h. The documented precision budget for this path is < 1.6×10⁻⁷.

This accuracy floor propagates directly into the VonMises normalisation constant (log(2π·I₀(κ))), producing a constant multiplicative offset in every PDF and log-PDF evaluation on macOS.

Evidence

Cross-machine benchmarks/scipy_comparison.py run (pylibstats v0.3.2, same scipy 1.18.0, κ=2.0):

Machine	SIMD	Bessel path	VonMises pdf error	VonMises log_pdf error
Asus TUF A16 (Windows, MSVC)	AVX-512	Tier 1: std::cyl_bessel_i	8.4×10⁻¹⁶ (≈ 1 ULP)	4.9×10⁻¹¹
MacBook Pro 2017 (macOS, AppleClang)	AVX2+FMA	Tier 2: A&S polynomial	2.3×10⁻⁹ (!)	3.3×10⁻⁹ (!)

All other distributions show bit-identical accuracy across the two machines, confirming vector_cos and vector_exp are not the source. The ~10⁻⁹ error on macOS is the gap between the A&S polynomial and scipy's Bessel reference for κ=2.0, well within the Tier 2 budget but visible against a 1-ULP reference.

The scipy version was also ruled out: upgrading from 1.17.1 to 1.18.0 on the Windows machine left the Zen4 result unchanged at 8.4×10⁻¹⁶.

Root cause

bessel.h lines 68–87 (Tier 2 bessel_i0): the A&S 9.8.1/9.8.2 polynomial has single-precision-level coefficients (7 significant digits), which limits the result to ~10⁻⁷ absolute error. For κ ≤ 3.75 (the common range including κ=2.0 used in benchmarks), log_bessel_i0 calls std::log(bessel_i0(x)), directly inheriting this error.

Potential fixes

In order of increasing implementation effort:

1. Replace A&S with a Chebyshev or minimax polynomial fitted to double precision. Cephes / Boost.Math both have double-precision I₀ implementations in the public domain. This would bring Tier 2 to < 1 ULP without changing the interface.

2. Implement log_bessel_i0 directly (without calling bessel_i0 first) using the log-space A&S expansion, similar to what Tier 1 already does for the x > 700 overflow guard. This avoids the intermediate exp→log round-trip but does not improve the polynomial coefficients.

3. Accept the current floor and document it explicitly in the benchmark notes (the VonMises row note already flags the CDF; the PDF/log_pdf rows do not currently mention the Bessel tier difference).

Affected platforms

• macOS (all versions) — AppleClang does not implement std::cyl_bessel_i
• Any future platform where LIBSTATS_HAS_CXX17_BESSEL is not detected

Not affected

• Windows (MSVC 2017 15.5+) — Tier 1 active
• Linux GCC 6.1+ / Clang 17+ — Tier 1 active (subject to platform libc++ providing the special functions)