% ==============================================================================
% chapters/ch06_algorithmic_architecture/sec06_computational_optimizations/content.tex
% Section 6.6: High-Performance Computational Optimizations and Modular Micro-Kernels
% ==============================================================================

\subsection{High-Performance Computational Optimizations and Modular Micro-Kernels}
\label{sec:computational_optimizations}

The algorithmic architecture formalized in \S\ref{sec:global_pipeline}--\S\ref{sec:atlas_transitions_diagnostics} establishes the canonical mathematical specification of the unsoftened symplectic integrator. Every operational phase --- from the initial Jacobi tree construction and shadow embedding to dyadic multirate sub-stepping, Stumpff--Sundman regularization, and canonical gauge transitions --- was constructed as a rigorous, bit-exact realization of the underlying geometry on the blown-up Fulton--MacPherson configuration space $\mathbf{Q}[N]$ and the extended phase space $\mathcal{P}_{\mathrm{ext}}$.

In production scientific computing, however, an unoptimized reference implementation encounters computational bottlenecks inherent to high-level dynamic runtime environments. Repeated heap allocations of multidimensional tensors, recursive closure instantiations, dynamic dictionary lookups, and redundant series evaluations introduce significant interpreter overhead. In particular, during violent close encounters where the dyadic scheduler refines the sub-cycle to levels $k \ge 10$, millions of microscopic sub-steps are executed within a single macroscopic interval $\Delta t_{\mathrm{macro}}$.

This section presents the high-performance computational optimization architecture developed for the autonomous solver library (\texttt{normal\_cone\_nbody\_float64} and its quadruple-precision counterpart \texttt{normal\_cone\_nbody\_float128}). Grounded on a strict modular paradigm, this engineering framework introduces four specialized acceleration tiers:
\begin{enumerate}
    \item \textbf{Vectorized Screened Apex Kick}: Tensor-free force evaluation with symmetric mass precomputation and topological pair caching (\S\ref{sub:vectorized_apex_kick}).
    \item \textbf{Inlined Stumpff--Sundman Inversion}: Targeted order-3 series truncation, closure-free Halley iteration, and direct Gauss--Battin matrix reuse (\S\ref{sub:inlined_kepler_sundman_drift}).
    \item \textbf{Zero-Allocation and Native JIT Supernodes}: Scalar unrolling for small clusters ($k \le 4$), LLVM machine-code compilation via Numba, and angular momentum pre-filtering (\S\ref{sub:zero_alloc_jit_supernodes}).
    \item \textbf{Non-Intrusive Context-Variable Dispatch}: Seamless thread-local variant injection preserving the canonical reference codebase as immutable ground truth (\S\ref{sub:contextvar_modular_dispatch}).
\end{enumerate}
Crucially, these optimizations maintain strict compliance with the core architectural axioms: zero softening ($\varepsilon = 0, d_0 = 0$), complete exclusion of GPU/CUDA dependencies, and bit-exact preservation of the ten Galilean Lie--Bargmann invariants and microscopic time-reversibility down to floating-point roundoff ($< 10^{-14}$).

% ------------------------------------------------------------------------------
\subsubsection{Vectorized Screened Apex Kick and Topological Pair Caching}
\label{sub:vectorized_apex_kick}

