ThreeBody and Rational Dimensions

In the specialized world of aerospace, flight dynamics and mission-control teams propagate multi-body trajectories every day. Cislunar nearshore halo orbits, station-keeping maneuver plans, resonant Jovian moon tours, and low-energy ballistic transfers all live in multi-body gravitational fields where no closed-form general solution exists. When popular accounts call ’the three-body problem’ unsolvable, they conflate a theorem of Heinrich Bruns (1887) and Henri Poincaré (1889) with numerical intractability. Bruns and Poincaré proved that beyond the ten classical algebraic integrals in three-dimensional space (total energy $H$, three components of linear momentum $\mathbf{P}$, three of angular momentum $\mathbf{L}$, and three center-of-mass integrals, leaving seven invariants in the barycentric frame and four in the planar restriction), no further global analytic integrals exist. Operational ground software integrates the equations routinely, provided the engineer is honest about force-model assumptions, close-approach singularities, step truncation, and finite-precision rounding.

Keep that last condition in mind: honesty about models, singularities, and rounding is a central takeaway of this piece.

At Braidpoint, we started the ThreeBody project as a concrete astrodynamics and compiler verification testbed for the Clef programming language and Composer compiler. With that, the ThreeBody project and the rational measure-exponent extension described below remain active design work. Writing the design down at full mathematical resolution forces our numerical, physical, and type-level commitments to meet on the same equations before we lay claim to any result.

This is a deep cut, as they say, written with two readers in mind: generalist software engineers curious how type systems and custom silicon handle nasty numerical edge cases, and specialists who wrestle with complex domains such as astrodynamics. Working through the hardest cases in the compiler and their result at runtime is how we make demanding numeric domains approachable, if not routine, across a wide variety of high-precision workloads. For the next several sections, the treatment is unapologetically mathematical with a few signposts along the way.

To begin, we consider $N$ point masses with positions $\mathbf{q}_i \in \mathbb{R}^d$ ($d \in \{2, 3\}$), momenta $\mathbf{p}_i = m_i \dot{\mathbf{q}}_i$, pairwise separations $\mathbf{d}_{ij} = \mathbf{q}_j - \mathbf{q}_i$, and Euclidean distances $r_{ij} = \|\mathbf{d}_{ij}\|$. Their motion is governed by the Newtonian Hamiltonian $H(\mathbf{q}, \mathbf{p})$ and the exact pairwise gravitational force $\mathbf{F}_i(\mathbf{q})$ on body $i$.

Our baseline ThreeBody configuration places three equal masses at the vertices of a right triangle in a two-dimensional plane (a symmetric-mass cousin of Carl Burrau’s 3-4-5 free-fall problem integrated numerically by Victor Szebehely and C. Frederick Peters in 1967) and releases them from rest, so $\mathbf{P}(0) = \mathbf{0}$ and $\mathbf{L}(0) = \mathbf{0}$. Choosing characteristic scales $L_0$, $M_0$, and dynamical time $T_0 = \sqrt{L_0^3 / (G M_0)}$ with momentum scale $P_0 = M_0 L_0 / T_0$ normalizes the dimensionless coupling to $\tilde{G} = 1$ while preserving explicit physical units at every interface. Repeated close passages in that free-fall sequence drive separations down by orders of magnitude before ejecting a single body and leaving a tight binary, exposing every shortcut in the force kernel and every lost bit in the arithmetic.

Anyone who has watched 3 Body Problem on Netflix will recognize dramatized variants of this chaotic dance. Three masses released from rest on a right triangle give us the clean laboratory version.

Softening: Vlasov Continua & Point Masses

Many published $N$-body examples replace the Newtonian $r_{ij}^{-3}$ denominator with a Plummer-softened potential $\Phi_{\varepsilon}$ and force kernel $\mathbf{F}_{ij}^{\text{soft}}$ controlled by a length parameter $\varepsilon > 0$.

Whether that modification is sound physics or a silent error depends on what a simulation particle represents. The two regimes split along physical scale.

In galactic and cosmological codes (GADGET, AREPO, GIZMO), discrete bodies do not represent individual stars or spacecraft. They are Monte Carlo macro-particles sampling a smooth six-dimensional phase-space mass density $f(\mathbf{q}, \mathbf{v}, t)$ governed by the collisionless Vlasov-Poisson equations.

When one computational particle carries $10^6$ solar masses of dark matter or stellar mass, hard two-body scattering between neighboring tracer particles is referred to as a discretization artifact of coarse phase-space sampling. Gravitational softening (or finite-support spline kernels, as analyzed by Walter Dehnen in 2001 and Chris Power et al. in 2003) suppresses artificial two-body relaxation so the discrete potential tracks the continuum mean field. Wayne Hayes demonstrated in 2003 through shadowing integrations that under-softening a collisionless ensemble below the inter-particle scale degrades numerical shadowability. In collisionless continuum astrophysics, softening is physically required.

