Fast statistics and fussy arithmetic in Rust

On a RefGene genome-annotation workload, Fastmash 0.1.0 calculates four statistics in 26.6 milliseconds. GNU datamash 1.9 takes 112.8 milliseconds on the same Intel Core i7-8550U. That is 4.2 times as fast, with matching output. The published measurements include both the wins and the jobs where GNU datamash wins instead.

The arithmetic makes this interesting. Fastmash computes with an 80-bit floating-point format implemented in software, including software log and exp. That is quite a bit of homework for a program with “fast” in its name.

I built Fastmash in Rust to run familiar commands faster and make their numerical results reproducible across supported machines. To do both, it shares work between statistics and makes each arithmetic fast path earn its answer. Precision, evaluation order and rounding all stay part of the deal.

Four statistics in 26.6 milliseconds

The command behind that result asks for quartiles of exon counts in RefGene annotations, using field nine:

fastmash q1 9 median 9 q3 9 iqr 9 < refGene.txt

q1 and q3 are the first and third quartiles; iqr is their difference. On an AMD Ryzen 7 9800X3D, the same command is 3.8 times as fast. The benefit also shows up in decimal accumulation and grouped statistics. All times below are milliseconds, and lower is better.

01 / MEASURED

Same jobs. Two CPUs. Both sides of the story.

GNU datamash 1.9Fastmash 0.1.0

Intel Core i7-8550U

Elapsed time in ms · lower is better

RefGene quartiles

GNU112.8

Fastmash26.6

Sum + mean · 1M decimals

GNU324.7

Fastmash122.5

Exon statistics by gene

GNU159.4

Fastmash83.1

Large sort · disk spill

GNU572.4

Fastmash694.1

0 ms750 ms

AMD Ryzen 7 9800X3D

Elapsed time in ms · lower is better

RefGene quartiles

GNU50.3

Fastmash13.1

Sum + mean · 1M decimals

GNU105.9

Fastmash46.8

Exon statistics by gene

GNU62.2

Fastmash36.5

Large sort · disk spill

GNU180.0

Fastmash309.7

0 ms350 ms

The setup matters. These measurements use the 0.1.0 release binary on two native Linux machines. Both programs got the same input and limits, and their outputs were checked before a timing counted. Each value averages the median from two sessions, with six measured runs per session and alternating program order.

Input came from regular files with a warm page cache. Each complete command had a 512 MiB memory limit, no swap and disk-backed temporary storage. The method and full results give the hardware, commands and individual comparisons.

Summing and averaging a million decimals is 2.7 times as fast on Intel and 2.3 times as fast on AMD. The disk-spill row is a useful reminder that GNU datamash still has a few tricks up its sleeve. More on that below.

Keep the last digit under control

GNU datamash 1.9 on x86-64 Linux uses 80-bit extended arithmetic and the C math library. Different builds can disagree in the last printed digit. That gives a byte-for-byte comparison something extra to complain about.

Fastmash's numerical contract fixes one answer for a given version, input and explicit settings on every supported machine. Getting there starts with the number format.

Rust's ordinary f64 has 53 bits of significand precision. The chosen 80-bit format has 64. Switching to f64 would change the computation. Add these values in order with binary64 arithmetic and the middle 1 disappears:

9007199254740992
1
-9007199254740992

Fastmash retains it:

printf '9007199254740992\n1\n-9007199254740992\n' | env LC_ALL=C fastmash sum 1
# 1

Precision is only half the problem. Floating-point additions round after each step, so rearranging a sum can change its answer. A parallel reduction or a different grouping strategy must preserve the promised evaluation order if it is to preserve the result. “Same numbers” does not automatically mean “same calculation.”

Fastmash represents the 80-bit values explicitly and computes their arithmetic in software. The implementation combines specialized integer arithmetic with a fallback using rustc_apfloat. Transcendental functions need their own treatment: Rust's f64::ln documentation, for example, does not promise identical precision across platforms. Calling a host math function would leave part of the numerical contract outside Fastmash's control.

Parse once and sort once

That RefGene command asks four questions about the same column. All four can share the expensive preparation: convert the numbers, collect the samples and sort them.

Fastmash builds a plan before consuming the records. Compatible operations share the work they need: sum and mean can use one accumulated sum, while q1, median, q3 and iqr can use one collection of samples and one sort. Asking for more statistics does not have to mean doing all the preparation again.

The conversion reuse fits in a few lines. The operation collector either takes an earlier operation's parsed value or converts the field itself:

let parsed = match plan.conversion_source {
    Some(source) => previous[source].parsed,
    None => record.number(field, line, options)?,
};
state.parsed = parsed;
let Some(value) = parsed else {
    continue;
};

The easy detail to miss is None. When a missing value is skipped, that absence must also replace the previous record's parsed value. Otherwise the next operation could help itself to a number left over from the previous record.