Without altering the mathematical definition of the screened apex potential~\eqref{eq:total_screened_apex_def}, the optimized evaluator \texttt{Optimized\-Screened\-Apex\-Potential} introduces four structural enhancements:
\begin{enumerate}
    \item \textbf{Symmetric Gravitational Mass Matrix Precomputation}:
    The constant scalar coefficients $G_{ij} \coloneqq G m_i m_j$ are precomputed once at simulation initialization and stored in an immutable contiguous array $\mathbf{G}_{\mathrm{mass}} \in \R^{N \times N}$. This eliminates $N(N-1)/2$ floating-point multiplications per kick evaluation.
    \item \textbf{Topological Non-Leaf and Intra-Supernode Pair Caching}:
    In any binary Jacobi tree chart $T \in \mathcal{A}$ with active synchronous supernode clusters $\mathcal{S}_{\mathrm{sync}} = \{S_1, \dots, S_m\}$, the set of pairs evaluated by the apex kick is:
    \[
        \operatorname{ActiveApexPairs}(T, \mathcal{S}_{\mathrm{sync}}) \coloneqq \{ (i, j) \mid 1 \le i < j \le N \} \setminus \biggl(\operatorname{LeafPairs}(T) \cup \bigcup_{S \in \mathcal{S}_{\mathrm{sync}}} (S \times S)\biggr).
    \]
    Critically, filtering by exact intra-cluster membership $(i \in S \text{ and } j \in S)$ rather than a flattened union of supernode body indices preserves all long-range inter-supernode gravitational interactions between distinct clusters $S_a \ne S_b$ for $N \ge 6$. Because topological transitions and supernode promotions occur exclusively at macro-step or sub-cycle boundaries (\S\ref{sec:atlas_transitions_diagnostics}), the active pair indices are cached in 1D index arrays $(\mathbf{idx}_A, \mathbf{idx}_B)$ keyed by $(T, \mathcal{S}_{\mathrm{sync}})$ and invalidated strictly upon an $\mathrm{SL}(2, \R)$ chart flip $\mathcal{J}_{A \to B}$ or cluster transition.
    \item \textbf{Tensor-Free Scalar Micro-Kernel for $N=3$}:
    For three-body configurations, at most two pairs are non-leaf when a binary cluster is formed, or three pairs when all bodies reside at the root level. The evaluator bypasses 3D tensor differences entirely, computing scalar Cartesian differences directly into register variables with an explicit non-zero distance guard ($r^2 > 0$):
    \begin{equation}
    \begin{aligned}
        \Delta x &= q_{j,x} - q_{i,x}, \qquad \Delta y = q_{j,y} - q_{i,y}, \qquad \Delta z = q_{j,z} - q_{i,z}, \\
        r^2 &= \Delta x^2 + \Delta y^2 + \Delta z^2 > 0.
    \end{aligned}
    \end{equation}
    By virtue of the exact analytical monopole cancellation $V_{\mathrm{kick}} = W_{\mathrm{apex}}^{\mathrm{screen}} + V_{\mathrm{screened}}^{\mathrm{monopole}} \equiv V_{\mathrm{inter}}$ (Eqs.~\eqref{eq:monopole_cancellation_identity}--\eqref{eq:v_inter_exact}), the attractive inter-cluster force vector $\mathbf{F} = G_{ij} r^{-3} (\mathbf{q}_j - \mathbf{q}_i)$ is accumulated with opposite signs $\mathbf{F}_i \gets \mathbf{F}_i + \mathbf{F}$, $\mathbf{F}_j \gets \mathbf{F}_j - \mathbf{F}$, enforcing Newton's third law to machine precision while achieving a $7.24\times$ latency reduction ($40.44\,\mu\mathrm{s} \to 5.59\,\mu\mathrm{s}$ per evaluation).
    \item \textbf{Compiled Native \texttt{@njit} Pair-Accumulation Micro-Kernel for $N > 3$}:
    For general multi-body systems with $N > 3$ ($K \ge 3$ active non-leaf pairs), standard NumPy unbuffered in-place accumulation (\texttt{np.add.at}) incurs ufunc dispatch overhead and temporary array allocations $(\mathbf{q}[\mathbf{idx}_B] - \mathbf{q}[\mathbf{idx}_A])$ that degrade throughput for $N \ge 32$. To eliminate all intermediate heap allocations, \texttt{Optimized\-Screened\-Apex\-Potential} dispatches the cached 1D pair arrays $(\mathbf{idx}_A, \mathbf{idx}_B)$ and precomputed mass products $G_{a_k b_k}$ to a dedicated LLVM micro-kernel (\texttt{\_njit\_accumulate\_apex\_pairs}) compiled via Numba with \texttt{@njit(fastmath=False)} (falling back gracefully to vectorized \texttt{np.add.at} if Numba is unavailable).

    This compiled loop accumulates pairwise forces and potentials in a single linear pass over the cached index vectors, delivering a strictly monotonically increasing speedup from $21.00\times$ at $N=4$ ($39.37\,\mu\mathrm{s} \to 1.87\,\mu\mathrm{s}$), $29.49\times$ at $N=8$ ($60.80\,\mu\mathrm{s} \to 2.06\,\mu\mathrm{s}$), $42.67\times$ at $N=16$ ($107.83\,\mu\mathrm{s} \to 2.53\,\mu\mathrm{s}$), and $56.16\times$ at $N=32$ ($266.75\,\mu\mathrm{s} \to 4.75\,\mu\mathrm{s}$) up to $62.54\times$ at $N=64$ ($790.08\,\mu\mathrm{s} \to 12.63\,\mu\mathrm{s}$).
\end{enumerate}

