Simplify {f16, f32, f64, f128}::midpoint()

[Jules-Bertholet]

(float_ty::MAX / 2) - (float_ty::MIN_POSITIVE * 2) equals (float_ty::MAX / 2) + (float_ty::MIN_POSITIVE * 2) equals (float_ty::MAX / 2). So these branches are pointless.

CC @Urgau who wrote the original implementation in #92048; does this seem right?

@rustbot label A-floating-point

[rustbot]

r? @jhpratt

rustbot has assigned @jhpratt.
They will have a look at your PR within the next two weeks and either review your PR or reassign to another reviewer.

Use r? to explicitly pick a reviewer

[tgross35]

This makes sense to me but I'll wait for Urgau to double check. I assume the original libcxx version may have been accounting for ftz or something.

One extra nice thing about f16 is it's feasible to test 2 to 3-dimensional operations exhaustively (passed for me):

#![feature(f16)]

pub const fn midpoint1(a: f16, b: f16) -> f16 {
    const LO: f16 = f16::MIN_POSITIVE * 2.;
    const HI: f16 = f16::MAX / 2.;
    let abs_a = a.abs();
    let abs_b = b.abs();
    if abs_a <= HI && abs_b <= HI {
        // Overflow is impossible
        (a + b) / 2.
    } else if abs_a < LO {
        // Not safe to halve `a` (would underflow)
        a + (b / 2.)
    } else if abs_b < LO {
        // Not safe to halve `b` (would underflow)
        (a / 2.) + b
    } else {
        // Safe to halve `a` and `b`
        (a / 2.) + (b / 2.)
    }
}

pub const fn midpoint2(a: f16, b: f16) -> f16 {
    const HI: f16 = f16::MAX / 2.;
    let abs_a = a.abs();
    let abs_b = b.abs();
    if abs_a <= HI && abs_b <= HI {
        // Overflow is impossible
        (a + b) / 2.
    } else {
        // Safe to halve `a` and `b`
        (a / 2.) + (b / 2.)
    }
}

fn main() {
    for i in 0..=u16::MAX {
        for j in 0..=u16::MAX {
            let a = f16::from_bits(i);
            let b = f16::from_bits(j);
            let m1 = midpoint1(a, b);
            let m2 = midpoint2(a, b);
            assert_eq!(m1.to_bits(), m2.to_bits());
        }
    }
}

View changes since this review

[GrigorenkoPV]

Tested f16 (took 2s), f32 (took 13s), and f64 (took 101s) with kani.

#![feature(f128, f16)]

type F = f64;

pub fn old(a: F, b: F) -> F {
    const LO: F = F::MIN_POSITIVE * 2.;
    const HI: F = F::MAX / 2.;

    let abs_a = a.abs();
    let abs_b = b.abs();

    if abs_a <= HI && abs_b <= HI {
        // Overflow is impossible
        (a + b) / 2.
    } else if abs_a < LO {
        // Not safe to halve `a` (would underflow)
        a + (b / 2.)
    } else if abs_b < LO {
        // Not safe to halve `b` (would underflow)
        (a / 2.) + b
    } else {
        // Safe to halve `a` and `b`
        (a / 2.) + (b / 2.)
    }
}

pub fn new(a: F, b: F) -> F {
    const HI: F = F::MAX / 2.;

    let abs_a = a.abs();
    let abs_b = b.abs();

    if abs_a <= HI && abs_b <= HI {
        // Overflow is impossible
        (a + b) / 2.
    } else {
        (a / 2.) + (b / 2.)
    }
}

#[cfg(kani)]
#[kani::proof]
fn equiv() {
    let a: F = kani::any();
    let b: F = kani::any();
    kani::assume(!(a == -b && a.abs() == F::INFINITY));
    let prev = old(a, b);
    let next = new(a, b);
    assert_eq!(prev.is_nan(), next.is_nan());
    if !prev.is_nan() && !next.is_nan() {
        assert_eq!(prev, next)
    }
}

It's telling me the implementations are equivalent. The only thing it complained about was "addition with NaN" (though I believe that is well-defined and legal) which occurs when the arguments are the opposite infinities. Though, both implementations return NaN in that case, from what I can tell.

Currently testing f128. I guess it is going to take around 1000s, so see y'all in 15min or so.

One thing I've noticed during the testing is that f128 seems to lack a Display impl, but it has a Debug impl. Is that intentional?

---

UPD: f128 succeeded too, in mere 807s.

[Jules-Bertholet]

One thing I've noticed during the testing is that f128 seems to lack a Display impl, but it has a Debug impl. Is that intentional?

It’s a known TODO. #116909 (“Figure out f128 which will not be straightforward.”)

[Urgau]

Looks good to me as well.

However, I also wonder why libcxx implementation was and still does that. Perhaps to save a division for performance reasons?

View changes since this review

[Jules-Bertholet]

Perhaps to save a division for performance reasons?

Perhaps. But in that case, they/we could skip the addition as well.

[quaternic]

I assume the original libcxx version may have been accounting for ftz or something.

I think it's pretty much this, or more generally supporting non-default floating point environments. The condition x.abs() < LO is exactly when x / 2.0 is either zero or subnormal. The original implementation avoids causing a spurious FE_UNDERFLOW in the latter case. In those cases the addition like a + (b / 2.0) is then used to nudge the result in the correct direction according to the current rounding mode and set FE_INEXACT as appropriate.

Rust doesn't support any of that so the patch seems fine to me.

[tgross35]

One thing I've noticed during the testing is that f128 seems to lack a Display impl, but it has a Debug impl. Is that intentional?

It’s a known TODO. #116909 (“Figure out f128 which will not be straightforward.”)

For a bit more on this, we can't use the generic implementation used for other types because it uses u64s and we'd pay a cost there. It's on my todo list but a much bigger change than f16 Display, which was a pretty big change :)

[tgross35]

@bors r+

[bors]

📌 Commit bc17bcd has been approved by tgross35

It is now in the queue for this repository.