Point-mass astrodynamics solves the opposite physical problem. Three discrete gravitating bodies in free fall, a spacecraft executing a powered perilune pass, or a resonant asteroid triplet are genuine point-mass or small-body systems rather than tracer samples of a fluid. Close encounters are the dynamics under study. Injecting $\varepsilon > 0$ integrates a perturbed Plummer Hamiltonian $H_{\varepsilon} = T + \sum_{i<j}\Phi_{\varepsilon}(r_{ij})$ instead of Newtonian gravity. Inside $r_{ij} \ll \varepsilon$, the true $r_{ij}^{-2}$ force magnification collapses into a linear harmonic core.

That core caps peak accelerations, alters deflection angles during flybys, and damps the local variational stretching that stresses finite-precision accumulators.

Hard-coding rSq + 1e-12 in source code also corrupts compiler verification. Because $r_{ij}^2 + \varepsilon^2 \ge \varepsilon^2 > 0$ holds everywhere on $\mathbb{R}^{dN}$, static interval propagation sees a strictly positive denominator and clears the expression without recording a singularity obligation in what emerges as our program semantic graph.

The pole at $r_{ij} = 0$ disappears behind an unverified constant. In general software terms, slipping + 1e-12 into a gravitational denominator is the numerical equivalent of catching a divide-by-zero exception and returning a plausible-looking fallback value. In mission-control trajectory code, it quietly bends a close perilune or planetary flyby onto the wrong outbound asymptote without raising a single diagnostic. Addressing that gap cleanly, with both precision and execution efficiency intact, motivates how we split simulation regimes below.

Separated Regimes and Collisional Trajectories

Because no single numerical formula is valid across smooth continuum sampling, separated orbital flight, singular point-mass flybys, and physical surface impacts, our simulation plan separates non-colliding and colliding regimes explicitly as opposed to simply hiding both behind softening.

Non-Colliding Point Masses and Error-Bounded Far Fields

For non-colliding point-mass runs (Phase 1 of ThreeBody), we evaluate the exact unsoftened $1/r_{ij}^3$ kernel over an explicit minimum-separation domain $\mathcal{D}(r_{\min})$.

Over unconstrained reals, the range of $r_{ij}^{-3}$ is unbounded above, and Composer’s planned range checker under Numeric Selection must refuse fixed-scale accumulator commitment until a lower separation bound is supplied. That bound can come from checked physical separation invariants (for example, bounded Hill-stability or orbital clearance premises) or from an explicit runtime domain guard on every force evaluation. If a trajectory crosses $r_{ij} < r_{\min}$, the step halts with a located domain violation or hands off to a close-encounter regime instead of inventing a denominator floor.

Scaling non-colliding trajectories to large $N$ follows the same discipline. Rather than Barnes-Hut octrees parameterized by heuristic opening angles $\theta = s/R$, our multi-body scaling design uses the Fast Multipole Method (Leslie Greengard and Vladimir Rokhlin, 1987, with adaptive tree refinement following J. Carrier, Greengard, and Rokhlin, 1988). For a source cell of total mass $M_{\text{cell}}$ enclosed in a sphere of radius $a$ and a well-separated target at center distance $R > a$, truncating the three-dimensional spherical-harmonic expansion at order $p$ bounds the potential remainder analytically.

An accuracy target chooses $p$ directly from the geometry ratio $a/R$. Well-separated multipole-to-local translations (FMM List V and adaptive Lists W and X) lower cleanly onto wide SIMD, GPU, or NPU tensor tiles, while FMM’s near-field List U (each leaf box and its immediate neighbors where $R/a$ violates separation) structurally isolates close pairs for direct unsoftened summation. Planar trajectories embedded in physical three-space retain the three-dimensional $1/r^2$ Newtonian force and must use the restricted three-dimensional multipole kernel rather than the two-dimensional logarithmic Laplace potential $\ln r$, which governs infinite line masses. Where a user genuinely intends collisionless Vlasov continuum smoothing, Plummer or spline potentials remain available as distinct force laws with their own typed parameter $\varepsilon : \text{float}\langle\text{m}\rangle$.

Software generalists can read FMM’s List V and List U split as a classic systems pattern: batch the predictable fast path onto vector hardware, and isolate dangerous edge cases in an exact handler. Ground-system engineers will recognize the same architecture from space situational awareness and conjunction screening, where bulk catalog propagation filters thousands of resident space objects down to the few close-approach pairs requiring high-fidelity evaluation.

Analytic Regularization and Finite-Radius Contact

A separation guard $r_{ij} \ge r_{\min}$ keeps a non-colliding integrator honest, yet stopping at $r_{\min}$ cannot propagate through a deep perilune, a resonant three-body scramble, or an actual collision. Supporting colliding and near-colliding simulations requires distinguishing two physical cases: singular point-mass encounters ($R_i = 0, r_{ij} \to 0$) and finite-radius physical contact ($r_{ij} = R_i + R_j > 0$).

