std::round and friends

Overview of rounding

Rounding modes (floating-point)

¹ new in C23, not in C++

60559 default rounding mode

• Avoids bias that any other rounding mode would introduce.
• Did you ever use anything else?

Rounding functions

The deal with nearbyint and rint

You get the dynamic rounding mode:

• Your function is dependent on global state (local to each CPU core).
• The result could be different in different threads of your program.

Upcoming rounding functions (roundeven)

It's already in C23 and GCC also implements it as a builtin (__builtin_roundeven)

Ideas for implementing rounding (1)

float r0(float x) {
  int i = x + .5f;
  return i;
}

r0(1.5f) = 2
r0(1.4f) = 1
r0(0.5f) = 1
r0(0.4f) = 0
r0(-0.4f) = 0
r0(-0.5f) = 0
r0(-1.4f) = 0
r0(-1.5f) = -1
r0(4e9f) = UB 

What rounding function is this?

Ideas for implementing rounding (2)

using L = std::numeric_limits<float>;
float r1(float x) {
  constexpr float imant = 1 << (L::digits - 1);
  return (x + imant) - imant;
}

r1(1.5f) = 2
r1(1.4f) = 1
r1(0.5f) = 0
r1(0.4f) = 0
r1(-0.4f) = -0.5
r1(-0.5f) = -0.5
r1(-1.4f) = -1.5
r1(-1.5f) = -1.5
r1(4e9f) = 4e9f

Rounding the way we learned in school

std::round(2.5)

But how can this be implemented?

No or little hardware support (x86)

• Floating-point rounding modes didn't support roundTiesToAway for the longest time
• Consequently, no (x + a) - a could ever implement std::round
• x87 has a round to integral instruction (basically rint)
• SSE4.1 added round?? instructions: all rounding functions except round
• glibc's libm implements round via bit-manipulation

Idea: trunc(x ± 0.5) should do it

• With -march=x86-64-v2, trunc becomes a single roundss instruction.
• Without it we can use float->int conversion (after an input range branch)

float round(float x) {
  float half = copysign(0.5f, x);
  return trunc(x + half);
}

Idea: trunc(x ± 0.5) almost did it, let's adjust

float round(float x) {
  constexpr float almost_half = std::nextafter(0.5f, 0.f);
  return trunc(x + copysign(almost_half, x));
}

Now:

• round(0x1p23f ) = 0x1p23f
• round(0x1p23f+1) = 0x1p23f+1
• round(0.5f) = trunc(0x1p-1f + 0x1.fffffep-2f)
  with roundTiesToEven we get trunc(1.f)

GCC doesn't emit trunc(x ± 0.4999…), why?

https://compiler-explorer.com/z/dGPhxMY3T

• Clang does, so what's up here?

• GCC assumes you care about floating-point exceptions. Clang assumes you don't care.

• I believe Clang has the better default.

• Anyhow, what do floating-point exceptions have to do with round?

FE_INEXACT — who cares?

• Floating-point addition x + 0x1.fffffep-2f is likely going to be inexact.
• Inexact result raise FE_INEXACT.
• std::round is specified to not raise any floating-point exceptions except for SNaN.

How does glibc implement roundf?

• It's integer-based bit-manipulation only:
  • independent from rounding modes
  • never going to trap
• It's actually surprisingly efficient.
  • My benchmarks: calling roundf is more efficient than addps; roundss

• The glibc implementation is fast thanks to branches.
• A branchless implementation is not the most efficient one.

Rounding without trapping (1)

Let's try a different approach:

return copysign(trunc(abs(x)) + fixup, x);

• All operations except for the addition won't trap.
• Can we generate a value for fixup that will avoid all traps?

Rounding without trapping (2)

t = trunc(abs(x));
fixup = abs(x) - t >= .5 ? 1 : 0;
return copysign(t + fixup, x);

This works and doesn't trap.

Except for :

• abs(x) and t are both ,
• leading to ,
• which traps and gives us NaN <= .5.

Rounding without trapping (3)

if (!isinf(x)) {
    t = trunc(abs(x));
    fixup = abs(x) - t >= .5 ? 1 : 0;
    return copysign(t + fixup, x);
} else {
    return x;
}

Not bad, but can we go faster and without a branch (for SIMD, anyway)?

Rounding without trapping (4)

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=127119

• I proposed a integer-based fixup after trunc(abs(x)).
• Hongyu Wang proposed another solution.

Integer-based fixup

template <__vec_builtin V>
  V
  __round(V __x)
  {
    using T = __vec_value_type<V>;
    using U = _UInt<sizeof(T)>;
    using I = __integer_from<sizeof(T)>;
    const auto abs_x = __fabs<_Traits>(__x);
    const auto t_abs = __trunc<_Traits>(abs_x);
    const auto int_abs_x  = __vec_bit_cast<I>(abs_x);
    const auto int_t_abs  = __vec_bit_cast<I>(t_abs);

    // a) for |x| without fractional mantissa bits: diff = 0
    // b) for |x| >= 1: XOR cancels exponent bits and mantissa bits signifying values >= 1
    //                  If the MSB is at the .5 position we need to add 1 to the result
    // c) for |x| < 1: trunc(|x|) is 0 => so diff = bit-pattern of |x|
    // => diff needs to be compared against
    // a) any positive non-zero integer
    // b) integer with 1 bit at .5 position
    // c) bit-pattern of .5
    const auto diff = int_abs_x ^ int_t_abs;

    constexpr int mant_width = numeric_limits<T>::digits - 1;
    constexpr I one_bits = std::bit_cast<I>(T(1));
    constexpr I one_half_bits = std::bit_cast<I>(T(.5));
    constexpr I bias = one_bits >> mant_width;

    const auto biased_exp = int_abs_x >> mant_width;
    const auto exponent = __vec_bit_cast<I>(__vec_bit_cast<U>(biased_exp - bias) % numeric_limits<U>::digits);
    const auto threshold = ((biased_exp > (mant_width + bias)) | (int_abs_x < one_bits))
        ? one_half_bits                                 // a) and c)
        : (I(1) << (mant_width - 1)) >> exponent;   // b)

    const V sign_bit = __vec_xor(abs_x, __x);
    const V r_abs = t_abs + std::bit_cast<V>(diff < threshold ? 0 : one_bits);
    return __vec_or(sign_bit, r_abs);
  }

Discuss