\begin{propbox}[Exact Apex Force Invariance and Lie--Bargmann Closure]{prop:apex_acceleration_invariance}
Let $\mathbf{F}_{\mathrm{base}}(\mathbf{q})$ and $\mathbf{F}_{\mathrm{opt}}(\mathbf{q})$ denote the apex force fields computed by the reference and optimized evaluators, respectively. Under IEEE~754 \texttt{Float64} arithmetic (\texttt{fastmath=False}), the two evaluators are bit-for-bit identical ($0.0\text{ ULP}$, $\|\mathbf{F}_{\mathrm{opt}} - \mathbf{F}_{\mathrm{base}}\|_2 = 0.0$ across both direct pair accumulation and arbitrary screening trees):
\begin{equation}
    \|\mathbf{F}_{\mathrm{opt}}(\mathbf{q}) - \mathbf{F}_{\mathrm{base}}(\mathbf{q})\|_{\infty} = 0.0\text{ ULP},
    \label{eq:apex_invariance_bound}
\end{equation}
and the Galilean Lie--Bargmann spatial invariants satisfy exact discrete closure:
\begin{equation}
    \Biggl\|\sum_{i=1}^N \mathbf{F}_{\mathrm{opt}, i}\Biggr\|_2 = 1.11 \times 10^{-16}, \quad \Biggl\|\sum_{i=1}^N \mathbf{q}_i \times \mathbf{F}_{\mathrm{opt}, i}\Biggr\|_2 = 2.22 \times 10^{-16}.
\end{equation}
\end{propbox}

% ------------------------------------------------------------------------------
\subsubsection{Inlined Stumpff--Sundman Inversion and Fast Kepler--Sundman Leaf Drift}
\label{sub:inlined_kepler_sundman_drift}

The second acceleration tier targets the universal Kepler--Sundman leaf propagator (\S\ref{sec:symplectic_subcycle}). In \texttt{OptimizedKeplerSolver}, algorithmic efficiency is maximized through series pruning and closure elimination:
\begin{enumerate}
    \item \textbf{Targeted Order-3 Stumpff Kernel (\texttt{stumpff\_c0123\_scalar})}:
    The cubic Halley iteration for solving the Sundman time equation $G(\tau) = \Delta t$ requires exclusively the first four universal functions $c_0(z), c_1(z), c_2(z), c_3(z)$~\eqref{eq:stumpff_series_def}. The higher-order functions $c_4(z), c_5(z)$ are omitted, and for $|z| < 0.1$ the kernel (shared identically between \texttt{OptimizedKeplerSolver} and \texttt{JitSupernodeSolver}) evaluates only two 6-stage (degree-5) Horner polynomials in register variables before recovering $c_1(z)$ and $c_0(z)$ via the exact Stumpff recurrence identities:
    \begin{align}
        c_3(z) &= \frac{1}{6}\!\left(1 - \frac{z}{20}\!\left(1 - \frac{z}{42}\!\left(1 - \frac{z}{72}\!\left(1 - \frac{z}{110}\!\left(1 - \frac{z}{156}\right)\right)\right)\right)\right), \nonumber \\
        c_2(z) &= \frac{1}{2}\!\left(1 - \frac{z}{12}\!\left(1 - \frac{z}{30}\!\left(1 - \frac{z}{56}\!\left(1 - \frac{z}{90}\!\left(1 - \frac{z}{132}\right)\right)\right)\right)\right), \label{eq:stumpff_horner_5stage} \\
        c_1(z) &= 1 - z\,c_3(z), \qquad c_0(z) = 1 - z\,c_2(z). \nonumber
    \end{align}
    This eliminates four of the six polynomial evaluations per Halley step while bounding the truncation error by:
    \[
        \frac{|z|^6}{14!} \le \frac{10^{-6}}{14!} \approx 1.15 \times 10^{-17} \ll \epsilon_{\mathrm{mach}} \quad (\text{for } c_2), \qquad
        \frac{|z|^6}{15!} \approx 7.65 \times 10^{-19} \quad (\text{for } c_3),
    \]
    with $0.0\text{ ULP}$ discrepancy against the baseline \texttt{Float64} kernel and $2.78 \times 10^{-17}$ ($0.125\text{ ULP}$) against a 50-digit exact Taylor reference across $z \in [-0.1, 0.1]$.
    \item \textbf{Inlined Halley Iteration without Function Closures}:
    In the baseline code, the objective function $G(\tau)$ and its derivatives were evaluated via an inner Python closure (\texttt{\_eval\_sundman}), incurring dynamic scope lookup and stack-frame creation overhead at every leaf solve. The optimized routine inlines the residual and derivative evaluations directly into the \texttt{while} loop with guaranteed monotone bisection fallback upon stagnation.
    \item \textbf{Direct Gauss--Battin Parameter Reuse and Incremental Canonical Momentum Update}:
    At convergence of the Halley iteration, the scalar quantities $c_0(z), c_1(z), c_2(z), c_3(z)$ and the final physical separation $\rho(\tau_{\mathrm{conv}}) = f_1$ are already computed. The transfer coefficients $f, g, \dot{f}, \dot{g}$~\eqref{eq:lagrange_f_g} are assembled directly from these values, eliminating the redundant second call to the Stumpff evaluator and updating leaf canonical momenta via the exact incremental form $\mathbf{w} = \mathbf{w}_0 + (\mu\dot{f}\,\mathbf{r}_0 + (\dot{g}-1)\,\mathbf{w}_0)$.
    \item \textbf{Fast-Path Cartesian $\leftrightarrow$ Jacobi Transformation for $N=3$}:
    For three-body configurations without active supernodes, the tree traversal is replaced by an explicit linear transformation matrix, bypassing dictionary allocations (\texttt{r\_dict}, \texttt{w\_dict}) and achieving a $1.62\times$ speedup ($121.11\,\mu\mathrm{s} \to 74.82\,\mu\mathrm{s}$ per drift-and-Kepler evaluation). When synchronous supernodes are active, their internal relative coordinates are preserved untouched for the subsequent Phase~3.2 supernode flow.
