Numerics

Need the whole and fractional parts of a float? The (x as i64) as f64 round-trip works — until the value doesn’t fit in an i64. trunc and fract split the float directly.

The cast detour everyone writes first:

1
2
3
4
5
6

let t = 7.75_f64;

let whole = (t as i64) as f64;
let frac = t - whole;
assert_eq!(whole, 7.0);
assert_eq!(frac, 0.75);

The direct version — no casts, no subtraction:

1
2
3

let t = 7.75_f64;
assert_eq!(t.trunc(), 7.0);
assert_eq!(t.fract(), 0.75);

Both parts keep the sign, and trunc() + fract() always reassembles the original — which also makes trunc the “toward zero” sibling in the rounding family from bite 303:

1
2
3
4

let t = -7.75_f64;
assert_eq!(t.trunc(), -7.0);
assert_eq!(t.fract(), -0.75);
assert_eq!(t.trunc() + t.fract(), t);

The real reason to drop the cast: as i64 silently saturates past ±2⁶³, so the round-trip hands back the wrong number without a peep. trunc has no range limit:

1
2
3
4
5

let big = 1e19_f64; // > i64::MAX

assert_eq!(big.trunc(), 1e19); // correct
assert_eq!((big as i64) as f64,
    9.223372036854776e18); // wrong, silently

The classic use case — splitting a quantity across two units:

1
2
3
4
5

let secs = 92.375_f64;
let mins = (secs / 60.0).trunc();
let rest = secs - mins * 60.0;
assert_eq!(mins, 1.0);
assert_eq!(rest, 32.375);

Note that fract on a negative number is negative — if you want “distance above the floor” instead (always in [0, 1)), that’s x - x.floor(), or x.rem_euclid(1.0).

round() sends every .5 away from zero — so summing rounded values drifts upward. round_ties_even rounds ties to the even neighbor and the bias cancels out.

Round a batch of exact midpoints with round() and watch the total inflate:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

let vals = [0.5_f64, 1.5, 2.5, 3.5];
// true sum: 8.0

let away: f64 = vals.iter()
    .map(|v| v.round())
    .sum();
assert_eq!(away, 10.0); // 1 + 2 + 3 + 4

let even: f64 = vals.iter()
    .map(|v| v.round_ties_even())
    .sum();
assert_eq!(even, 8.0); // 0 + 2 + 2 + 4

round() uses “half away from zero”: every tie moves in the same direction, so the errors all point the same way and accumulate. round_ties_even — banker’s rounding — sends ties to whichever neighbor is even, so roughly half go down and half go up:

1
2
3
4

assert_eq!(0.5_f64.round_ties_even(), 0.0);
assert_eq!(1.5_f64.round_ties_even(), 2.0);
assert_eq!(2.5_f64.round_ties_even(), 2.0);
assert_eq!((-2.5_f64).round_ties_even(), -2.0);

Only exact ties behave differently — 2.4999 and 2.5001 round the same either way:

1
2

assert_eq!(2.4999_f64.round_ties_even(), 2.0);
assert_eq!(2.5001_f64.round_ties_even(), 3.0);

This is the rounding mode IEEE 754 uses by default for a reason: it’s statistically unbiased over many operations. It’s also what Python’s built-in round does — if a value “rounds wrong” when porting Python code to Rust, this is why. Stable since Rust 1.77.

Every negative number has a real cube root — -27 is just -3. But (-27.0).powf(1.0 / 3.0) returns NaN. cbrt is the method that actually knows what a cube root is.

The naive formula falls over the moment the input goes negative:

1
2
3
4
5
6
7

let x = -27.0_f64;

let naive = x.powf(1.0 / 3.0);
assert!(naive.is_nan());

let root = x.cbrt();
assert_eq!(root, -3.0);

Why the NaN? powf computes x^y roughly as exp(y * ln(x)), and ln of a negative number doesn’t exist. powf can’t know that 0.3333333333333333 was meant to be ⅓ — it’s just a float that isn’t 1/3, so there’s no “odd root of a negative” special case to fall into. Mathematically fine operation, wrong tool.

The exponent being off also costs accuracy on positive inputs. 1.0 / 3.0 isn’t exactly one third, so powf computes a slightly different power:

1
2
3
4

let approx = 729.0_f64.powf(1.0 / 3.0);
assert_ne!(approx, 9.0); // 8.999999999999998

assert_eq!(729.0_f64.cbrt(), 9.0);

cbrt also keeps the edge cases sane where the formula produces nonsense: (-0.0).cbrt() is -0.0, and f64::NEG_INFINITY.cbrt() is NEG_INFINITY — not NaN.

Same story as this morning’s ln_1p (bite 301): when std ships a dedicated method for the operation you’re approximating with a formula, the method wins on both correctness and precision. cbrt has been stable since Rust 1.0 — and there’s a sqrt next to it you were already using.

(1.0 + x).ln() looks innocent, but for tiny x the 1.0 + step rounds away most of x’s digits before ln even sees them. ln_1p and exp_m1 do the + 1 internally, where nothing is lost.