For point-mass close encounters and mathematical collisions ($r_{ij} \to 0$), celestial mechanics removes the Keplerian pole analytically without altering Newton’s law. If you have never worked with regularized astrodynamics, the core idea is surprisingly visual: instead of blunting gravity near zero, we re-clock the integration so physical time steps shrink automatically in proportion to distance ($dt = r\,ds$), while lifting coordinates into a squared space that unfolds a sharp $180^\circ$ gravitational slingshot into a smooth pass straight through the origin. Tullio Levi-Civita solved the planar case in 1920 via complex squaring, and Paul Kustaanheimo and Eduard Stiefel extended the construction to three-dimensional space in 1965 by lifting relative Cartesian coordinates $\mathbf{q} = \mathbf{q}_2 - \mathbf{q}_1 \in \mathbb{R}^3$ (embedded in $\mathbb{R}^4$ with fourth component zero) to a four-parameter vector $\mathbf{u} \in \mathbb{R}^4$ paired with Karl Sundman’s fictitious time $s$.

The Kustaanheimo-Stiefel (KS) linear operator $\mathcal{L}(\mathbf{u})$ takes the $4 \times 4$ matrix form below.

This matrix satisfies $\mathcal{L}(\mathbf{u})^T \mathcal{L}(\mathbf{u}) = r I_4$. Along trajectories obeying the bilinear gauge constraint $u_4 u_1' - u_3 u_2' + u_2 u_3' - u_1 u_4' = 0$ (which reduces in the planar $u_3 = u_4 = 0$ slice to Levi-Civita’s $2 \times 2$ matrix), it maps differentials via $d\mathbf{q} = 2\mathcal{L}(\mathbf{u})\,d\mathbf{u}$.

Write $(\cdot)' = d/ds$, let $\mu = G(m_1 + m_2)$, and define the pair’s osculating specific two-body energy $h_2 = \frac{1}{2}\|\dot{\mathbf{q}}\|^2 - \frac{\mu}{r} = (2\|\mathbf{u}'\|^2 - \mu)/r$. The vector $\mathbf{P}$ collects the differential tidal acceleration exerted on the pair by all external bodies $k \ge 3$.

Substituting $\mathbf{q} = \mathcal{L}(\mathbf{u})\mathbf{u}$ and $dt = r\,ds$ into the singular relative equation $\ddot{\mathbf{q}} = -\mu\mathbf{q}/r^3 + \mathbf{P}$ cancels the $r^{-3}$ pole algebraically and yields the regular KS equations of motion.

For a bound pair ($h_2 < 0$) in isolation ($\mathbf{P} = \mathbf{0}$), the equation is a linear harmonic oscillator with constant real frequency $\omega = \sqrt{-h_2/2}$. Its solution passes smoothly right through a head-on collision at $\mathbf{u} = \mathbf{0}$ ($r = 0$) with finite derivatives in $s$. In our planned multi-body architecture (following Sverre Aarseth’s NBODY6/NBODY7 neighbor-switching strategy and three-body KS extensions such as Elena Lega and Massimiliano Guzzo, 2024), actor-scheduled neighbor tracking switches pairs crossing $r_{ij} < r_{\text{reg}}$ into regularized coordinates and restores Cartesian coordinates once they separate past hysteresis radius $r_{\text{exit}} > r_{\text{reg}}$.

Finite-radius bodies ($R_i > 0$), such as spacecraft hard-body volumes in conjunction assessment, orbital debris fragments, or accreting planetesimals, follow a different collisional branch because physical contact or keep-out violation occurs at $r_{ij} = R_i + R_j > 0$ before the point-mass pole is ever reached. Within a step $[t_n, t_{n+1}]$, dense output and root location isolate the earliest contact epoch $t_c$ where the signed clearance $g_{ij}(t) = \|\mathbf{q}_j(t) - \mathbf{q}_i(t)\| - (R_i + R_j)$ crosses zero with approaching normal speed $\dot{g}_{ij}(t_c) < 0$. At $t_c$, state transition rules resolve either a mission-level conjunction event, inelastic accretion ($M = m_i + m_j$, $\mathbf{p}_M = \mathbf{p}_i + \mathbf{p}_j$ at the pair barycenter), or normal impulse exchange $J$ with restitution coefficient $e \in [0, 1]$ along $\hat{\mathbf{n}} = (\mathbf{q}_j - \mathbf{q}_i)/r_{ij}$.

Any partially inelastic impact ($0 \le e < 1$) or merger dissipates mechanical kinetic energy $\Delta T$.

That energy loss breaks time-reversal invariance physically rather than numerically. In our verification design, live time-reversal diagnostics apply across conservative segments (or strictly elastic $e = 1$ event-resolved rebounds), while dissipative collisions verify conservation of total linear momentum $\sum_i \mathbf{p}_i$, angular momentum accounting for spin transfer, and exact energy-decrement bookkeeping.