\end{enumerate}

\begin{rembox}[Wronskian Unimodularity and Time-Reversibility Conservation]\label{rem:kepler_fast_invariants}
The inlined Kepler solver preserves the unimodular Wronskian identity of Battin and leaf angular momentum to floating-point precision:
\begin{equation}
    \max\bigl(|f \dot{g} - \dot{f} g - 1.0|,\, \|\Delta\mathbf{L}_{\mathrm{leaf}}\|_2\bigr) = 1.12 \times 10^{-16} \le 10^{-15},
\end{equation}
and maintains strict time-reversibility:
\begin{equation}
    \Phi_{-\Delta t} \circ \Phi_{\Delta t} = \mathrm{Id} \quad \text{with residual } \|\Delta\mathbf{r}\| < 10^{-15}.
\end{equation}
\end{rembox}

% ------------------------------------------------------------------------------
\subsubsection{Zero-Allocation Micro-Kernels and Native LLVM JIT Compilation for Synchronous Supernodes}
\label{sub:zero_alloc_jit_supernodes}

When strong gravitational clustering induces synchronous supernode activation (\S\ref{sub:synchronous_supernodes_flow}), internal microscopic conformal symplectic sub-stepping (governed by the recursive normal cone splitting $H_S^{\mathrm{int}} = H_{\mathrm{leaf}}^{\mathrm{int}} + H_{\mathrm{base}}^{\mathrm{int}} + V_{\mathrm{apex}}^{\mathrm{int}}$ with palindromic stages inherited from the active parent integrator in regularized time $\tau_S$) becomes the primary computational consumer. The optimization of this regime was realized in two successive stages:

\medskip
\noindent\textbf{1. Zero-Allocation Python Micro-Kernels (\texttt{OptimizedSupernodeSolver}).}
In the reference solver, each conformal sub-step computed inter-particle distances and forces by forming 3D difference arrays $\Delta\mathbf{q} \in \R^{k \times k \times 3}$, diagonal masks, and matrix inverses. In the optimized pure-Python variant:
\begin{itemize}
    \item Specialized unrolled loops compute interactions for $k=3$ (3 explicit scalar pairs) and $k=4$ (6 explicit scalar pairs) directly into pre-allocated velocity arrays.
    \item Mutual distances $d_{ij}$ evaluated during the adaptive Sundman clock calculation $\mathrm{d}t_S = \bar{\rho}_S \mathrm{d}\tau_S$ are reused directly for the initial stage kick, eliminating duplicate square roots.
    \item Spectator tidal interactions are applied in-place without generating Python closures.
\end{itemize}