Say you’re turning a tiny growth rate into a log-return:

1
2
3
4
5
6

let r = 1e-12_f64;

// 1.0 + r rounds away most of r
let naive = (1.0 + r).ln();
// ~1.0000889e-12 — 4 digits, then noise
assert!((naive - r).abs() > 1e-17);

The problem is spacing. Around 1.0, consecutive f64 values are about 2.2e-16 apart, so 1.0 + 1e-12 snaps to the nearest representable value — and the snap error (~9e-17 here) is enormous compared to r itself. ln then faithfully computes the log of the wrong number. You get 1.0000889e-12: four correct digits, then noise.

ln_1p computes ln(1 + x) without ever materializing 1 + x:

1
2
3

let precise = r.ln_1p();
// true value is r - r²/2 ≈ r - 5e-25
assert!((precise - r).abs() < 1e-20);

Full precision — the true answer differs from r only in the 25th decimal place, and ln_1p nails it.

The same trap exists in the other direction. e^x - 1 for tiny x:

1
2
3
4
5
6
7

// exp(r) ≈ 1.0000000000010000...
// subtracting 1.0 exposes the rounding
let naive = r.exp() - 1.0;
assert!((naive - r).abs() > 1e-17);

let precise = r.exp_m1(); // e^r - 1
assert!((precise - r).abs() < 1e-20);

And since they’re exact inverses, they round-trip cleanly:

1
2

let back = r.ln_1p().exp_m1();
assert!((back - r).abs() < 1e-26);

If x is comfortably large — say 0.1 and up — the naive forms are fine. But interest rates, probabilities, and per-step deltas live exactly in the tiny range where they aren’t. Both methods have been stable since Rust 1.0; they’re a rename away.

(dy / dx).atan() happily returns an angle — in the wrong quadrant half the time, and NaN at the origin. atan2 takes both components and gets every case right.

This morning’s bite 299 covered hypot, the robust way to get the length of a vector. Its companion problem is the angle — and the naive formula fails in a sneakier way.

atan alone can only return angles between −π/2 and π/2, because the division throws away the signs of the inputs:

1
2
3
4
5

let (dx, dy) = (-3.0_f64, 4.0);

// (-3, 4) is up-and-left: second quadrant
let naive = (dy / dx).atan();
assert!(naive < 0.0); // ...a negative angle?!

4.0 / -3.0 is the same number as -4.0 / 3.0, so atan can’t tell the second quadrant from the fourth. atan2 takes dy and dx separately, keeps both signs, and returns the full −π to π range:

1
2
3
4
5
6

let angle = dy.atan2(dx);
assert!((angle - 2.2142974355881810).abs() < 1e-12);

// naive answer is off by exactly π
let err = angle - naive;
assert!((err - std::f64::consts::PI).abs() < 1e-12);

It also survives the edge the division dies on — the origin:

1
2

assert!((0.0_f64 / 0.0).atan().is_nan());
assert_eq!(0.0_f64.atan2(0.0), 0.0);

Together with hypot, this is the complete cartesian-to-polar kit, and the round trip closes to within an ulp:

1
2
3
4

let r = dx.hypot(dy);
let theta = dy.atan2(dx);
assert!((r * theta.cos() - dx).abs() < 1e-12);
assert!((r * theta.sin() - dy).abs() < 1e-12);

Argument order trips everyone up once: it’s y.atan2(x) — the vertical component is self. If your angles keep coming out mirrored across the diagonal, you’ve swapped them.

(dx*dx + dy*dy).sqrt() looks harmless — until the squares overflow to infinity even though the answer would fit in an f64 just fine. hypot computes the same distance without ever squaring your inputs.

The textbook distance formula squares first, adds, then takes the root:

1
2
3
4

let dx = 3.0_f64;
let dy = 4.0_f64;
let dist = (dx * dx + dy * dy).sqrt();
assert_eq!(dist, 5.0);

Works great — until the components get big. f64::MAX is about 1.8e308, so squaring anything past ~1.3e154 blows straight through it:

1
2
3
4
5
6

let dx = 3.0e200_f64;
let dy = 4.0e200_f64;

// dx*dx is 9e400 — that's infinity in f64
let naive = (dx * dx + dy * dy).sqrt();
assert!(naive.is_infinite());

The final answer, 5.0e200, fits in an f64 with room to spare. It’s only the intermediate squares that overflow. hypot is built for exactly this — it computes sqrt(x² + y²) using a rescaling algorithm that never materializes the squares:

1
2
3

let dist = dx.hypot(dy);
assert!(dist.is_finite());
assert!((dist / 5.0e200 - 1.0).abs() < 1e-15);

(Note the relative comparison — hypot is accurate to within an ulp or so, not bit-exact.)

The same trick saves you at the other end of the scale: for tiny components the squares underflow to zero, and the naive formula returns 0.0 for a distance that plainly isn’t zero. hypot gets that case right too:

1
2
3