Rational Dimensions & Astrodynamics

Look closely at the physical dimensions inside the regularized KS equations above. In unregularized Cartesian coordinates, computing Euclidean distance sqrt (dx*dx + dy*dy + dz*dz) maps $\mathsf{L}^2$ back to $\mathsf{L}$, staying entirely within integer unit exponents. The moment a trajectory switches into Levi-Civita or Kustaanheimo-Stiefel coordinates, fractional dimensions enter the state vector itself.

Over base dimensions of length $\mathsf{L}$, mass $\mathsf{M}$, and time $\mathsf{T}$, the identity $r = \|\mathbf{u}\|^2$ with $\dim(r) = \mathsf{L}$ assigns both the KS state vector $\mathbf{u}$ and the entries of $\mathcal{L}(\mathbf{u})$ the fractional dimension $\mathsf{L}^{1/2}$. Sundman’s differential relation $dt = r\,ds$ assigns fictitious time $s$ dimension $\mathsf{T}\mathsf{L}^{-1}$. Differentiating with respect to $s$ gives regularized velocity $\mathbf{u}' = d\mathbf{u}/ds$ dimension $\mathsf{L}^{3/2}\mathsf{T}^{-1}$ and regularized acceleration $\mathbf{u}'' = d^2\mathbf{u}/ds^2$ dimension $\mathsf{L}^{5/2}\mathsf{T}^{-2}$. Every term of the perturbed KS oscillator balances cleanly in $\mathbb{Q}$-exponents.

Fractional dimensions show up across everyday astrodynamics even before regularization:

• Gravitational root factors and vis-viva scaling: Standard gravitational parameter $\mu = GM$ has dimension $\mathsf{L}^3\mathsf{T}^{-2}$, so $\sqrt{\mu}$ carries $\mathsf{L}^{3/2}\mathsf{T}^{-1}$. Specific angular momentum $h = \sqrt{\mu p}$ multiplies $\mathsf{L}^{3/2}\mathsf{T}^{-1}$ by semi-latus rectum root $p^{1/2}$ ($\mathsf{L}^{1/2}$) to land on $\mathsf{L}^2\mathsf{T}^{-1}$.
• Keplerian period and mean motion: Orbital period $T = 2\pi a^{3/2}\mu^{-1/2}$ divides $\mathsf{L}^{3/2}$ by $\mathsf{L}^{3/2}\mathsf{T}^{-1}$ to yield $\mathsf{T}$, while mean motion $n = \mu^{1/2}a^{-3/2}$ yields $\mathsf{T}^{-1}$.
• Universal variables: Battin’s universal anomaly $\chi$, defined through $\sqrt{\mu}\,dt = r\,d\chi$, carries dimension $\mathsf{L}^{1/2}$, with Stumpff-series terms scaling in half-integer powers of length and energy.
• Natural-unit nondimensionalization: Forming dynamical time $T_0 = L_0^{3/2}(G M_0)^{-1/2}$ and momentum scale $P_0 = M_0^{3/2} G^{1/2} L_0^{-1/2}$ from separate typed factors steps directly through half-integer exponents of $\mathsf{L}$ and $\mathsf{M}$.

Calculating the Costs of Bare Floats

After wading through spherical harmonics and four-dimensional coordinate lifts, this is a good place to come up for air and look at what all of those fractional powers mean for someone actually writing and maintaining software.

Anyone who has worked on mission software stacks knows where unit bugs bite in production. Rarely does someone design to add meters directly to seconds inside a five-line helper. Instead, mismatches hide at subsystem handoffs where units were erased to bare double[] state vectors: $\mu$ passed in $\text{km}^3/\text{s}^2$ into a maneuver targeter expecting $\text{m}^3/\text{s}^2$, or auxiliary universal-variable and KS state arrays whose fractional dimensions could not be expressed in the host language’s unit library.

Software generalists run into the same wall in language ecosystems like C#, C++ dimensional template libraries, or Rust unit crates. In Andrew Kennedy’s classical units-of-measure calculus (and Clef’s current Units of Measure specification), $k$ primitive base measures generate the free abelian group $\mathbb{Z}^k$. Every measured type carries an integer exponent vector $\mathbf{z} \in \mathbb{Z}^k$, and sqrt : F<'U^2> -> F<'U> requires every unit exponent under the radical to be divisible by two. The moment a formula takes the square root of an odd power, such as $\mu : \text{float}\langle\text{m}^3/\text{s}^2\rangle$ or separation $r : \text{float}\langle\text{m}\rangle$, integer divisibility fails.

When the type checker rejects a valid physical intermediate, developers do the only thing an integer-exponent language lets them do: strip the units off, run the math on raw floats, and re-attach the expected unit at the end by assertion. Consider computing Keplerian orbital period from semi-major axis a and gravitational parameter mu:

[<Measure>] type m
[<Measure>] type s

// Today under integer (Z) exponents: strip units, trust yourself in the middle
let orbitalPeriodStripped (a: float<m>) (mu: float<m^3 / s^2>) : float<s> =
    let rawA = float a
    let rawMu = float mu
    let rawT = 2.0 * Math.PI * (rawA ** 1.5) / sqrt rawMu
    rawT * 1.0<s>

// Proposed rational (Q) exponents: units stay attached through every intermediate
let orbitalPeriod (a: float<m>) (mu: float<m^3 / s^2>) : float<s> =
    let rootA3 = pow_q<3, 2> a     // inferred: float<m^(3/2)>
    let rootMu = sqrt mu           // inferred: float<m^(3/2) / s>
    2.0 * Math.PI * (rootA3 / rootMu)

Look at the failure mode hidden inside orbitalPeriodStripped. If a tired engineer accidentally writes sqrt rawMu / (rawA ** 1.5) on a Friday afternoon, multiplying by 1.0<s> stamps seconds onto a value that is actually inverse seconds (1/s), and the compiler waves it through. In our proposed $\mathbb{Q}$-exponent version, flipping rootA3 / rootMu produces float<s^-1>, which fails immediately against the declared float<s> return type without a single annotation on rootA3 or rootMu.

Why ℚ Calms the Compiler

Choosing which number system governs unit exponents is an engineering trade-off between expressiveness and what the compiler can decide automatically without pestering the programmer.

Integers ($\mathbb{Z}$) handle multiplication, division, and reciprocals (s^-1), yet break on odd roots. Jumping all the way to arbitrary real exponents ($\mathbb{R}$) would allow nonsensical units like $\text{m}^{\pi}$ or $\text{s}^{\sqrt{2}}$, pulling irrational equality checks into type inference for powers that never appear in dimensional physics. The rational numbers ($\mathbb{Q}$), ordinary fractions $a/d$ of integers reduced by greatest common divisor, sit in the sweet spot: closed under every root and fractional power we encountered above, yet represented inside the compiler as exact integer pairs with zero floating-point rounding of their own.

Extending measure exponents from $\mathbb{Z}$ to $\mathbb{Q}$ upgrades the measure sort from a $\mathbb{Z}$-module ($\mathbb{Z}^k$) to a finite-dimensional vector space $\mathbb{Q}^k$ over the rationals.

Here is the pleasant surprise for language tooling: moving from $\mathbb{Z}^k$ to $\mathbb{Q}^k$ makes unification algebraically cleaner. Over integers, solving $2\mathbf{z} = \mathbf{w}$ requires Hermite-normal-form divisibility checks because $2$ has no multiplicative inverse in $\mathbb{Z}$. Over the field $\mathbb{Q}$, every nonzero fraction $p/d$ is invertible. Under multiplication ($\mathbf{q} + \mathbf{r}$), division ($\mathbf{q} - \mathbf{r}$), and rational power primitives pow_q<p, d> : F<'U> -> F<'U^(p/d)> ($\frac{p}{d}\mathbf{q}$), dimensional checking reduces to standard rational Gaussian elimination (reduced row echelon form over exact normalized fractions). Unification stays decidable, unitary, and principal: every satisfiable set of equations has a single most-general measure solution up to fresh $\mathbb{Q}$-variables.

We are putting real care into how that algebra runs across design time and build time so developers get stronger guarantees without extra ceremony:

• At design time (instant inference, zero annotation tax): Because $\mathbb{Q}^k$ unification has principal solutions, Clef’s checker (CCS) can infer fractional units on local bindings, closures, and helper functions automatically. Developers write ordinary code without sprinkling unit casts or manual proof hints on intermediate steps. When a formula disagrees dimensionally, editor diagnostics report normalized human-readable exponents directly on the offending expression rather than dumping raw solver matrices.
• At build time (Baker saturation and zero-cost codegen): Most unit libraries erase measures right after parsing and throw that knowledge away. Clef’s Units of Measure architecture instead carries resolved measures through Baker elaboration into the Program Semantic Graph (Fidelity.PSG), where downstream range analysis, cvc5 domain checks, and accumulator-adequacy obligations can consume them directly without asking the programmer to restate invariants in a sidecar spec language. Once those build-time obligations are discharged, lowering strips the measure metadata completely: CPU, GPU, NPU, and FPGA backends emit flat hardware registers and vector instructions with zero wrapper objects, zero boxing, and zero runtime overhead.

You get rich fractional units in your editor while you type, automated invariant plumbing during compilation, and bare-metal register speed when the binary runs.

Decidable Obligations and Numeric Selection

A reader versed in formal methods might reasonably ask where the hidden proof burden sits. In general dependent type systems, encoding value predicates inside types makes type inference and arbitrary proof search undecidable, forcing engineers to supply manual proof terms, fuel limits, or solver tactics (The Decidability Sweet Spot). How do we verify an unsoftened gravitational kernel and select hardware representations without reintroducing that tax?