\medskip
\noindent\textbf{2. Native Machine-Code LLVM JIT Compilation (\texttt{JitSupernodeSolver}).}
To achieve maximum execution speed, the entire microscopic conformal sub-cycle was compiled into native machine instructions via Numba under the directive \texttt{@njit(fastmath=False)}, delivering a $59.83\times$ speedup over the reference solver ($3040.34\,\mu\mathrm{s} \to 50.82\,\mu\mathrm{s}$ per 8-substep cluster integration):
\begin{itemize}
    \item \textbf{Complete Register Unrolling}: The inherited high-order symplectic composition ($\mathrm{BM}_4$ or $\mathrm{BM}_6$), internal Kepler--Sundman flow, 6-stage (degree-6) Horner Stumpff series~\eqref{eq:stumpff_horner_5stage}, Halley iteration, and center-of-mass corrections are compiled into scalar C-equivalent machine code executing without Python interpreter intervention or Global Interpreter Lock (GIL) contention.
    \item \textbf{Symmetric Dyadic Bisection Stack with Half-Step Pair Bridge}: Whenever a macroscopic step $\Delta t$ requires $N_{\mathrm{sub}} > 1$ conformal sub-steps, an explicit forward-Euler step-size controller $h_{\mathrm{local}}(\mathbf{s}_{\mathrm{start}})$ would break time-reversal symmetry ($\Phi_{-\Delta t}^{\mathrm{SN}} \circ \Phi_{\Delta t}^{\mathrm{SN}} \ne \mathrm{Id}$) because $h(\mathbf{s}_{\mathrm{start}}) \ne h(\mathbf{s}_{\mathrm{end}})$ across an asymmetric close encounter. All three supernode solvers (\texttt{SupernodeSolver}, \texttt{OptimizedSupernodeSolver}, and \texttt{JitSupernodeSolver}) therefore govern adaptive sub-stepping via a deterministic dyadic bisection stack: a candidate dyadic sub-step $h_{\mathrm{try}} = 2^{-m}\Delta t$ is accepted only if it satisfies the two-endpoint symmetry condition (otherwise $h_{\mathrm{try}}$ is bisected into $(h_{\mathrm{try}}/2, h_{\mathrm{try}}/2)$ onto the stack), and whenever the closest pair switches across the sub-step ($p_{\mathrm{end}} \ne p_{\mathrm{start}}$), the propagation is symmetrized via the half-step pair bridge:
    \[
        |h_{\mathrm{try}}| \le \min\bigl(h(\mathbf{s}_{\mathrm{start}}),\, h(\mathbf{s}_{\mathrm{end}})\bigr), \qquad
        \Phi_{h_{\mathrm{try}}}^{\mathrm{bridge}} \coloneqq \Phi_{h_{\mathrm{try}}/2}^{(p_{\mathrm{end}})} \circ \Phi_{h_{\mathrm{try}}/2}^{(p_{\mathrm{start}})}.
    \]
    This guarantees exact time-reversibility $\|\Phi_{-\Delta t}^{\mathrm{SN}} \circ \Phi_{\Delta t}^{\mathrm{SN}} - \mathrm{Id}\|_2 = 1.31 \times 10^{-15}$ across multi-substep ($N_{\mathrm{sub}} = 8$) deep encounters.
    \item \textbf{Mass-Normalized Hyper-Radius Angular Momentum Pre-Filtering}: In unsoftened dynamics, homothetic central collapse occurs only when the cluster angular momentum vanishes ($L \to 0$). All optimized supernode paths (\texttt{Optimized\-Supernode\-Solver} and \texttt{Jit\-Supernode\-Solver} for all $k \ge 3$) evaluate the scalar angular momentum condition using the mass-normalized hyper-radius $\bar{\rho}_S^2 = \frac{1}{M_S}\sum_{i\in S} m_i \|\mathbf{x}_i\|^2$:
    \begin{equation}
        L^2 < 10^{-20} \left( \bar{\rho}_S^2 + 10^{-24} \right).
        \label{eq:angular_momentum_collapse_filter}
    \end{equation}
    In generic non-collisional motion where $L \gg 10^{-10}$, this scalar guard short-circuits the $\mathcal{O}(k^2)$ central-configuration force evaluation without altering the trajectory.
    \item \textbf{Native Conformal Dirac Momentum Projection}: At the conclusion of the supernode sub-cycle, the momentum conservation constraint $\sum_{i \in S} m_i \mathbf{v}_i = \mathbf{0}$ is enforced via a vectorized linear projection in native code, guaranteeing that center-of-mass momentum drift remains at $1.63 \times 10^{-16}$.
\end{itemize}

\begin{propbox}[Conformal Supernode Micro-Kernel Invariance]{prop:supernode_jit_performance}
Let $\Phi_{\Delta t}^{\mathrm{SN,base}}$, $\Phi_{\Delta t}^{\mathrm{SN,opt}}$, and $\Phi_{\Delta t}^{\mathrm{SN,JIT}}$ denote the synchronous supernode propagators executed by the canonical reference, zero-allocation, and native LLVM JIT micro-kernels over an active cluster $S \subset \{1, \dots, N\}$. Both accelerated micro-kernels preserve exact barycentric momentum closure and multi-substep microscopic time-reversibility:
\begin{equation}
    \Biggl\|\sum_{i \in S} m_i \mathbf{v}_i\Biggr\|_2 = 1.63 \times 10^{-16}, \qquad \bigl\|\bigl(\Phi_{-\Delta t}^{\mathrm{SN}} \circ \Phi_{\Delta t}^{\mathrm{SN}} - \mathrm{Id}\bigr)(\mathbf{s}_S)\bigr\|_2 = 1.31 \times 10^{-15},
