float rounding is slow

[raphlinus]

Update (2025-02-03)

The semantics of round will not be changed. However, the docs for round should point out that on most hardware, round_ties_even is faster.

Original issue

The scalar fallback for the sinewave benchmark in fearless_simd is very slow as of the current commit, and the reason is the f32::round() operation. When that's changed to (x + 0.5).floor() it goes from 1622ns to 347ns, and 205ns with target_cpu=haswell. With default x86_64 cpu, floorf() is a function call, but it's an efficient one. The asm of roundf() that I looked at was very unoptimized (it moved the float value into int registers and did bit fiddling there). In addition, round() doesn't get auto-vectorized, but floor() does.

I think there's a rich and sordid history behind this. The C standard library has 3 different functions for rounding: round, rint, and nearbyint. Of these, the first rounds values with a 0.5 fraction away from zero, and the other two use the stateful rounding direction mode. This last is arguably a wart on C and it's a good thing the idea doesn't exist in Rust. In any case, the default value is FE_TONEAREST, which rounds these values to the nearest even integer (see Gnu libc documentation and Wikipedia; the latter does a reasonably good job of motivating why you'd want to do this, the tl;dr is that it avoids some biases).

The implementation of f32::floor is usually intrinsics::floorf32 (but it's intrinsics::floorf64 on msvc, for reasons described there). That in turn is llvm.floor.f32. Generally the other round functions are similar, til it gets to llvm. Inside llvm, one piece of evidence that "round" is special is that it's not listed in the list of instrinsics that get auto-vectorized.

Neither the C standard library nor llvm intrinsics have a function that rounds with "round half to even" behavior. This is arguably a misfeature. A case can be made that Rust should have this function; in cases where a recent Intel CPU is set as target_cpu or target_feature, it compiles to roundps $8 (analogous to $9 and $a for floor and ceil, respectively), and in compatibility mode the asm shouldn't be any slower than the existing code. I haven't investigated non-x86 architectures though.

For signal processing (the main use case of fearless_simd) I don't care much about the details of rounding of exactly 0.5 fraction values, and just want rounding to be fast. Thus, I think I'll use the _mm_round intrinsics in simd mode (with round half to even behavior) and (x + 0.5).floor() in fallback mode (with round half up behavior). It's not the case now (where I call f32::round) that the rounding behavior matches the SIMD case anyway. If there were a function with "round half to even" behavior, it would match the SIMD, would auto-vectorize well, and would have dramatically better performance with modern target_cpu.

[m1el]

Edit: reworded If you use (x+0.5).floor() to work as a replacement for round function, it will produce incorrect results, not only for exact 0.5 fraction values.

You can test this by enumerating all f32 values (the same works for f64, but not as easily enumerable): https://gist.github.com/m1el/dc18a85b82357589825d9ffa4532b0a6

For example, 0.49999997 + 0.5 == 1.0 and 8388609.0 + 0.5 == 8388610.0, due to loss of precision.

[bstrie]

@raphlinus , I'm actually surprised to discover that round doesn't do half-to-even, I guess this is a case of LLVM's choices inadvertently dictating Rust's semantics (though, reading LLVM's docs, it seems that they chose away-from-zero in order to have identical output to libm (path dependence!)). It's probably too late to change the semantics of round now (though it seems minor enough that if there were very vocal community demand, people might consider it a bugfix), which is a shame since e.g. Python's stdlib rounds half-to-even and I've never heard of anyone complaining about it, and the bias prevention is certainly nice. And if we did have a hypothetical round_unbiased (or w/e, let the bikeshed begin!), it would be a shame to effectively penalize people for reaching for the slower-and-more-obvious round (maybe we ought to consider deprecating round for round_away or smth).

At the same time, people may also find this a good opportunity to implement the rest of the five rounding modes specified by IEEE 754: https://en.wikipedia.org/wiki/IEEE_754#Rounding_rules . Would be ideal if LLVM could be the one handling all of this.

[bstrie]

@m1el I'm afraid I don't see how that's relevant to this issue. Please file a separate bug if you think that the behavior of floor is incorrect.

[m1el]

@bstrie as far as I can tell, the behavior of floor is correct. My point is that the behavior of (x+0.5).floor() as a round function is not.

[raphlinus]

@bstrie Right. I do think there are strong benefits to llvm implementing this; it's a lot easier to make auto-vectorization work, and then all of (round, floor, ceil) will be similar.

@m1el I do appreciate the heads-up, I hadn't considered the fact that the addition would round up on precision loss. I agree it's technically not relevant to this issue but is a good warning to people who might come across this issue. It's also a signal of the importance of fixing this problem, as the obvious workaround is flawed.

Personally I would be fine with the rounding-on-half behavior to change. The number of people who will get hit by the performance and the inconsistency with SIMD is large. The number of people who actively depend on 0.5f32.round() == 1.0 is small; I can imagine it maybe breaking a golden master in one test or two. But I understand any non-backwards compatible changes can be a sensitive point.

[bstrie]

Poking around, there may be some tentative LLVM support for this after all, though it's marked "experimental":

https://llvm.org/docs/LangRef.html#constrained-floating-point-intrinsics

https://llvm.org/docs/LangRef.html#id1888

https://llvm.org/docs/LangRef.html#id1939

@raphlinus , I don't know a lot about how rounding modes work at a low-level, but does that look like it might be a promising start?

[raphlinus]

@bstrie On a first glance, these don't seem to be what we're looking for, it seems like these experimental intrinsics affect rounding mode for the precision loss in other floating point operations. Also, I didn't specifically see round-half-to-even behavior listed.

One thing we should probably do before nailing down new semantics for rounding is research what other chips implement performantly (I'm personally mostly concerned about arm). Want me to look into that?

[varkor]

On a first glance, these don't seem to be what we're looking for, it seems like these experimental intrinsics affect rounding mode for the precision loss in other floating point operations.

llvm.experimental.constrained.rint seems to be the appropriate intrinsic. The LLVM "round.tonearest" mirrors C's FE_TONEAREST, but this appears not to specify behaviour under ties (though most comments I found indicated it was round-to-even). It should be possible to make this explicit in LLVM if we wanted.

[raphlinus]

I'm still not convinced the llvm.experimental.constrained.rint intrinsic is the one we want. Reading the docs, I believe it still compiles to roundps $4, which means it gets the rounding mode from the environment. What I think the extra arguments do is let llvm do optimizations assuming that the environment matches those arguments. I see the point in that but it's a very thin win. I will experiment to make sure though. The phabricator issue introducing these intrinsics makes very good reading.

I verified with llc and hand-editing IR code that llvm.experimental.constrained.rint compiles to roundps $4 and llvm.experimental.constrained.nearbyint compiles to roundps $12. (These are decimal, I had hex above for ceil, apologies.)

We already have intrinsics::nearbyintf32 which generates roundps $12 (or call nearbyintf in fallback mode, which is ok). Btw, I do think nearbyint is the better choice than rint here; as I read it, the only difference is that the latter generates exceptions, and I'm pretty sure we don't want the exceptions.

So it looks like the way forward is to push on llvm for a new intrinsic. Having started reading their stuff, this is pretty deep water; it has to be supported through many optimization and codegen components, and there needs to be a fallback implementation for all targets, because this is not something that's represented anywhere in libm.

Another possibility to explore is to expose nearbyint, with the understanding that it has an implicit dependence on the FP state. This is certainly not much work, but it might be a bad idea.

Lastly, I checked availability of these operations on arm64, and it's much nicer than Intel - the FRINTN instruction exists in both scalar and vector, with "ties to even" behavior, as do instructions for the all the libm rounding functions. So I believe this is a small amount of additional support for the idea that it should be exposed as an intrinsic in llvm.

[raphlinus]

I looked also at wasm; here's what I found. Short story is that current codegen is both wrong (for emscripten) and inefficient.

The code generated for wasm32-unknown-emscripten for round is effectively if x >= 0 { (x + 0.5).floor() } else { (x - 0.5).ceil() }. For f32, the value is converted to f64 then back to f32, and this is correct but slow. For f64, this code will produce incorrect values for 0.49999999999999994 (correct 0.0, actual 1.0) and 4503599627370497.0 (correct same as input, actual 4503599627370498.0) among others. I think the source of this code may be src/library.js from kripken/emscripten, but I'm not 100% sure.

(Note: edited, I thought I had a clever fix, but it was a testing error).

For wasm32-unknown-unknown, the generated code seems to match japaric/libm. I haven't checked this code for correctness, but it looks pretty carefully written and very very slow (it's something like 92 wasm instructions and lots of branches). (Edit: I did verify that this produces identical results as libm, with the exception that std produces negative zero for negative inputs near zero and japaric/libm produces positive zero).

In any case, that's the current situation. The designers of wasm are also a fan of round-half-to-even behavior, and so wasm includes the nearest instruction in both f32 and f64 flavors for that. They don't have an instruction for round-away-from-zero behavior.

Lastly, I found that there is an architecture-specific intrinsic for this already in llvm: int_aarch64_neon_frintn. Based on what I've found, I think the argument that this should be promoted into an llvm-wide intrinsic is strong.