Faster floating point math with Rust’s new API

Floating point math is often slower than integer math because the compiler is being conservative about how it optimizes your code. While some programming languages already had solutions of a sort, until now Rust did not have a good stable way to deal with this limitation. But now, starting in version 1.98, Rust will allow telling the compiler it can optimize your code further—but with extra control so that you can still write numeric algorithms with minimal rounding errors.

In this article you will learn:

Why by default the compiler won’t optimize floating point math as much as it does integer math.
Rust’s new API to solve this limitation.
Examples of using this new API, its speed impact, and how you can control where it is used.

Summing integers is fast

I’m going to start with an example using integers, as a baseline of what sort of performance is possible.

To get the fastest code generation, I’m telling Rust that it’s not 2004
anymore,
and that it can generate CPU instructions that require modern hardware,
namely x86-64 machines from the past 10 years or so. Specifically, all
the code in this article is being compiled with
RUSTFLAGS="-C target-cpu=x86-64-v3". (For maximum compatibility, in
real-world usage you could provide a fallback implementation for older
computers.)

Here’s a Rust function to sum a slice of int64 numbers:

fn naive_sum_i64(values: &[i64]) -> i64 { let mut total = 0; for value in values { total += value; } total }
I’ll omit the code to expose this to Python, but it’s a variant of the Rust/Python code in a previous article.

To benchmark it, I’ll create an array of integers in NumPy:

import numpy as np DATA_INT = np.ones((1_000_000,), dtype=np.int64) assert naive_sum_i64(DATA_INT) == 1_000_000
And now I can measure the speed of summing this array:

| Code | ➘ Elapsed µ-seconds | ➘ CPU instructions per value | |
|---|---|---|---|
| naive_sum_i64(DATA_INT) | 168.1 | 0.5 |

➘ Lower numbers are better

That’s 0.5 CPU instructions per value! How does that even work?

Probably the compiler is using specialized Single Instruction, Multiple Data (SIMD) CPU instructions, that do batch operations on multiple values at once. The i7-12700K CPU I’m using here has 256-bit SIMD instructions, meaning it can do some specific operations on four 64-bit integers at a time. If there’s a specialized SIMD summing CPU instruction, the CPU would only need to loop 250,000 times and then sum 4 integers in each iteration.

And in fact:

| Code | ➘ Elapsed µ-seconds | ➘ CPU instructions | 256-bit SIMD integer instructions | |
|---|---|---|---|---|
| naive_sum_i64(DATA_INT) | 156.8 | 521,280 | 250,003 |

➘ Lower numbers are better

In short, by using a specialized SIMD instruction, my CPU can sum integers very quickly.

Summing floats is slow?!

But what about floats—are they fast too?

Again, I’ll create a million floating point values:

```

Array of 1M float64 values between 0 and 1.

DATA = np.random.random((1_000_000,))
```
I’ll implement a simple floating point sum function:

fn naive_sum(values: &[f64]) -> f64 { let mut total = 0.0; for value in values { total += value; } total }
And compare the performance of summing integers and floats:

| Code | ➘ Elapsed µ-seconds | ➘ CPU instructions | 256-bit SIMD integer instructions | 256-bit SIMD float instructions |
|---|---|---|---|---|
| naive_sum_i64(DATA_INT) | 151.9 | 521,214 | 250,003 | 0 |
| naive_sum(DATA) | 595.2 | 1,458,269 | 0 | 0 |

➘ Lower numbers are better

The floating point sum is much slower than the integer sum, and the compiler didn’t use SIMD float operations. Why the difference?

Floating point operations aren’t commutative

Like most compilers, when Rust compiles your code in release mode it
optimizes your code, transforming it in a variety of ways to (hopefully)
make it faster. But there’s a promise compilers make when they do this:
the optimized code will behave exactly the same as the unoptimized
code.

