Floating point round function gives some inaccurate results

[raphlinus]

This bug is branched from rust-lang/rust#55107.

The heart of the implementation of round() in src/library.js is:

   return d >= +0 ? +Math_floor(d + +0.5) : +Math_ceil(d - +0.5);

This produces incorrect values for 0.49999999999999994 (correct 0.0, actual 1.0) and 4503599627370497.0 (correct same as input, actual 4503599627370498.0) among others. The latter is disturbing because it means the function does not reliably preserve integers.

I worked out a fix, which is basically this in the positive branch:

        Math_floor((d + (0.25 - 0.5 * F64_EPSILON)) + (0.25 + 0.5 * F64_EPSILON))

where F64_EPSILON is is the difference between 1.0 and the next largest representable number, or 2.2204460492503131e-16. At least this gives correct results on my x86_64; I'm not sure to what extent correctness is guaranteed across all target platforms.

Also, the ternary might profitably be replaced by (this is C syntax, because I know the libm names of things but not necessarily the relevant library functions/intrinsics in Emscripten):

    copysign(floor(((abs(d) + (0.25 - 0.5 * F64_EPSILON)) + (0.25 + 0.5 * F64_EPSILON)), d)

This will probably be faster in wasm because there is a copysign intrinsic. But if the intrinsic doesn't get compiled properly it's probably an overall lose.

I'm willing to work on this, but wanted to give people a heads-up, and maybe could get some guidance on the copysign question.

[kripken]

Interesting how many corner cases there are with this stuff ;)

May be good to look at musl's implementation,

https://github.com/kripken/emscripten/blob/incoming/system/lib/libc/musl/src/math/round.c

We can do basically what they are, probably. In fact, we may want to just build round from musl - right now we avoid doing that,

https://github.com/kripken/emscripten/blob/incoming/tools/system_libs.py#L168

We avoided it because the JS version is smaller, but if it's incorrect, maybe we should just use musl.

[raphlinus]

So the previous version did use musl and this was changed to an "own" implementation for speed, see fecf842. Basically, that commit was inadequately attentive to correctness. On the other hand, the musl version is pretty slow and leads to significant code bloat - I think it's something like 92 wasm instructions, compiled.

I'm proposing that my version is better. I checked the 32 bit version against all possible 32-bit bit patterns, and since filing the bug read up on IEEE 754 rounding guarantees some more, so not am confident it's correct on all compliant implementations. It's only one 32-bit addition more than the current master. And with copysign, the generated code can be even shorter, and eliminate the branch.

[kripken]

Interesting. Sounds promising, but we'd need careful review for code like that.

About the copysign issue, to benefit from it we may write it in C and use an intrinsic (or, write it in C and write inline wasm).

[raphlinus]

Yes, I would welcome more eyes on this. Already the Julia people are looking at the idea (JuliaLang/julia#29700), and I've also proposed it to https://github.com/japaric/libm, so I think it's likely it will get analyzed more carefully.

I pushed a tentative PR with the idea, which seems like it will use the llvm intrinsic where appropriate, but I haven't verified that, or written tests. I have some questions, such as whether it's possible to have the computations done in 32-bit floats (the ones in library.js seem to basically imply 64 bits). I'd be quite willing to believe that writing it in C is the solution to that.

I doubt inline wasm is necessary, if we get the intrinsic selected, I'm sure it'll compile to good code.

[VirtualTim]

Apparently pretty much every javascript engine made the same mistake at some point (plus the implementation in the spec did not implement the spec). https://ridiculousfish.com/blog/posts/JavaScript-Rounding.html

[simonbyrne]

A faster option might might be trunc(x + copysign(prevfloat(0.5), x)) (this is Julia code, but hopefully you get the idea).

[raphlinus]

@simonbyrne Indeed, thanks for that suggestion! I updated the PR and also fixed a dumb typo.

[stale]

This issue has been automatically marked as stale because there has been no activity in the past year. It will be closed automatically if no further activity occurs in the next 7 days. Feel free to re-open at any time if this issue is still relevant.

[sbc100]

Seems like this might still be an issue for people?

[stale]

This issue has been automatically marked as stale because there has been no activity in the past year. It will be closed automatically if no further activity occurs in the next 30 days. Feel free to re-open at any time if this issue is still relevant.