In our Decidable By Construction architecture, we avoid a monolithic dependent-type solver by stratifying verification into distinct stages, each confined to a decidable or terminating mathematical fragment. Insisting on dimensional types and Numeric Selection in Clef keeps those stages connected. We can imagine how those pieces meet in a Clef pairwise encounter check.

[<Measure>] type kg
[<Measure>] type N = kg * m / s^2

type Vec3<[<Measure>] 'u> = { X: float<'u>; Y: float<'u>; Z: float<'u> }

type EncounterRegime =
    | SeparatedForce of Vec3<N>
    | SwitchToRegularized of uNorm: float<m^(1/2)>

// Unsoftened point-mass pair check: no "+ 1e-12" hiding the pole.
// Measured quantities trigger range-coverage checks, and the r >= rMin
// guard bounds 1 / r^3 away from zero before representation selection.
let pairGravity
    (rMin: float<m>)
    (m1: float<kg>)
    (m2: float<kg>)
    (d: Vec3<m>) : EncounterRegime =
    let g = 6.67430e-11<m^3 kg^-1 s^-2>
    let rSq = d.X * d.X + d.Y * d.Y + d.Z * d.Z     // inferred: float<m^2>
    let r = sqrt rSq                                 // inferred: float<m>
    if r >= rMin then
        let scale = (g * m1 * m2) / (r * rSq)        // inferred: float<N / m>
        SeparatedForce { X = scale * d.X; Y = scale * d.Y; Z = scale * d.Z }
    else
        SwitchToRegularized (sqrt r)                 // inferred: float<m^(1/2)>

Notice that no concrete bitwidth or posit format (f64, bposit32) appears anywhere in pairGravity. Under our Numeric Selection specification, float<'u> denotes an abstract physical real quantity in source code, while concrete hardware encodings belong to Composer’s target catalog and platform boundary declarations. Three bounded compiler mechanisms turn pairGravity into checked machine code without handwritten proofs:

1. Measured types trigger the range requirement: Under Numeric Selection, an unmeasured float with an unobservable range may default to IEEE f64, whereas any quantity carrying a physical measure (float<m>, float<kg>, Vec3<N>) requires a bounded interval [a, b] before Composer commits a hardware format. Insisting on dimensional types prevents unbounded physical intermediates from slipping silently into fixed-point or quire-backed arithmetic, while CCS settles all unit exponents in pairGravity (including float<m^(1/2)> on sqrt r) via rational Gaussian elimination.
2. Branch guards and closed-form intervals bound the pole: Given caller or platform preconditions establishing rMin > 0.0<m> alongside bounded domain extents and masses, the if r >= rMin guard supplies r * rSq >= rMin^3 > 0. Inside that branch, cvc5 discharges the non-zero denominator obligation and Composer propagates outward-rounded interval images across closed-form expressions within decidable quantifier-free real fragments (QF_LRA and bounded polynomial intervals), avoiding undecidable trigonometric proof search. When r < rMin, execution branches to SwitchToRegularized (sqrt r), passing the pair to the pole-free KS oscillator before 1 / (r * rSq) is evaluated.
3. Representation and quire commitment is a finite argmin: Once every measured intermediate in the SeparatedForce branch carries a closed interval [a, b], selecting a concrete representation requires no open-ended theorem proving. Composer filters the target’s finite format catalog (f64, fixed-point, or bposit32 paired with its 800-bit quire) to formats whose dynamic range covers [a, b], checks that products and partial sums fit within the 800-bit signed quire capacity, and picks the lowest-cost admissible format under a ULP-floored error metric.

Keeping adjacent parts of our type system distinct prevents confusion as that pipeline comes together:

• Reciprocal units such as $\text{s}^{-1}$ already belong to $\mathbb{Z}^k$ and need no extension.
• Rational measure exponents ($\mathbb{Q}^k$) check dimensional homogeneity across fractional roots, KS coordinates, and universal-variable formulas, carrying measured types through Baker saturation into the program semantic graph (Proofs from Dimensional Types).
• Negative and fractional types, discussed in A Path Less Traveled, type operational resource duals and directed channels as opposed to physical measure exponents or numerical reversibility.

The Lyapunov Horizon

Once the force law and coordinate domain are fixed, numerical error enters a gravitational propagation at three distinct points, and no single arithmetic representation eliminates all three.

First comes pairwise coordinate subtraction $\mathbf{d}_{ij} = \mathbf{q}_j - \mathbf{q}_i$. By Sterbenz’s lemma, when two IEEE floating-point coordinates satisfy $\frac{1}{2} \le q_{j,\alpha} / q_{i,\alpha} \le 2$, their subtraction introduces zero new rounding error on the stored inputs. Cancellation of leading bits during a close pass nonetheless exposes relative errors already accumulated in $\mathbf{q}_i$ and $\mathbf{q}_j$ across thousands of preceding steps.

