Environments, Dependencies, and Precision

The same source and input files need not produce the same bits—or even the same trustworthy conclusion—when the compiler, numerical libraries, precision, hardware, or reduction order changes. A reproducible environment makes those dependencies visible. A credible experiment then does something more demanding: it deliberately varies the environmental choices most likely to expose numerical fragility while holding the physical problem fixed. Pinning and perturbing the environment are complementary tests, not rival philosophies.

Required background. Calculus and Linear Algebra for Relativity supplies conditioning, norms, and floating-point scale estimates; Differential Equations for Relativity supplies truncation error, stiffness, and tolerance control.

Helpful background. Formulations of General Relativity connects continuum variables to their implementations; Probability and Inference for Gravity helps distinguish numerical variation from ensemble uncertainty.

What belongs to the environment

An environment record should be sufficient for a competent reader to reconstruct the numerical semantics of a run. A package-name list is not enough. At minimum, record the following:

Layer	Details to preserve	Why it can matter
source	revision, local changes, generated code, build recipe	selects the actual equations and algorithms
toolchain	compiler or interpreter, version, flags, linker	changes optimization, contraction, and math-library calls
dependencies	exact versions, build variants, solver backends	changes tolerances, ordering, and implementation details
arithmetic	format, rounding mode, fused operations, mixed precision	changes accumulated roundoff and exceptional values
execution	operating system or container digest, CPU or accelerator model	changes instruction sets and library kernels
parallelism	process and thread counts, decomposition, scheduler settings	changes reduction order and sometimes the mesh itself
run inputs	canonical configuration, data versions, checksums, random seeds	fixes the finite problem actually solved

Some entries must be pinned for a clean rerun; others must be varied for a robustness check. Pin a dependency graph so that an old result can be reconstructed. Vary precision, optimization, thread count, and a genuinely independent solver when they could change the conclusion. A container can freeze user-space software, but it does not freeze the processor, accelerator driver, parallel schedule, external data service, or undocumented host assumptions.

The IEEE floating-point standard specifies formats and operations, including rounding and exceptional values, but it does not promise that two algebraically equivalent programs execute the same operation sequence; see IEEE 754-2019, clauses 3–10. Record compiler transformations that relax those semantics. Flags commonly described as “fast math” may reassociate expressions, ignore signed zeros, contract operations, or assume that NaNs and infinities do not occur. Such a build may be useful, but its assumptions are scientific inputs whenever they affect a stopping condition or reported observable.

Precision, conditioning, and error budgets

For a basic rounded operation, the standard model is

$$\operatorname{fl}(a\circ b)=(a\circ b)(1+\delta), \qquad |\delta|\le u,$$

when the exact result is in range and the operation is one to which the model applies. Here $u=2^{-24}\simeq5.96\times10^{-8}$ for binary32 and $u=2^{-53}\simeq1.11\times10^{-16}$ for binary64. This is a local bound, not a promise that a long calculation has relative error $u$. Cancellation, ill-conditioning, a long reduction, unstable recurrences, underflow, overflow, and branch decisions can magnify the effect. Goldberg develops the formats and the main failure mechanisms in Goldberg 1991, §§ 1–4, while backward error and conditioning are treated systematically in Higham 2002, Chapters 1–4.

Three errors should be separated:

$$\begin{aligned} \|Q_{h,\tau,u}-Q\|\le{}& \underbrace{\|Q_h-Q\|}_{\text{discretization}}\\ &+\underbrace{\|Q_{h,\tau}-Q_h\|}_{\text{solver and time tolerances}}\\ &+\underbrace{\|Q_{h,\tau,u}-Q_{h,\tau}\|}_{\text{roundoff and implementation}}. \end{aligned}$$

The split is conceptual because the terms can interact. Raising precision cannot repair an inconsistent boundary condition or an under-resolved mesh. Refining the mesh can make roundoff worse when a derivative divides nearly cancelling numbers by a smaller spacing. Tightening a nonlinear tolerance below the attainable arithmetic floor wastes work and may cause stagnation rather than accuracy.

A precision study therefore refines the discretization and arithmetic independently. Hold $u$ fixed while varying $h$ to locate the truncation regime; hold $h$ fixed while moving from binary32 to binary64 and, when necessary, a multiprecision reference. If the answer changes at binary64 but stabilizes at higher precision, diagnose conditioning before merely adopting the expensive format. If it does not stabilize, the dominant failure may be mathematical or algorithmic rather than arithmetic.

Scalar-wave environment benchmark

Use the periodic scalar wave from the preceding experiment design: $u(t,x)=\sin(2(x-t))$ on $0\le x<2\pi$, the fourth-order centered derivative, fourth-order Runge–Kutta, maximum Courant factor $\lambda_{\max}=1/4$, and $N=32,64,128,256$. Reach $T=1$ with $m=\lceil T/(\lambda_{\max}\Delta x)\rceil$ equal steps of size $\Delta t=T/m$, so the actual Courant factor is at most $1/4$. The continuum energy is $\mathscr H=4\pi$, and the exact semidiscrete solution—including the stationary contribution from the sampled continuum initial data—supplies an analytic reference independent of the evolution implementation.

The finite benchmark has four environment arms:

1. a scalar binary64 build with contraction disabled;
2. an otherwise identical optimized binary64 build with documented contraction and reassociation settings;
3. binary32 storage and arithmetic; and
4. a reference evaluation of the modified wavenumber and full semidiscrete field formula with at least 80 decimal digits.