If I add two integers a and b, a + b == b + a. That gives the
compiler plenty of scope to optimize how the code runs, for example by
using SIMD operations that might slightly change the order of additions.

Floating point numbers are different. For example, because floating point numbers span such a range of values, from tiny to huge, adding a sufficiently large number to a sufficiently small number results in that same large number:

print( "Does adding a small number do nothing?", 1e16 + 1.0 == 1e16 )
Does adding a small number do nothing? True
More broadly, for floating point numbers, a + b is not always the same
as b + a, at least once you’re adding multiple numbers in a row. Let’s
say I have an array that starts with 1e16 followed by many 1.0
values, and another that is the reverse. Summing these arrays will give
different results:

import math HIGH_VALUE_FIRST = np.ones((1_000_000,), dtype=np.float64) HIGH_VALUE_FIRST[0] = 1e16 HIGH_VALUE_LAST = np.ones((1_000_000,), dtype=np.float64) HIGH_VALUE_LAST[-1] = 1e16 print( "Is the sum the same?", naive_sum(HIGH_VALUE_FIRST) == naive_sum(HIGH_VALUE_LAST) )
Is the sum the same? False
Because the order of summing impacts the result, the compiler assumes,
as it should, that I asked for this particular order for a reason. As a
result, the compiler will not reorder these operations. Nor will it
apply any other optimization that might change the results, even if the
resulting code is slower.

Rust’s new algebraic operators: telling the compiler when to be flexible

While conservatism on the part of the compiler is the right default,
sometimes you as the programmer know that re-ordering operations isn’t a
problem. In that situation, it would be good to be able to tell the
compiler that while over here it should not change the order of
operations, over there it’s actually fine.

Starting in Rust 1.98 there is a new feature that allows just that. In addition to the normal arithmetic operations you can do on floating point numbers, there are new set of so-called “algebraic” arithmetic operators that per the documentation “allow the compiler to optimize floating point operations using all the usual algebraic properties of real numbers”, including changing the order of operations.

Rust 1.98’s release date is August 20, 2026. Since I wrote the article before then, I ran the code in this article with the

"beta"channel, which includes the same functionality.

An example: Optimized pairwise summation

Let’s see these operators in action, and how they allow for fast code.

Because summing floating points can have unexpected behavior—remember
how 1e16 + 1.0 == 1e16—if you want to sum an a large number of
floating point numbers, there are various algorithms you can use to
minimize the resulting rounding errors. The trade-off between them is
usually about how fast they are vs how accurate.
numpy.sum()
mostly uses an algorithm called pairwise
summation that has a
good balance between speed and reducing accumulated errors. Certainly it
has less error accumulation than just adding one float at a time in
order, as I did with naive_sum() above.

The basic idea in pairwise summation is to split the array into two, sum each side recursively with the same algorithm, and then sum the resulting two floating point numbers. Below a certain threshold of array size—128 in NumPy’s case—summation is done normally… and that can be done in any order.

Let’s implement this algorithm in Rust. When adding those top two floating point numbers, I’m going to use normal addition, where the compiler can’t reorder operations or otherwise do anything that would change the output. Once the function hits the threshold of doing normal summing, I’ll switch to algebraic adds because at this point I don’t care what the order of addition is, I just want speed.

