Tracking Issue for num_midpoint

[Urgau]

Feature gate: #![feature(num_midpoint)] and #![feature(const_num_midpoint)]

This is a tracking issue for the midpoint function to {u,i}{8,16,32,64,128,size}, NonZeroU{8,16,32,64,size} and f{32,64}.

The midpoint function calculates the middle point of lhs and rhs.

Public API

impl {u,i}{8,16,32,64,128,size} {
    pub const fn midpoint(self, rhs: Self) -> Self;
}

impl NonZeroU{8,16,32,64,size} {
    pub const fn midpoint(self, rhs: Self) -> Self;
}

impl f{32,64} {
    pub fn midpoint(self, rhs: Self) -> Self;
}

Steps / History

• Implementation: Add midpoint function for all integers and floating numbers #92048
• Final comment period (FCP)¹
• Stabilization PR

Unresolved Questions

• Rounding of negative signed integer (current behavior is round towards negative infinity): Add midpoint function for all integers and floating numbers #92048 (comment)
• Whether it's commutative (some people want midpoint(a, b) == midpoint(b, a), but others want a non-commutative definition)

Footnotes

1. https://std-dev-guide.rust-lang.org/feature-lifecycle/stabilization.html ↩

[barryfam]

If this is still open for comments, I'd like to support rounding towards the first argument.

In binary search, a.midpoint(b) will eventually be called where a and b are adjacent integers. If b is an exclusive range bound, it is helpful to know the function will always return a and not b. [insert Zeno's Paradox meme]

Another reason is for symmetry of step size; see CppCon 2019 talk on C++'s std::midpoint, question at 1:00:00:

Q: [Why is midpoint(a, b) different from midpoint(b, a)?]

A: Partially it's, when you think about it, you're going midway from a to b or midway from b to a. But — I suspect that this comes from the people, who were interested in the beginning — if you're doing simulations and you're stepping around a space where you're measuring things as you step, you start at 1,000 and you step halfway over and over again, you want those steps to be consistent. It doesn't matter which direction you're coming from. If you're coming from +1,000 toward 0 or -1,000 up to 0, you want those steps to be the same size whether you're going counting down or counting up.

In other words, they want:

midpoint(7, 0) = 4    // vs. "round to -∞" would be 3
midpoint(-7, 0) = -4

[scottmcm]

In binary search, a.midpoint(b) will eventually be called where a and b are adjacent integers.

I don't find this part persuasive, as ever binary search I've ever seen has had a ≤ b as an invariant, and thus round-to-first and round-to-negative-infinity are the same.

you want those steps to be consistent

Can you elaborate on the value of "consistent" here? I would have thought that if one needed predictable and consistent steps, they'd use a power-of-two number of quanta, since otherwise each step isn't half the size of the previous one.

---

I guess this boils down to midpoint(a, b) == midpoint(b, a) and midpoint(-a, -b) == -midpoint(a, b) both being "obvious" expectations, but they're incompatible with each other. (Kind of like how -a / -b = -(a / b) and a % b ε [0, |b|) are also "obvious" but incompatible.)

I'm personally still of the opinion that being consistent with (a + b) / 2 (for unsigned values where that doesn't overflow) is the way to go. After all, even P0811R3 used that as the example of how people would write it if that didn't overflow. And it sure would be nice to let clippy say "replace (a + b) / 2 with midpoint" without needing to warn that it might change the behaviour in non-overflow cases too.

---

And since it was in the PR, not the tracking issue, cc #92048 (comment) where I argued for rounding to -∞ instead of to first.

---

Also, FWIW, cranelift has https://docs.rs/cranelift/latest/cranelift/prelude/trait.InstBuilder.html#method.avg_round for (a + b + 1) >> 1 without overflow, which is rounding to +∞. I'm not sure why it picked that one, but if it has it from WASM that might be an interesting proposal to explore for rationale too.

---

EDIT: I happened to stumble on this post accidentally today https://devblogs.microsoft.com/oldnewthing/20220207-00/?p=106223, which is also about being consistent with (a + b)/2.

[bluebear94]

For x86 targets, wouldn’t it be possible to implement {u32, u64}::num_midpoint using the RCR instruction?

[scottmcm]

@bluebear94 Exactly what assembly instructions is more up to LLVM than to us. For now the SHLD approach it gives us is pretty good -- only one instruction more -- someone could always open an LLVM bug for a different instruction later.

[RustyYato]

Why doesn't f32::midpoint upcast to f64 to calculate the midpoint? This should give the same results (excluding NaNs, which may have different bit values).
This should enable a very simple branchless implementation, Godbolt

fn midpoint(a: f32, b: f32) {
    ((a as f64) + (b as f64) / 2.0) as f32
}

I checked this with kani, and for all finite values it is bit-wise identical to the current implementation and they both output NaNs at the same time.

NOTE: that there aren't any calls to kani::assume, so this works for every pair of f32.

#[kani::proof]
fn test_midpoint() {
    let a = kani::any::<f32>();
    let b = kani::any::<f32>();

    // my proposed implementation
    let c = (a as f64 + b as f64) / 2.0;
    let c = c as f32;

    let d = a.midpoint(b); // the current std implementation

    if c.is_nan() {
        assert!(d.is_nan())
    } else {
        assert_eq!(c, d)
    }
}

Run kani with this command
NOTE: this needs --no-overflow-checks because by default kani checks that addition doesn't cause NaNs, but we don't care about that here and there isn't any possibility of overflow anyways.

cargo kani --no-overflow-checks

[Urgau]

@RustyYato As noted in the implementation PR, f64 and f32 follows the libcxx implementation.

This should enable a very simple branchless implementation, Godbolt

Feel free to send a PR changing the implementation (code source).

[RustyYato]

Alright, I've put up #121062, 🤞

[tgross35]

rust/library/core/src/num/f64.rs

Lines 1052 to 1071 in e35364a

	const LO: f64 = f64::MIN_POSITIVE * 2.;
	const HI: f64 = f64::MAX / 2.;
	let (a, b) = (self, other);
	let abs_a = a.abs_private();
	let abs_b = b.abs_private();
	if abs_a <= HI && abs_b <= HI {
	// Overflow is impossible
	(a + b) / 2.
	} else if abs_a < LO {
	// Not safe to halve a
	a + (b / 2.)
	} else if abs_b < LO {
	// Not safe to halve b
	(a / 2.) + b
	} else {
	// Not safe to halve a and b
	(a / 2.) + (b / 2.)
	}

I think the first condition could have an || a.is_sign_positive() !- b.is_sign_positive(), since this is another never overflowing case.

Also, MIN_POSITIVE is the min positive normal number, like in the C++ implementation, but we might be able to adjust this to f32::From_bits(0b10) (2 * min subnormal) since Rust shouldn't flush subnormals (cc rust-lang/rfcs#3514).

The default branch's comment is incorrect too, that operation is the safe one. I have a commit to update it in #127027.

[ronnodas]

I guess this boils down to midpoint(a, b) == midpoint(b, a) and midpoint(-a, -b) == -midpoint(a, b) both being "obvious" expectations, but they're incompatible with each other. (Kind of like how -a / -b = -(a / b) and a % b ε [0, |b|) are also "obvious" but incompatible.)

Round-towards-zero (or round-away-from-zero) satisfies both of those properties so they're not incompatible.