Sharing sorted samples has a similar correctness detail. The sample collection tracks whether it has already been sorted. Repeated sorts can change the traversal of NaNs, so sorting once also preserves the intended behavior for those inputs.

Sharing work saves conversion, accumulation and sorting. The end-to-end timings measure the combined effect with the rest of the program; they do not tell us how much of the speedup to assign to each choice.

Make the arithmetic fast path earn its answer

Rust's u128 gives the multiplication fast path a useful starting point: the exact product of two 64-bit significands fits inside it. For normal finite operands, Fastmash can keep the discarded bits until it makes one rounding decision.

The relevant part of the multiplication fast path is:

let mut significand = product >> shift;
let remainder = product & ((1u128 << shift) - 1);
let half = 1u128 << (shift - 1);
if remainder > half || (remainder == half && significand & 1 != 0) {
    significand += 1;
}

shift separates the retained significand from the discarded low bits. A remainder above halfway rounds up. At exactly halfway, the last retained bit chooses the even result. Surrounding code checks the operand and result ranges, handles carry, and sends cases outside this path to the fallback. The shortcut still follows the specified precision and rounding rule.

Then there are log and exp. Their exact mathematical answers generally cannot be represented with a finite number of bits, so Fastmash aims for the correctly rounded 80-bit result: the representable value nearest to the mathematical answer, with ties resolved consistently.

Its fast exponential reduces the input to a smaller range, combines a table lookup with a polynomial, and bounds the approximation error. The final decision is just four lines:

let m = multiply(entry as i128, polynomial);
let error = n.abs() + 64;
let lower = scaled_binary80(m - error, n)?;
(scaled_binary80(m + error, n)? == lower).then_some(lower)

Here m is the computed approximation in fixed-point units, and n supplies the power-of-two scaling. scaled_binary80 rounds an endpoint to the target format. If both ends of the error interval round to the same value, every answer inside the interval rounds there too. The uncertainty is too small to affect the answer. That agreement is the certificate allowing the fast path to return Some(value).

02 / CONCEPTUAL

One interval. Two possible outcomes.

A The shortcut is certified

Both endpoints round to value A.
Return Some(value).

B The fallback has work to do

The endpoints round to different values.
Return None; use the fallback.

If the interval crosses a rounding boundary, or the input or result is outside this path's supported range, it returns None. The caller then uses its arbitrary-precision path. Near the format's range boundaries, it also checks wider results rounded downward and upward before accepting a conversion. Inputs that cannot be settled within the numerical engine's limits are refused with exit status 77.

Rust's Option makes the handoff explicit: either the fast path has a certified answer, or the caller has more work to do. The error bound still needs mathematical justification and testing. The compiler is helpful, but it has not volunteered to prove the polynomial.

Check the answers independently

Fastmash's testing layers check command behavior and numerical accuracy separately, then require correct output in the performance measurements.

The command corpus contains more than 3,000 cases with exact expected standard output, standard error and exit status. For GNU-compatible behavior, those expectations were observed from GNU datamash 1.9, twice per case. Intentional differences have their own documented expectations. The corpus runs with piped input and repeats applicable cases with regular-file input, because those transports exercise different execution paths.

Numerical tests use independent high-precision reference values. Logarithm and exponential fixtures come from MPFR, including difficult inputs close to rounding boundaries and the limits of the format. Agreement between two paths in the same program is useful, but both could share a mistake. Bugs are perfectly capable of teamwork.

Reproducibility and compatibility are related but distinct. Fastmash can differ from a local GNU build in the last printed digit of some transcendental results, and it has explicit NaN rules. The differences guide records those choices. Numerical output changes are versioned and listed in the changelog.

Where GNU datamash still wins

GNU datamash wins the large sort in the table on both hosts: Fastmash took 694.1 ms against 572.4 ms on Intel, and 309.7 ms against 180.0 ms on AMD. The full results also retain a slower many-key job and a geometric mean over 200,000 distinct values that cost 7.4 ms extra on Intel. Some other comparisons were inconclusive across the two sessions. None of those establish a win.

The released 0.1.0 package targets Linux x86-64 with glibc, including Linux under WSL2. Portable numerical semantics provide a consistent target for supported machines; they do not by themselves provide operating-system support or prove performance on unmeasured hardware.

Try it on your own machine

With Fastmash installed, a small version of the quartile command is:

printf '1\n2\n3\n4\n' | env LC_ALL=C fastmash q1 1 median 1 q3 1 iqr 1
# 1.75    2.5    3.25    1.5

The output fields are separated by tabs. The benchmark kit runs eight representative jobs, downloads the public datasets, generates synthetic data and checks output agreement before reporting times. It uses a simpler timing setup than the controlled measurements above, so its numbers need not match the table.

If you are building a numerical tool in Rust, pin down the precision, rounding and evaluation order first. Then look for work that several operations can share, and common cases whose answers can be certified cheaply. There is plenty of room to get clever about the implementation. The last digit should be boring.