Deléglise–Rivat $O(X^{2/3})$ Mertens Engine

A self-contained, header-only C++20 implementation of the Deléglise–Rivat combinatorial method with inclusion–exclusion reductions and micro-architectural optimizations for evaluating the Mertens function: $$ M(X) = \sum_{n=1}^X \mu(n) $$

Based on and extending the mathematical framework described in Greg Hurst's preprint ("Practical Computations of the Mertens Function: $M(10^{24})$ and $M(10^{25})$", arXiv:2607.07566).

1. Inclusion–Exclusion Reductions:

   • Outer Parity Split & Mod-6 Wheel Reduction: Evaluates $S(y, u) - S(y/2, u) - S(y/3, u) + S(y/6, u)$ across square-free $k$ coprime to 6, cutting outer terms by $>69%$.
   • $S_1$ Parity Cancellation: Cancels even quotient lookups across surviving odd/even intervals (stride-2 queries).
   • $S_2$ Piecewise Reduction: Partitions $j \le A, (j, 6)=1$ into 8 exact branchless arithmetic intervals.

2. Compressed Hierarchical Mertens Table (CHMT) Algorithm (Unified Engine for all $X > 10^5$):

   • Seamless Multi-Tier Scaling: Dynamically sizes u_dense = std::min(u_total, 600M). For $X \le 10^{14}$ ($u \le 600\text{M}$), it operates as a 100% dense table with 0 compressed chunks (zero overhead). For massive $X \ge 10^{15}$, it seamlessly engages Tier 2 (2-bit bit-plane chunks), scaling $u$ up to 8 Billion entries in ~3.1 GB RAM.
   • 13.3× Outer Term Reduction: Slashes the outer summation workload for $10^{16}$ from 16.7 million terms down to 1.25 million terms ($N = X / u$).
   • Constant-Time Hardware Popcount Reconstruction: Evaluates $M(q)$ in single-digit CPU clock cycles via 4 branchless 64-bit hardware popcount instructions (__builtin_popcountll(nz & ~sg) and __builtin_popcountll(nz & sg)), requiring zero memory decompression.
   • Branchless 8-Element Bit-Plane Packing: Exploits $\mu(8k+4)=\mu(8k+8)=0$ and $\mu(8k+2)=-\mu_{\text{odd}}[2k+1]$ to pack 8 numbers per iteration directly from odd sieve bytes, achieving a 13× speedup in table compression.
   • Compile-Time Reciprocal Prime Sieving: Replaces runtime 8-cycle UDIV instructions for small primes ($p \in {3, 5, \dots, 71}$) with 1-cycle UMULH compile-time fixed-point reciprocals (sieve_prime_const<P>).
   • Two-Tier $S_1$ Execution Engine: Resolves initial large quotients via compressed popcount reconstruction and the remaining 99.9% of quotients ($q \le u_{\text{dense}}$) through an 8-way unrolled direct read from M_dense.

3. Micro-Architectural & SIMD Acceleration:

   • Zero-Allocation Compact Sieve (FastSieve): Generates $M(n)$ prefix sums directly from odd-sieved segments via sequential 64-bit quad writes (int16_t quads), capping the $\mu$ table allocation to $<40\text{ MB}$ ($A_{\max}$) instead of $u$ (600 MB) and saving $>1.2\text{ GB}$ of DRAM traffic.
   • 4-Way Pipelined NEON $S_2$ Runner: Unrolled 4-way (8 coprime-to-6 values per iteration), maintaining $j$ strictly in SIMD registers with vector step increments (vaddq_f64), pipelining 4 vdivq_f64 instructions, and using vcvtq_s64_f64 hardware conversions.
   • Pure SIMD $S_1$ and Even-$n$ Correction: Eliminates division overhead from software prefetching and evaluates $S_1$ arithmetic with SIMD vector additions and hardware integer truncation.
   • Unified Dynamic Parallel Dispatch: Merges all 4 OpenMP ranges into a single unified term vector with dynamic chunk scheduling (schedule(dynamic, 8)), eliminating barrier syncs and fixing thread starvation on small $k$.
   • Flat 16-Bit Cache (int16_t): Uses single-cycle ldrh lookups without multi-level cache table indirections.
   • Dynamic Hardware-Aware $c_x(X)$ Tuning: Automatically balances $S_1$ memory traffic and $S_2$ streaming throughput ($c_x = 0.95$ for $10^{15..16}$, $0.95$ for $10^{11..14}$, $0.70$ for $X \le 10^{10}$).

• Clang / GCC supporting C++20
• OpenMP (libomp via Homebrew on macOS: brew install libomp)

make

make test
# Or:
./test_dr

# Evaluate a single value:
./bench_dr 1000000000000000 8

# Benchmark sweep from 10^1 to 10^15:
./bench_dr --benchmark 15 8

# Comparative benchmark vs O(X^(3/4)) Dirichlet DP Engine:
./bench_dr --compare 14 8

The engine provides CUDA backend support with GPU-accelerated $S_1$ and $S_2$ evaluations:

• __ldg Read-Only Caching: Directly streams lookups of $M(n)$ from VRAM (~900MB).
• __constant__ Memory LUT Evaluators: Fast broadcast caching for $S_2$ piecewise summands with zero register pressure.
• Branchless Piecewise $S_2$ Modulo-6: Evaluates arithmetic ranges without branch divergence.
• Warp-Level Parallel Reduction: Uses __shfl_down_sync and shared memory tree reductions for high-occupancy 64-bit atomic accumulation.

# Build CUDA binaries targeting NVIDIA T4 (sm_75):
make cuda

# Or specify a custom architecture (e.g. sm_80 for A100, sm_89 for RTX 4090, sm_90 for H100):
make cuda CUDA_ARCH=sm_80

# Run CUDA unit test suite:
make test-cuda

# Evaluate a single value on GPU with detailed timing breakdown:
./bench_cuda 1000000000000000

# Run CUDA benchmark sweep:
./bench_cuda --benchmark 16

# Compare CPU vs CUDA GPU performance:
./bench_cuda --compare 16

Measured 2026-09-20 on Apple M1 Mac mini (8T, 16GB) with the current engine (scalar double-division quotient loops by default with opt-in NEON via MERTENS_USE_NEON_DIV, small-prime {3,5,7,11} sieve stencil, tunable OpenMP chunk MERTENS_DYN_CHUNK=64). MERTENS_FX/MERTENS_CX env overrides select the sieve balance fx and split cx; defaults are fx=0.32 (cap 1.5G) and cx=1.35 ($\ge 10^{15}$) / 1.15 ($\ge 10^{13}$) / 0.95 ($\ge 10^{11}$) / 0.70 (below). Best 2026-09-20 tuning on this box for $10^{15}$ is MERTENS_FX=0.2 MERTENS_CX=1.15. Expect $\pm 10$–$15%$ run-to-run variance (frequency ramp / heat soak on Apple Silicon); figures below are single-sweep runs with ./bench_dr --benchmark 15 8.

Small inputs (same 2026-09-20 sweep, default vs tuned — tuning is a no-op below $u = X$, shown for completeness):

Retained from the earlier T4 run — not re-measured here (no CUDA-capable device on the Apple M1 host used above).