Neither compensated summation nor a wide posit quire can reconstruct coordinate precision that vanished before the subtraction began.

Second comes non-associative reduction across squared norms $\mathbf{d}_{ij} \cdot \mathbf{d}_{ij}$, multi-body force sums $\sum_{j \ne i} \mathbf{F}_{ij}$, and energy or momentum tallies. In IEEE 754 binary64 ($53$ significand bits, unit roundoff $u = 2^{-53} \approx 1.11 \times 10^{-16}$), rounding after every addition breaks associativity.

Reordering force contributions across threads or tree nodes changes low-order bits on every step. In distributed ground-control pipelines, that non-associativity is why two compute nodes replaying the same initial ephemeris with different thread schedules can diverge on downstream maneuver targets or close-approach miss distances.

Third comes state writeback: even if an inner force sum or dot product is accumulated exactly, updating $\mathbf{q}_{i,n+1}$ and $\mathbf{p}_{i,n+1}$ rounds back to the working state width at the step boundary.

To separate intermediate-reduction error from state-width error rather than conflating them, our arithmetic construction and placement design specifies four paths on identical initial conditions and step sizes:

• Uncompensated IEEE binary64: standard double-precision evaluation with per-operation round-to-nearest-even semantics.
• Compensated IEEE binary64 (TwoSum / Kahan-Babuška-Neumaier): double-precision state storage paired with error-free residual tracking across force, norm, and diagnostic sums.
• Bounded posit bposit32 (eS = 5, rS = 6) with an 800-bit quire: 32-bit bounded-regime posit state storage paired with hardware exact accumulation on an Arty A7-100T FPGA sidecar over raw Layer 2 Ethernet via an eBPF interlink.
• Convergence-checked high-precision reference: independent integration in MPFR / arb ball arithmetic (or float128 under step-halving convergence checks) providing reference states $(\mathbf{q}_i^{\text{ref}}(t), \mathbf{p}_i^{\text{ref}}(t))$.

Compensated binary64 serves as a mandatory experimental control. Donald Knuth and Ole Møller’s error-free transformation (analyzed for compensated dot products by Takeshi Ogita, Siegfried Rump, and Shin’ichi Oishi in 2005) captures the exact rounding residual of IEEE addition inside ordinary registers.

Measuring compensated binary64 alongside uncompensated binary64 isolates how much trajectory drift disappears on commodity IEEE hardware before we attribute any improvement to custom silicon.

On the FPGA path, our ThreeBody plan uses the bounded-posit design, which caps regime width at six bits for parallel multiplexor decoding and specifies an 800-bit fixed-point quire ($25 \times 32$-bit words, distinct from the Posit Standard 2022 $16n = 512$-bit quire for posit32). Under our posit arithmetic and Numeric Selection quire-adequacy contract, an accumulator carrier $Q$ with real denotation $D : Q \to \mathbb{R}$ obeys four exact monoid and finalization laws.

Those laws hold whenever every product $x_k y_k$ is exact at the quire’s fixed-point quantum and every reachable partial sum and parallel reduction merge stays within the 800-bit signed capacity. Posits concentrate fraction bits near magnitude $1$ and taper toward lower precision near $0$ and $\infty$, which is why scaling ThreeBody to natural units ($\tilde{G} = 1$) places characteristic separations and forces inside the format’s peak-significance band.

For conservative segments, our baseline time-stepper is the fixed-step symmetric Störmer-Verlet (leapfrog) map $\Phi_{\Delta t}$.

Two geometric properties of $\Phi_{\Delta t}$ in exact real arithmetic do different jobs here (Rounding on Real Hardware). Symplecticity ($\Phi_{\Delta t}^* \omega = \omega$) bounds secular drift in total energy over exponentially long intervals via backward error analysis on a nearby shadow Hamiltonian, provided $\Delta t$ stays uniform or adapts through a time-reversible transformation such as $dt = r\,ds$. Time-symmetry under momentum negation $S(\mathbf{q}, \mathbf{p}) = (\mathbf{q}, -\mathbf{p})$ gives an exact step inverse.

Propagating $K$ steps forward to $T = K\Delta t$, negating momenta via $S$ (an ordinary value negation -p preserving float<kg m / s>), propagating $K$ steps forward with the same $\Delta t$, and applying $S$ again retraces the trajectory back toward $(\mathbf{q}(0), \mathbf{p}(0))$ live, without recording a checkpoint tape. Mission-control engineers perform a close cousin of this maneuver during orbit determination, back-propagating state estimates across tracking passes to tie new range and Doppler residuals to an earlier epoch. Finite-precision rounding breaks exact involution closure ($S \circ \tilde{\Phi}_{\Delta t} \circ S \ne \tilde{\Phi}_{\Delta t}^{-1}$), turning the dimensionless return residual $R_{\text{rev}}(T)$ into a direct sensor of cumulative irreversible rounding.

