Fix testcases on x87

[plugwash]

The tests for the libm crate are failing on debian i386 and this is preventing the rust-libm package from migrating to testing. Debian i386 uses the x87 floating point unit.

---- math::ceil::tests::sanity_check stdout ----
thread 'math::ceil::tests::sanity_check' panicked at 'assertion failed: `(left == right)`
  left: `1.10009765625`,
 right: `2.0`', src/math/ceil.rs:50:9
stack backtrace:
   0: rust_begin_unwind
   1: core::panicking::panic_fmt
   2: libm::math::ceil::tests::sanity_check
             at ./src/math/ceil.rs:50
   3: libm::math::ceil::tests::sanity_check::{{closure}}
             at ./src/math/ceil.rs:49
   4: core::ops::function::FnOnce::call_once
             at /usr/src/rustc-1.47.0/library/core/src/ops/function.rs:227
note: Some details are omitted, run with `RUST_BACKTRACE=full` for a verbose backtrace.

---- math::floor::tests::sanity_check stdout ----
thread 'math::floor::tests::sanity_check' panicked at 'assertion failed: `(left == right)`
  left: `0.10009765625`,
 right: `1.0`', src/math/floor.rs:49:9
stack backtrace:
   0: rust_begin_unwind
   1: core::panicking::panic_fmt
   2: libm::math::floor::tests::sanity_check
             at ./src/math/floor.rs:49
   3: libm::math::floor::tests::sanity_check::{{closure}}
             at ./src/math/floor.rs:48
   4: core::ops::function::FnOnce::call_once
             at /usr/src/rustc-1.47.0/library/core/src/ops/function.rs:227
note: Some details are omitted, run with `RUST_BACKTRACE=full` for a verbose backtrace.

---- math::j1f::tests::test_y1f_2002 stdout ----
thread 'math::j1f::tests::test_y1f_2002' panicked at 'assertion failed: `(left == right)`
  left: `-0.10703231`,
 right: `-0.10703229`', src/math/j1f.rs:370:9
stack backtrace:
   0: rust_begin_unwind
   1: core::panicking::panic_fmt
   2: libm::math::j1f::tests::test_y1f_2002
             at ./src/math/j1f.rs:370
   3: libm::math::j1f::tests::test_y1f_2002::{{closure}}
             at ./src/math/j1f.rs:369
   4: core::ops::function::FnOnce::call_once
             at /usr/src/rustc-1.47.0/library/core/src/ops/function.rs:227
note: Some details are omitted, run with `RUST_BACKTRACE=full` for a verbose backtrace.

---- math::sincosf::tests::rotational_symmetry stdout ----
thread 'math::sincosf::tests::rotational_symmetry' panicked at 'assertion failed: (s - s_plus).abs() < TOLERANCE', src/math/sincosf.rs:148:13
stack backtrace:
   0: rust_begin_unwind
   1: core::panicking::panic_fmt
   2: core::panicking::panic
   3: libm::math::sincosf::tests::rotational_symmetry
             at ./src/math/sincosf.rs:148
   4: libm::math::sincosf::tests::rotational_symmetry::{{closure}}
             at ./src/math/sincosf.rs:138
   5: core::ops::function::FnOnce::call_once
             at /usr/src/rustc-1.47.0/library/core/src/ops/function.rs:227
note: Some details are omitted, run with `RUST_BACKTRACE=full` for a verbose backtrace.

The x87 unlike most modern FPUs uses a different representation internally (in floating point registers) than externally. This means that intermediates stored in floating point registers can have more precision than the programmer intended ("excess precision"). Most of the time this just results in slightly more accurate than intended results, but sometimes when code uses clever tricks it can lead to major breakage. I believe that this excess precision is the cause of the test failures.

This pull request represents a first stab at dealing with these issues, I am a distribution maintainer who deals with issues in a whole bunch of languages not a rust specialist, so it's probable that the code does not follow best practices.

The failure in j1f does not look serious to me and I merely adjusted the test to allow the marginally different result.

The failures in floor and ceil look far more serious, the values returned from those functions should always be integers.

The test errors for the rotational_symmetry test did not reveal whether the failure was serious or just a marginal difference, so I adjusted the test code a bit and confirmed that it was also a major failure.

The root cause of all of these failures is a trick used in the foor, ceil, round, roundf and rem_pio2f (but not floorf or ceilf) functions. A large number is added and subtracted (or subtracted and added in some cases) to force rounding to an integer value. I am 99% sure that lack of rounding in this process is the cause of the faliures.

For the floor and ceil functions I modified the functions to use an alternative implementation as proposed in #219

For the round, roundf and rem_pio2f I did not have any alternative implementations readilly available, so I tried converting the floating point numbers to bit patterns and back. This seems to work correctly, though I cannot seem to find any documententation on the rust compiler's behaviour concerning excess precision.

[Lokathor]

This is mostly good.

There's a bunch of CI errors, I think you just need to make the use statements conditional.

[Lokathor]

This could be added on a file by file basis, but it shouldn't be at the crate root.

[plugwash]

I've updated the pull request to move the unreachable code declarations into the individual source files, fix the unused import and cargo fmt issues and added comments to explain why the code is using to_bits followed immediately by from_bits.