let tiny = 1.0e-200_f64;
assert_eq!((tiny * tiny).sqrt(), 0.0); // underflow
assert!(tiny.hypot(tiny) > 0.0);

If your coordinates are sane screen pixels, the naive formula is fine. But the moment the magnitudes are user-supplied or physics-scaled, reach for hypot — same one-liner, no cliff at 1e154.

a * b + c rounds twice and may compile to two separate operations. a.mul_add(b, c) computes the whole thing at full precision, rounds once, and on a CPU with FMA folds into a single instruction.

Two roundings vs. one

Floating-point math rounds after every operation. Writing a * b + c first rounds the product a * b to the nearest f64, then rounds again after adding c. That intermediate rounding throws away bits before they ever reach the sum:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

let a = 1.0_f64;
let b = 2.0_f64;
let c = 3.0_f64;

// Two roundings: (a*b) rounded, then (+ c) rounded
let separate = a * b + c;

// One rounding: a*b + c evaluated at full precision, rounded once
let fused = a.mul_add(b, c);

assert_eq!(separate, 5.0);
assert_eq!(fused, 5.0);

For tidy values they agree. The difference shows up when the product needs more bits than an f64 can hold and the addition would have recovered them — exactly the catastrophic-cancellation cases that wreck numeric code.

Where it pays off: Horner’s method

Polynomial evaluation is built out of multiply-then-add, so it’s the textbook home for mul_add. To evaluate 2x² + 3x + 4 with Horner’s method you nest the operations, and each step is one fused step:

1
2
3
4
5
6

fn poly(x: f64) -> f64 {
    // ((2*x) + 3)*x + 4
    2.0_f64.mul_add(x, 3.0).mul_add(x, 4.0)
}

assert_eq!(poly(2.0), 18.0); // 2*4 + 3*2 + 4

Each mul_add keeps full precision through the chain instead of rounding at every * and +. The same pattern carries dot products, which are a sum of products:

1
2
3
4
5
6
7

fn dot(a: &[f64], b: &[f64]) -> f64 {
    a.iter()
        .zip(b)
        .fold(0.0, |acc, (&x, &y)| x.mul_add(y, acc))
}

assert_eq!(dot(&[1.0, 2.0, 3.0], &[4.0, 5.0, 6.0]), 32.0);

The honest caveat

mul_add is faster only when the target CPU has a hardware fused-multiply-add instruction (modern x86-64 with FMA, AArch64, most others you’ll deploy to). Where it doesn’t, the standard library has to emulate the exact single-rounding semantics in software, and that emulation is slower than a plain a * b + c. So this is a hot-path tool: reach for it in tight numeric loops on hardware you control (or behind target-feature/target-cpu flags), and lean on it freely when you need the extra accuracy regardless of speed.

The bottom line

mul_add gives you a single rounding step and, on capable hardware, a single instruction for a * b + c. Use it in polynomial evaluation, dot products, and any multiply-accumulate loop where precision or throughput matters — and remember it can regress on FMA-less targets, so measure when speed is the goal.

Computing the average of two integers sounds trivial — until it overflows. The midpoint method gives you a correct result every time, no wider types required.

The classic binary search bug lurks in this innocent-looking line:

1

let mid = (low + high) / 2;

When low and high are both large, the addition wraps around and you get garbage. This has bitten production code in every language for decades.

The textbook workaround avoids the addition entirely:

1

let mid = low + (high - low) / 2;

This works — but only when low <= high, and it’s one more thing to get wrong under pressure.

Rust’s midpoint method handles all of this for you:

1
2
3
4
5
6
7
8
9

let a: u64 = u64::MAX - 1;
let b: u64 = u64::MAX;

// This would panic in debug or wrap in release:
// let avg = (a + b) / 2;

// Safe and correct:
let avg = a.midpoint(b);
assert_eq!(avg, u64::MAX - 1);

It works on signed integers too, rounding toward zero:

1
2
3

let x: i32 = -3;
let y: i32 = 4;
assert_eq!(x.midpoint(y), 0);  // rounds toward zero, not -∞

And on floats, where it’s computed without intermediate overflow:

1
2
3

let a: f64 = f64::MAX;
let b: f64 = f64::MAX;
assert_eq!(a.midpoint(b), f64::MAX);  // not infinity

Here’s a binary search that actually works for the full u64 range:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

fn binary_search(sorted: &[i32], target: i32) -> Option<usize> {
    let mut low: usize = 0;
    let mut high: usize = sorted.len();
    while low < high {
        let mid = low.midpoint(high);
        match sorted[mid].cmp(&target) {
            std::cmp::Ordering::Less => low = mid + 1,
            std::cmp::Ordering::Greater => high = mid,
            std::cmp::Ordering::Equal => return Some(mid),
        }
    }
    None
}

let data = vec![1, 3, 5, 7, 9, 11];
assert_eq!(binary_search(&data, 7), Some(3));
assert_eq!(binary_search(&data, 4), None);

Available on all integer types (u8 through u128, i8 through i128, usize, isize) and floats (f32, f64). No crate needed — it’s in the standard library.