Near close approaches in the three-body problem, trajectories diverge exponentially at a rate governed by the continuous system’s maximal Lyapunov exponent $\lambda > 0$ (The Lyapunov Window). Changing arithmetic formats never changes $\lambda$. If each step amplifies existing phase-space error by a factor $A \approx e^{\lambda \Delta t} > 1$ and injects local truncation plus rounding error bounded by $\rho$, step error grows geometrically.

Against a dimensionless tolerance $\delta > 0$, the useful prediction horizon $T_{\delta}$ scales logarithmically with perturbation amplitude.

Each additional bit of effective precision extends $T_{\delta}$ additively by roughly $(\ln 2)/\lambda$ dynamical times rather than defeating chaos. We measure that horizon directly against the high-precision reference trajectory rather than inferring accuracy from conserved integrals.

Tracking barycentric momentum drift $\|\sum_i \mathbf{p}_i(t)\| / P_0$, angular momentum drift, and energy residual $|H(t) - H(0)| / |H(0)|$ catches non-symmetric rounding bias, yet flight-dynamics teams know the classic operational trap: along-track (in-track) phase drift. An integrator with slight period error can keep energy and momentum plots looking flat while placing a spacecraft kilometers ahead or behind in its orbital slot. Reference comparison $E(t)$, tape-free return residual $R_{\text{rev}}(T)$, and invariant drift stay separate curves in the planned ThreeBody experiments.

Once we have the foundational measurements and software tooling that demonstrate this in a 2D plane, the next steps include expanding the count, size, and placement of bodies and moving the simulation into three dimensions. That step will trigger its own round of follow-on work to keep live visualization smooth and responsive across dramatic scale changes, placing coordinated demand on GPU, NPU, and FPGA accelerator resources that we look forward to exploring, for this problem and for many others.

Physical Fidelity Without Stripping the Types

Stepping back from the equations, the engineering goal behind all of this machinery is straightforward. Teams building ground-control and mission-control software, from ephemeris propagation and orbit determination to maneuver planning and conjunction screening, should be able to keep physical laws, singularity boundaries, accumulator guarantees, and dimensional units attached to executable code through the hardest parts of the math. Stripping units down to raw double arrays at a Kustaanheimo-Stiefel transformation or hiding a Newtonian pole behind an ad hoc softening constant trades away the exact invariants that catch silent operational bugs.

ThreeBody gives us a concrete astrodynamics laboratory where every piece of that commitment has to hold together at once: unsoftened point-mass domains alongside FMM far-field error bounds, KS regularization and event-resolved contact for colliding trajectories, four controlled arithmetic paths including compensated binary64 and FPGA bposit32 quires, and $\mathbb{Q}$-dimensional unification across fractional powers. As each component moves from specification into running Clef and Composer code, its curves and diagnostics will have to stand on measured reference residuals.

Three-body astrodynamics is a highly specialized test case, yet that reach holds practical lessons across a wide array of commercial and industrial domains. None of the compiler mechanisms stressed here stay locked inside spaceflight software. Once rational dimensions, range-checked denominators, and verified accumulator placement exist in the toolchain, every Clef program inherits them:

• Autonomous robotics, guidance, and sensor fusion: Inertial navigation and Kalman filtering routinely combine coordinate-frame rotations with fractional noise spectral densities such as $\text{m}/\sqrt{\text{Hz}}$ and angle random walk in $\text{rad}/\sqrt{\text{s}}$, keeping units intact across state estimators where frame and scaling bugs otherwise hide in untyped matrices.
• Quantitative finance, risk, and settlement: Multi-leg derivatives pricing, volatility surfaces scaling with $\sqrt{t}$, and parallel portfolio reductions face the same twin hazards ThreeBody isolates: non-associative rounding drift across millions of sums and unit mismatches across currencies, rates, and tenors.
• Energy grids, industrial control, and multiphysics digital twins: Power-system state estimation, battery electrochemistry (where diffusive impedance scales with $t^{1/2}$), and coupled fluid-structure solvers mix stiff local events with bulk field calculations across heterogeneous processors.
• AI, quantum-adjacent simulation, and everyday software: Narrow-bitwidth neural and neuromorphic inference on NPUs and FPGAs depends on wide-accumulator discipline so dot products do not collapse into rounding noise, quantum circuit simulations lean on controlled phase precision, and everyday backend or embedded tooling gets zero-cost dimensional safety, supervised actor concurrency, and clean hardware targeting out of the same compiler.

Pushing against an unforgiving physical problem forces us to solve precision, typing, and hardware coordination cleanly in the Fidelity Framework. We consider this to be the honest work of advanced design made practical, and it stands in contrast to the old approach of simply leaving application developers to patch around these concerns by hand. We are excited to see ThreeBody take shape, both for what its trajectories will measure directly and for the everyday engineering confidence that same foundation brings to ordinary software.