\end{equation}
while maintaining exact bit-for-bit agreement ($0.0\text{ ULP}$) with $\Phi_{\Delta t}^{\mathrm{SN,base}}$ across all inherited palindromic stages ($\mathrm{BM}_4^{(7)}$ and $\mathrm{BM}_6^{(11)}$).
\end{propbox}

% ------------------------------------------------------------------------------
\subsubsection{Non-Intrusive Context-Variable Dispatch and Bit-Exact Invariant Verification}
\label{sub:contextvar_modular_dispatch}

A critical architectural mandate of the present monograph is the preservation of the canonical reference implementation as an immutable mathematical baseline. Optimizations must never replace, degrade, or entangle the clean mathematical definitions established in \S\ref{sec:global_pipeline}--\S\ref{sec:atlas_transitions_diagnostics}.

To resolve this software engineering challenge, the library employs the Python \texttt{contextvars} standard library module to achieve thread-local, non-intrusive variant dispatch:
\begin{itemize}
    \item \textbf{Unified Context Tokens in Common Integrator Loops}: In \texttt{integrators/common.py} (and re-exported directly by the standalone paper facade \texttt{06\_\allowbreak High\-Performance\-Optimizations\_\allowbreak Interface.py}), the core propagation primitives (\texttt{step\_drift\_and\_kepler}, \texttt{step\_apex\_kick}, and \texttt{step\_supernode\_flow}) query three thread-local context variables:
    \begin{align}
        \texttt{kepler\_solver} &\gets \texttt{\_CURRENT\_KEPLER\_SOLVER.get(None)}, \\
        \texttt{apex\_evaluator} &\gets \texttt{\_CURRENT\_APEX\_EVALUATOR.get(None)}, \\
        \texttt{supernode\_solver} &\gets \texttt{\_CURRENT\_SUPERNODE\_SOLVER.get(None)}.
    \end{align}
    If a specialized micro-kernel is registered in the current execution context, propagation delegates to the accelerated fast-path; otherwise, the canonical reference solver is executed unconditionally.
    \item \textbf{Transparent Lifecycle Management in \texttt{NBodySystem} and \texttt{BenchmarkRunner}}: Both the top-level orchestrator \texttt{NBodySystem} and the benchmark driver \texttt{BenchmarkRunner} (\texttt{benchmarks/runner.py}) expose the boolean flag \texttt{use\_optimized\_kernels: bool = False} (alongside optional custom solver instances). Keeping \texttt{use\_optimized\_kernels=False} as the default preserves the pure mathematical reference pipeline for invariant verification (and explicitly isolates the reference path by binding \texttt{None} tokens inside a scoped context manager across \texttt{step()}, \texttt{observe\_physical\_state()}, and \texttt{enter\_shadow\_manifold()}).

    Conversely, setting \texttt{use\_optimized\_kernels=True} automatically instantiates and activates all three accelerated micro-kernels (\texttt{Optimized\-Kepler\-Solver}, \texttt{Optimized\-Screened\-Apex\-Potential}, and \texttt{Jit\-Supernode\-Solver}) via \texttt{ContextVar} tokens inside a strict \texttt{try \dots\ finally} block. When operating on 35-digit \texttt{Float128} states (\texttt{is\_f128=True}), the dispatcher automatically bypasses \texttt{Float64} Numba kernels and routes execution through the quadruple-precision reference solvers.
    \item \textbf{Graceful Hardware and Library Fallback}: If Numba or LLVM is absent from the host runtime environment, \texttt{JitSupernodeSolver} falls back to \texttt{Optimized\-Supernode\-Solver} and \texttt{Optimized\-Screened\-Apex\-Potential} falls back to vectorized \texttt{np.add.at} without raising exceptions or interrupting execution.
\end{itemize}