fn pairwise_sum(values: &[f64]) -> f64 { let n = values.len(); if n > 128 { // Precise addition of two recursive applications // of the algorithm: let half = n / 2; pairwise_sum(&values[0..half]) + pairwise_sum(&values[half..n]) } else { // 😎 Normal addition, where the order doesn't matter // as far as the algorithm is concerned. Thus using a // algebraic add is fine—and that gives the compiler // permission to optimize aggressively. let mut total: f64 = 0.0; for value in values { total = total.algebraic_add(*value); } total } }
This implementation of pairwise implementation is fast, in fact it’s faster than NumPy’s implementation:

| Code | ➘ Elapsed µ-seconds | ➘ CPU instructions | 256-bit SIMD float instructions | |
|---|---|---|---|---|
| naive_sum(DATA) | 563.1 | 1,458,279 | 0 | |
| np.sum(DATA) | 190.7 | 2,191,767 | 0 | |
| pairwise_sum(DATA) | 144.5 | 🏆 | 1,298,028 | 270,336 |

➘ Lower numbers are better

And it has the same (general) precision as NumPy’s implementation, as it implements the same algorithm:

```
from math import fsum

fsum() uses an accurate summation algorithm that gets rid of any

avoidable floating-point error.

assert fsum(HIGH_VALUE_FIRST) == fsum(HIGH_VALUE_LAST)
correct_sum = fsum(HIGH_VALUE_FIRST)
print(
"naive_sum() error: ",
naive_sum(HIGH_VALUE_FIRST) - correct_sum
)
print(
"np.sum() error: ",
np.sum(HIGH_VALUE_FIRST) - correct_sum
)
print(
"pairwise_sum() error:",
pairwise_sum(HIGH_VALUE_FIRST) - correct_sum
)

naive_sum() error: -1000000.0
np.sum() error: -14.0
pairwise_sum() error: -6.0
```
The difference in error between

np.sum()andpairwise_sum()is not meaningful, it’s just luck; the key point is they have a similarly low-order-of-magnitude error compared tonaive_sum().

As an added bonus, the code structure pairwise_sum() uses to encourage
generation of SIMD code is much more succinct than
NumPy’s.

Another example: Sum of squared differences

You can do more than just add numbers with algebraic operators. In the following examples I will calculate the sum of squared differences between two arrays.

I’m going to need two arrays:

DATA1 = np.random.random((1_000_000,)) DATA2 = np.random.random((1_000_000,))
Here’s an implementation with normal arithmetic operators:

fn ssd_normal(arr1: &[f64], arr2: &[f64]) -> f64 { assert_eq!(arr1.len(), arr2.len()); let mut total = 0.0; for (val1, val2) in arr1.iter().zip(arr2) { total += (val1 - val2).powi(2); } total }
And here’s what it looks like with algebraic operators; it’s not clear
to me if f64.powi() will use algebraic operators or not, so I did so
explicitly:

fn ssd_optimized(arr1: &[f64], arr2: &[f64]) -> f64 { assert_eq!(arr1.len(), arr2.len()); let mut total: f64 = 0.0; for (val1, val2) in arr1.iter().zip(arr2) { // 😎 Algebraic operations to allow more compiler // optimizations: let diff = val1.algebraic_sub(*val2); let squared_diff = diff.algebraic_mul(diff); total = total.algebraic_add(squared_diff); } total }
First, a quick test to make sure the results are similar:

print(ssd_normal(DATA1, DATA2)) print(ssd_optimized(DATA1, DATA2))
166770.0055951995 166770.00559520238
Next, I’ll measure performance:

| Code | ➘ Elapsed µ-seconds | ➘ CPU instructions per value | |
|---|---|---|---|
| ssd_normal(DATA1, DATA2) | 628.7 | 4.5 | |
| ssd_optimized(DATA1, DATA2) | 371.1 | 🏆 | 1.0 |

➘ Lower numbers are better

Using algebraic operations enables the compiler to generate code that runs twice as fast on my computer.

Go speed up some numeric code!

The pairwise summation algorithm is a great example of why you want both kinds of operators, strict and lenient:

If you only use strict in-order operations, the code is slower.
If you only use lenient algebraic optimize-as-you-like operations, the compiler might optimize the algorithm out of existence, losing the accuracy the algorithm aims for.

The implementation I showed above benefits from using both normal addition—for accuracy—and algebraic addition—for speed—in different parts of the algorithm.

If you are writing numeric code with Rust, your code may benefit too—give it a try once you can use Rust 1.98. And if you’re not using Rust yet, this is another good reason to switch.