Run each compiled arm at one and four threads, record the compiler and flags, and repeat it in a clean environment from the same checksum-bearing input. The evolution itself may be deterministic and serial for this small problem; the thread variation still tests whether library initialization or analysis reductions introduce a hidden dependency. Report the root-mean-square field error, energy drift, constraint $\Phi-D_4u$, and SHA-256 digest of the canonical input beside every result.

The predeclared numerical criteria are:

• binary64 must recover fourth-order convergence over the two finest pairs and give $E_{256}<10^{-6}$;
• its fields must approach the 80-digit semidiscrete prediction under time-step refinement; use the full analytic Runge–Kutta amplification factor when checking a finite time step to $10^{-12}$, since a single halving only estimates, rather than removes, temporal error;
• scalar and threaded analysis paths must agree within $10^{-13}$ in $E_N$ and relative energy drift, even if their files are not bitwise identical;
• binary32 is expected to leave the clean fourth-order regime near the finest error, because the predicted $E_{256}\simeq2.45\times10^{-7}$ is only a few units of its machine precision; that breakdown is a diagnostic, not a failed binary64 result.

These are acceptance thresholds for an implementation to measure, not reported measurements from this prose page. The independent check is analytic: compute

$$\widetilde{k}=\frac{8\sin(k\Delta x)-\sin(2k\Delta x)}{6\Delta x}$$

at high precision and insert it into the full semidiscrete field formula on the experiment-design page. A second run of the same evolution routine is useful for reproducibility, but it is not independent verification.

Reduction order and nondeterminism

Floating-point addition is not associative. Parallel sums can therefore depend on how a runtime partitions work:

$$\operatorname{fl}\!\bigl(\operatorname{fl}(a+b)+c\bigr) \ne \operatorname{fl}\!\bigl(a+\operatorname{fl}(b+c)\bigr).$$

This matters for norms, conserved quantities, elliptic inner products, and global time-step estimates. A change of the last few bits is not automatically a scientific failure. It becomes consequential when it is amplified by an ill-conditioned solve, changes an adaptive branch, or exceeds the error budget of the observable.

Choose the reproducibility level before running:

• bitwise: appropriate for regression tests on a fixed supported platform when the operation order is controlled;
• numerical tolerance: appropriate when different legal evaluation orders should approximate the same deterministic quantity;
• statistical: appropriate for explicitly random ensembles, with a seed policy and distributional comparison.

If bitwise sums are required across decompositions, use a documented reproducible reduction rather than hoping the scheduler repeats itself. Algorithms for exact or reproducible accumulation make the trade between speed, storage, and determinism explicit; see Demmel and Nguyen 2013, §§ 2–5. Regardless of method, preserve process counts and decomposition because they can alter more than reduction order—for example, the placement of mesh-refinement boundaries.

Rerun sequence and failure criteria

A complete environment check can be reconstructed through the following sequence:

1. Resolve a lockfile or immutable image and verify the source, input, and dependency checksums before building.
2. Emit compiler or interpreter versions, flags, linked numerical libraries, arithmetic format, rounding mode, host and accelerator identifiers, thread and process counts, and source revision into a machine-readable run record.
3. Execute the scalar-wave resolution sweep for each precision and build arm; refuse to combine products carrying different input digests.
4. Recompute norms with the independent high-precision formula, halve the time step at fixed resolution, and compare one-thread with four-thread results at the declared tolerance.
5. Recreate the pinned environment in a clean directory or host, repeat the sweep without manual intervention, and preserve both the numerical comparison and any bitwise differences.

Reject the claim if the environment cannot be reconstructed, a dependency or input is unversioned, the binary64 refinement sequence does not reach its expected regime, the result changes beyond tolerance with threads or supported platforms, the high-precision check exposes conditioning not included in the uncertainty, or exceptional values are silently suppressed. Also reject a precision comparison that changes the algorithm at the same time: it cannot isolate arithmetic effects.

Check your understanding

In binary64, evaluate $f(x)=\sqrt{1+x}-1$ at $x=10^{-16}$. A direct implementation returns zero, while a higher-precision calculation gives approximately $5\times10^{-17}$. Is binary64 intrinsically too imprecise for the answer?

Solution

No. The direct expression subtracts two nearly equal rounded numbers. Rationalizing it gives

$$f(x)=\frac{x}{\sqrt{1+x}+1},$$

which evaluates to approximately $5\times10^{-17}$ in binary64 without the destructive subtraction. The function itself is well-conditioned here: $x f'(x)/f(x)\to1$ as $x\to0$. The direct algorithm loses accuracy through cancellation; the rationalized algorithm is a more stable evaluation of the same function. Higher precision is a valuable independent check. The environment record must still preserve the math-library and arithmetic assumptions used by either form.

References

• Demmel, James, and Hong Diep Nguyen. “Fast Reproducible Floating-Point Summation.” In 2013 IEEE 21st Symposium on Computer Arithmetic, 163–172. IEEE, 2013. DOI.
• Goldberg, David. “What Every Computer Scientist Should Know About Floating-Point Arithmetic.” ACM Computing Surveys 23, no. 1 (1991): 5–48. DOI.
• Higham, Nicholas J. Accuracy and Stability of Numerical Algorithms. 2nd ed. Philadelphia: Society for Industrial and Applied Mathematics, 2002. DOI.
• IEEE. IEEE Standard for Floating-Point Arithmetic. IEEE Std 754-2019. New York: IEEE, 2019. DOI.