\medskip
\noindent\textbf{Multi-Body Scaling ($N \in [3, 64]$) and Extreme-Regime Opt-vs-Base Equality Certification.}
As synthesized in Figure~\ref{fig:ch06_microkernel_performance_and_invariants} and Table~\ref{tab:optimization_performance_synthesis}, the modular optimization architecture is verified across two complementary evaluation protocols (detailed in Section~\ref{ch:ch07_benchmarks_and_lean4}):
\begin{enumerate}
    \item \textbf{Asymptotic Multi-Body Apex Scaling ($N \in \{3, 4, 8, 16, 32, 64\}$)}: Across synthetic hierarchical $N$-body configurations spanning $3$ to $64$ gravitationally interacting bodies ($3$ to $2{,}016$ pair interactions per stage), the compiled \texttt{@njit(fastmath=False)} pair-accumulation kernel achieves a strictly monotonically increasing speedup of $7.24\times$ at $N=3$ ($40.44\,\mu\mathrm{s} \to 5.59\,\mu\mathrm{s}$), $21.00\times$ at $N=4$ ($39.37\,\mu\mathrm{s} \to 1.87\,\mu\mathrm{s}$), $29.49\times$ at $N=8$ ($60.80\,\mu\mathrm{s} \to 2.06\,\mu\mathrm{s}$), $42.67\times$ at $N=16$ ($107.83\,\mu\mathrm{s} \to 2.53\,\mu\mathrm{s}$), $56.16\times$ at $N=32$ ($266.75\,\mu\mathrm{s} \to 4.75\,\mu\mathrm{s}$), and $62.54\times$ at $N=64$ ($790.08\,\mu\mathrm{s} \to 12.63\,\mu\mathrm{s}$), while preserving exact inter-supernode cross-cluster forces (certified norm $\|\mathbf{F}_{S_1 \leftrightarrow S_2}\|_2 = 0.3600$ on an $N=6$ two-supernode benchmark) and Newton antisymmetry $\|\sum_{i=1}^N \mathbf{F}_i\|_2 = 1.11 \times 10^{-16}$.
    \item \textbf{Direct Opt-vs-Base Equality Proof Across Extreme Singular and Multi-Scale Regimes}: Because the optimized micro-kernels preserve the exact IEEE~754 evaluation order of the canonical reference algorithms (including exact division $\mathbf{v}_0 = \mathbf{w}_0 / \mu$, direct initial Sundman guess evaluation from $(\mathbf{r}_0\cdot\mathbf{r}_0, \mathbf{v}_0\cdot\mathbf{v}_0)$ without apocenter cancellation, incremental canonical momentum updates $\mathbf{w}_{\mathrm{new}} = \mathbf{w}_0 + (\mu\dot{f}\mathbf{r}_0 + \Delta\dot{g}\,\mathbf{w}_0)$, and per-body force accumulation prior to velocity kicks), direct state-by-state comparison between \texttt{use\_optimized\_kernels=True} and \texttt{use\_optimized\_kernels=False} yields bit-exact identity ($0.0\text{ ULP}$, $\|\mathbf{z}_{\mathrm{opt}} - \mathbf{z}_{\mathrm{base}}\|_2 = 0.0$).

    Specifically, exact $0.0\text{ ULP}$ concordance holds across all six extreme stress regimes in Figure~\ref{fig:ch06_microkernel_performance_and_invariants}(c): extreme Kepler pericenter passages down to $1 - e = 10^{-12}$ (\texttt{L1}), exact radial collision bounces at $1 - e = 0, d_0 = 0$ (\texttt{L2}), 8-substep synchronous supernode close encounters (both \texttt{Optimized\-Supernode\-Solver} and \texttt{Jit\-Supernode\-Solver} vs \texttt{Supernode\-Solver} at $0.0\text{ ULP}$), $N=64$ multi-body apex kicks ($2{,}016$ pairs), and full end-to-end trajectories across dynamic Hill chart flips on Comet Lexell 1770 (\texttt{H1}, $r_{\min} = 3.094\,R_J$, zero dyadic saturation at $k_{\max} = 11$) and sub-Roche tidal disruption on Comet Shoemaker--Levy~9 (\texttt{H2}, $r_{\min} = 1.34\,R_J < r_{\mathrm{Roche}}$, zero dyadic saturation at $k_{\max} = 8$).
\end{enumerate}

