Numerics

(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.