However there still seem to be test failures being reported by the CI that don't seem to have anything to do with my changes, are those just expected?

[Lokathor]

I'd like to say that this is unrelated test failure and thus it's an acceptable test failure, however i am the junor-most member of the maintainer team here.

The two maintainers who are most likely to know are @alexcrichton and @japaric , but they are both very busy people, so I'm gonna bring this up to T-libs in Zulip as well.

[m-ou-se]

Is the goal of this hack just to print the value when this test fails? (If so, adding the message to the assert macro would be better: assert!((s - s_plus).abs() < TOLERANCE, "|{} - {}| = {} > {}", s, s_plus, (s - s_plus).abs(), TOLERANCE);.)

[m-ou-se]

How about x86 targets where the x87 fpu is not used? Should those platforms use the original or the new code?

(Maybe this should be something like #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))] instead.)

[m-ou-se]

Tip: You can put the use statement inside the {} block below. That way, you only need one #[cfg(..)] and don't need to keep them in sync.

[m-ou-se]

I'd like to say that this is unrelated test failure and thus it's an acceptable test failure

Definitely looks unrelated. The same failures happen locally both with and without this change.

The failure is potentially concerning though. Maybe a bug got introduced (or fixed!) in musl or rustc or llvm or something like that.

To unblock this PR, we could for now just merge a PR that comments out the failing test cases, as they were already failing anyway. Then open an issue to investigate those asap.

[plugwash]

Thanks for the review m-ou-se

Yes the goal of the stuff around the assert was to find out how badly wrong the results were, I originally tried to manually add a message to the assert, but floundered due to lack of rust knowledge. I've changed it to use the technique you proposed.

I agree that it makes sense to use the same implementations on x86 with sse2 as on x86-64, whether the old or new implementation is a better choice than the current one for platforms other than x87 is not a decision i am in a position to make.

Thanks also for the simplification tip on the import stuff.

I've updated the branch with the changes.

[plugwash]

Unfortunately I have further bad news.

First off I discovered from a testcase in num_complex that rem_pio2 is also broken. I made a patch for that too https://salsa.debian.org/rust-team/debcargo-conf/-/blob/master/src/libm/debian/patches/fix-rem-pio2.patch

But what is worse I discovered that the trick of converting the float to bits and back again to force conversion to the storage format doesn't work in release mode.

[Amanieu]

You can try using the force_eval! macro to force the compiler to write a floating-point value to memory (as a f32/f64) and then read it back. This uses a volatile operation which guarantees that the compiler won't optimize it away.

[plugwash]

Ok, I had another go at this.

The first thing I did was rebase my PR onto master. In the process I dropped the changes to the round and roundf functions as these have now been re-implemented in terms of floor and floorf.

I then replaced to_bits/from_bits with force_eval in rem_pio2f and confirmed that it was sufficient to force rounding to the storage format, even in release mode. I do feel this is a bit of a hack, but until/unless rust gets well-defined floating point semantics it may be the best we can do. I placed the force_evals behind an architecture check so that they do not cause pessimisation of the code on architectures that don't suffer from excess precision.

I then moved on to the issue with rem_pio2, which had been revealed by a test in num_complex. I applied a fix to rem_pio2 based on the one I had already made to rem_pio2f. I also added tests to sincos based on those from sincosf.

I then went through several tests which in release mode were reporting inequalities between numbers that looked equal when printed. I added some calls to force_eval! to force rounding to the storage format before comparison and these tests passed.

That left all tests passing with a plain "cargo test" and one test failing with "cargo test --release". The one remaining test was
math::rem_pio2::tests::test_near_pi , it failed with different values in the third component of the result.

What puzzled me is that the test passed before I made my changes, despite the fact that there was clearly horrible breakage in rem_pio2, I started to suspect that the calls to rem_pio2 in test_near_pi were being evaluated entirely at compile time and hence not suffering from excess precison. To test this I modified test_near_pi to use force_eval on the values before passing them to rem_pio2. I placed the commit with this change first in the commit sequence (despite it being the last commit written) so it could easily be tested with both the old and the new rem_pio2. With this change, the test failed in the same way both before and after my changes.

I couldn't really understand the code in rem_pio2, so I got someone I know who is much better at maths than me to look at it. After some discussion we came to the conclusion that it was likely doing a form of "double double" arithmetic, and that this was being affected by excess presision. It's not clear to me though if this is of wider consequence, it may well be that once rem_pio2 is inlined into a sin, cos or sincos function and optimized that the "double double" artithmetic essentially becomes "double extended" arithmetic which works just fine.

Overall while this PR may not be perfect, I believe it represents a substantial improvement over the status quo on x87. On other architectures it should be a no-op.

Strangely there seems to be a CI failure on powerpc64el though, I can't see how that can be caused by my changes given the way I put conditionals round them but when I look for test results for the master branch they seem to also be failures, though on different architectures (annoyingly the CI seems to abort testing when any one architecture fails which does not make comparisions easy).

[Amanieu]

I have a pending PR #259 which fixes CI. Once that is merged, rebase your PR so CI passes.

[Amanieu]

Thanks!

[lopopolo]

@Amanieu @plugwash with this PR merged, does that mean #243 is fixed? Should that issue be closed?