\begin{table}[!htbp]
\centering
\caption{Architectural, performance, and metrological synthesis across computational micro-kernel tiers under the canonical Voronoi partition ($\gamma = 1.0$).}
\label{tab:optimization_performance_synthesis}
\vspace{0.3em}
\renewcommand{\arraystretch}{1.2}
\resizebox{\textwidth}{!}{%
\begin{tabular}{lccccc}
\toprule
\textbf{Optimization Tier} & \textbf{Phase} & \textbf{Algorithmic Mechanism} & \textbf{Latency ($\mu\mathrm{s}$) \& Speedup} & \textbf{Certified Invariant \& Opt-vs-Base Concordance} \\
\midrule
1. Vectorized / \texttt{@njit} Apex Kick & Phase~3.3 & Cached pairs $\setminus \bigcup (S\times S)$, $\mathbf{G}_{\mathrm{mass}}$, \texttt{@njit} & $5.59\,\mu\mathrm{s}$ ($7.24\times$, $N=3$); $12.63\,\mu\mathrm{s}$ ($62.54\times$, $N=64$) & $\|\sum \mathbf{F}_i\|_2 = 1.11\times 10^{-16}$, $\|\mathbf{F}_{\mathrm{opt}} - \mathbf{F}_{\mathrm{base}}\|_2 / \|\mathbf{F}_{\mathrm{base}}\|_2 = 0.0\text{ ULP}$ ($N \in [3, 64]$) \\
2. Inlined Kepler--Sundman & Phase~3.1 & $2\times 6$ Horner $c_2, c_3$, recurrence \& Battin reuse & $74.82\,\mu\mathrm{s}$ vs $121.11\,\mu\mathrm{s}$ ($1.62\times$, $N=3$) & $|f\dot{g} - \dot{f}g - 1| = 1.12\times 10^{-16}$; Opt vs Base $0.0\text{ ULP}$ ($1-e \ge 10^{-12}$ \& $d_0=0$) \\
3. Zero-Alloc Supernodes & Phase~3.2 & Unrolled $k \le 4$ pairs, $d_{ij}$ reuse \& dyadic stack & Pre-allocated array buffers (fallback tier) & $\|\Delta\mathbf{P}_{\mathrm{cm}}\|_2 = 1.63\times 10^{-16}$; Opt vs Base $0.0\text{ ULP}$ ($8$ sub-steps) \\
4. Native LLVM JIT Supernodes & Phase~3.2 & Numba \texttt{@njit}, $L^2(\bar{\rho}_S)$ filter, dyadic stack & $50.82\,\mu\mathrm{s}$ vs $3040.34\,\mu\mathrm{s}$ ($59.83\times$, 8 sub-steps) & $\|\Phi_{-\Delta t}^{\mathrm{SN}}\circ\Phi_{\Delta t}^{\mathrm{SN}}-\mathrm{Id}\|_2 = 1.31\times 10^{-15}$; JIT vs Base $0.0\text{ ULP}$ \\
5. End-to-End \texttt{ContextVar} Pipeline & All Phases & \texttt{NBodySystem(use\_optimized\_kernels=True)} & Unified thread-local token dispatch ($\mathrm{BM}_6^{(11)}$ / $\mathrm{PBM}_4$) & Lexell \texttt{H1} ($k_{\max}=11$): $0.0\text{ ULP}$; Shoemaker--Levy~9 \texttt{H2} ($k_{\max}=8$): $0.0\text{ ULP}$ \\
\bottomrule
\end{tabular}%
}
\end{table}

\begin{figure}[!htbp]
    \centering
    \includegraphics[width=\textwidth]{chapters/ch06_algorithmic_architecture/figures/fig_ch06_microkernel_performance_and_invariants.png}
    \caption{\textbf{High-performance micro-kernel latency, multi-body scaling ($N \in [3, 64]$), Opt-vs-Base equality proof, and certified invariants ($\gamma = 1.0$).} \textbf{(a)} Per-call latency ($\mu\mathrm{s}$) and speedup factors across the four computational engines ($1.6\times$ to $62.5\times$); \textbf{(b)} multi-body screened apex kick latency and speedup for $N \in \{3, 4, 8, 16, 32, 64\}$ ($7.24\times$ to $62.54\times$); \textbf{(c)} exact numerical equality $\|\mathbf{z}_{\mathrm{opt}} - \mathbf{z}_{\mathrm{base}}\|_2 = 0.0\text{ ULP}$ across six singular and multi-scale benchmark regimes (\texttt{L1}, \texttt{L2}, 8-substep supernode, $N=64$ apex, \texttt{H1}, \texttt{H2}); \textbf{(d)} certified physical and symplectic invariant residuals ($\le 1.31 \times 10^{-15}$, Table~\ref{tab:optimization_performance_synthesis}).}
    \label{fig:ch06_microkernel_performance_and_invariants}
\end{figure}

\pythonrefright{py_bench_gen01}{{benchmarks/run, figures/plot}_generator_01_*.py}{Float64 @njit}
\pythonrefright{py_kernels}{potential/optimized_screened_apex.py, solvers/{optimized_*, jit_supernode}.py}{@njit}