Test against MPFR

[tgross35]

Test against MPFR via rug for infinite precision results in cases where musl might not be accurate.

[tgross35]

I think this is ready for a review too @Amanieu. This adds tests against MPFR for almost all functions, which gives us infinite precision results in cases where we want to be more accurate than musl.

The few untested functions are at

libm/crates/libm-test/tests/multiprecision.rs

Lines 47 to 61 in 4d51e47

	skip: [
	// FIXME: MPFR tests needed
	frexp,
	frexpf,
	ilogb,
	ilogbf,
	ldexp,
	ldexpf,
	modf,
	modff,
	remquo,
	remquof,
	scalbn,
	scalbnf,
	],

.

This targets my test-refactoring branch so it can't be merged directly (requires #300).

[tgross35]

Cc also other @tspiteri (rug maintainer) if you get a chance to take a look at this.

[beetrees]

I don't think subnormalize_ieee_round with Round::Nearest will work correctly for next_toward: e.g. for an f32 with bit repr 0x1 (the smallest representable positive non-zero f32), the correct result for next_toward positive infinity would be the f32 with bit repr 0x2. However, the rug next_toward method will calculate the next value as if the input is a normal number, meaning that subnormalize will round away the difference, meaning the MPFR result will be incorrect. The easiest way to fix this that I can think of is to have subnormalize_ieee_round always round the result in the direction next_toward is trying to go in (so use Round::Up if self.1 > self.0 and Round::Down if self.1 < self.0).

[tgross35]

This is a good catch. That sounds reasonable but I think I should probably add a unit test that the MPFR behavior is correct. I think maybe it is better to have that as a followup, so I just removed this test and left a FIXME.

[beetrees]

AFAIK, since MPFR discards NaN bits, this will also fail when input.0 is NaN (and similarly for fabs). This won't have been causing any tests to fail yet due to #300 (comment), but will do once that is fixed.

[tgross35]

Even testing locally with a much higher iteration count, I haven't seen this fail at all. I guess I will just leave it as-is, unless it turns out to be a problem later.

[tgross35]

The docs say this works "even when op1 or op2 is a NaN" https://www.mpfr.org/mpfr-current/mpfr.html and the code says:

   The computation of z with magnitude of x and sign of y:
   z = (-1)^signbit(y) * abs(x), i.e. with the same sign bit as y,
   even if z is a NaN.

The code does say https://gitlab.inria.fr/mpfr/mpfr/-/blob/0b0ce34f301bf75f0a4160350a179de19d9a5146/src/copysign.c#L24-30

The result does still seem surprising.

[tgross35]

The MPFR jn implementation on i686 seems to be hitting this assertion https://gitlab.inria.fr/mpfr/mpfr/-/blob/0b0ce34f301bf75f0a4160350a179de19d9a5146/src/jn.c#L250-251

[tgross35]

That assertion might actually be the better behavior, seems like large numbers of iterations can't complete in reasonable time with MPFR's jn. I assume carrying full precision through large numbers of iterations is significantly costlier than whatever approximation we are doing.

Will skip these for now.

[tgross35]

I think the current state of this PR should be a good enough start to getting a better test basis, even if there are some rough edges still. I am going to merge in order to unblock any refactoring that builds on it, as well as to help out with the implementation changes that would benefit from the testing.

Please feel free to add further reviews here, or open issues/PRs if there are problems.