Exact reductions

include/eon/ExactReduction.h provides two kinds of reduction.

• eonc::exact::dot, sum, sumStrided, squaredNorm, norm and normalized return the exact value rounded once to the nearest double, ties to even. The optimizers (L-BFGS, CG, FIRE, Quickmin, steepest descent) and the norm and rms convergence tests use them for every scalar they branch on. sumStrided is sum for a column of a row-major atom matrix.

• eonc::exact::orderedDot and orderedNorm add rounded products in a tree fixed in the source: four lanes, then (l0 + l1) + (l2 + l3), compiled without contraction. They are not exact, but their bits do not depend on the compiler, the Eigen release or the CPU. The L-BFGS two-loop recursion (which sets the step vector) and the per-atom lengths of the step caps and push-apart use them.

Why

Floating-point addition does not associate, so a sum depends on the order of its terms. Eigen chooses that order, and the order changes with the release and the SIMD width: for dot and squaredNorm, Eigen 3.4 adds two packet accumulators and Eigen 5.0 four (Eigen/src/Core/InnerProduct.h). With 39 coordinates (LJ13) under SSE2 the two releases build different summation trees.

Optimizers amplify a one-ulp difference. The basin-hopping LJ13 case in test_jobs first diverged at the finite-difference curvature of the first L-BFGS step (C in LBFGS::getStep differed by one ulp). The run then ended at 1640 force calls with Eigen 3.4 and 1700 with Eigen 5.0, under gcc 14 and gcc 15 alike. A plain dot product of length n carries an error of up to n u sum |x_i y_i| (u = 2^-53), so two orders can differ by twice that; near a convergence threshold or an L-BFGS curvature test, that difference flips a decision.

The mean force removed when removeNetForce is set is that same exact sum of each column, divided by the atom count. Eigen’s colwise().sum() follows the packet width, so the mean differed between x86 and NEON and the LJ13 basin hop took 1671 force calls on Apple clang instead of 1640. An exact mean left 1658. Rewriting the L-BFGS curvature probe and the two-loop step as a multiply followed by an add left the count at 1658, so that step was not the remaining gap.

The gap is the LJ pair evaluation. On x86 without FMA in the compile ISA, rgpot compiles the pair loop and the inlined kernel in accumulateHw: contraction off, and an explicit fma only for the force update. ARM has FMA in the base ISA, so that clone is not built. Apple clang’s default contraction then fuses dx*dx + dy*dy + dz*dz and the 12-6 products (b*(a-1) - shift, 2*a - 1). The LJ objects are compiled with contraction off, the same flag the minimum-image objects already use, and the force update stays __builtin_fma. That is the rounding of the Linux clone.

How

Every finite double is an integer significand m < 2^53 times 2^e. x_i y_i is therefore the 106-bit integer m_x m_y times 2^(e_x + e_y), with no rounding. The products go into a fixed-point integer accumulator that spans the whole exponent range, 135 words of 32 bits in signed 64-bit storage. Integer addition associates, so the accumulator contains the same exact sum for any order. One final rounding gives RN(sum x_i y_i).

This is the integer form of the error-free splitting of the Ozaki scheme (Ozaki, Ogita, Oishi and Rump 2012; Ozaki, Uchino and Imamura 2025). There, each operand splits into slices whose products are exact and whose sums are exact in any order. Compensated algorithms such as Dot2 (Ogita, Rump and Oishi 2005) improve accuracy but stay order-dependent. Reproducible summation in the sense of Demmel and Nguyen 2015 is order-independent without being exact. A correctly rounded result is both.

dot first tries a faster route to the same number. It runs a compensated dot product (Dot2) in four lanes, with TwoProduct by FMA where the CPU has it and Veltkamp-Dekker splitting where it does not; both give the same exact error term. It also bounds the error of the compensation sum a posteriori. When the exact value then lies strictly between the two midpoints around r = RN(s + c), r is the correctly rounded result. Near a tie, or when a range guard fails, the integer accumulator decides. Both routes return RN(sum x_i y_i), so the route never shows in the bits.

Non-finite inputs fall back to the plain IEEE sum, so NaN and infinity propagate unchanged. An exact zero returns +0.

Cost

Measured on one Xeon Platinum 8358 node (Elja, gcc 15.3, Eigen 5.0.1), with run directories on local disk, one warm-up run, then five interleaved runs; the table gives medians.

workload	before (s)	after (s)	change
LJ13 basin hopping, L-BFGS, 200 steps	0.221	0.241	+9 %
basin-hopping example, CG	2.58	2.18	within noise
Al NEB, 7 images, L-BFGS	0.957	0.998	+4 %
Pt saddle search, dimer	0.267	0.269	+1 %

LJ13 has the cheapest force call of the four, so it shows the largest relative cost. The trajectories differ from the old ones (14469 against 14348 force calls for LJ13), so part of each change is a different path, not overhead. In isolation, dot costs 3.4 ns per term at n = 39 and 1.0 ns per term at n = 1008, against 0.6 and 0.2 ns for Eigen; orderedDot matches Eigen at n = 39 and costs 1.9 times as much at n = 1008.

Checks

• client/validation/exact_reduction_cases.py writes reference cases from exact rational arithmetic: cancellation, ties, sticky bits, subnormal results, the overflow threshold, and optimizer-shaped vectors. test_exact_reduction checks each case bitwise and under 16 permutations.

• client/validation/exact_reduction_derivations.py (sympy) derives the Ozaki slice condition beta >= (p + log2 n) / 2, the portable 128-bit product and the accumulator piece split. It also lists the Eigen 3.4 and 5.0 trees for n = 39.

• client/validation/exact_reduction_bounds.sollya bounds the plain-sum error, the norm error (1.5 u) and the accumulator range.

• formal/ (Lean 4 and Mathlib, lake build) checks four contracts. First, floating-point sums on a common grid are exact in every parenthesisation within the exponent range, so any two orders agree. Second, the Ozaki slice-product condition. Third, the final rounding step is round to nearest, ties to even. Fourth, accumulator words cannot overflow between carry normalisations.