f32::FRAC_PI_6.to_degrees() is now less accurate

[cuviper]

On stable, f32 FRAC_PI_6.to_degrees() returns exactly 30.0. On beta and nightly, it returns 30.000002.

This is probably related to #47919, which intended to increase accuracy. cc @varkor

Here's a full test of PI constants, asserting errors < 1e-6:

use std::f32::consts::*;

const DEG_RAD_PAIRS: [(f32, f32); 7] = [
    (0.0, 0.),
    (22.5, FRAC_PI_8),
    (30.0, FRAC_PI_6),
    (45.0, FRAC_PI_4),
    (60.0, FRAC_PI_3),
    (90.0, FRAC_PI_2),
    (180.0, PI),
];

fn main() {
    for &(deg, rad) in &DEG_RAD_PAIRS {
        let rad2 = deg.to_radians();
        let deg2 = rad.to_degrees();
        
        println!("degrees: {:?} {:?}", deg, deg2);
        assert!((deg2 - deg).abs() < 1e-6);
        
        println!("radians: {:?} {:?}", rad, rad2);
        assert!((rad2 - rad).abs() < 1e-6);
    }
}

playground

It passes on stable 1.24, but fails on beta and nightly.

[varkor]

Ah, this isn't on 😕
It seems like the improvements were not uniform. I'll get this fixed.

[cuviper]

FWIW, I think that 1.0f32 test case, which spurred the change, is so trivial that it's misleading. This is a multiplicative identity, of course, but it's also not representative -- who ever uses exactly 1 radian? Accuracy with the PI constants will be more generally beneficial, even if it sacrifices the unit case. But hey, maybe there's a way to do better at both!

[varkor]

@cuviper: I agree — the identities with π are more important than the 1 radian case. The intention was to provide increased accuracy everywhere, and that was just one demonstrative test case, but clearly that hasn't worked out. I think I have an implementation that works for both these cases, though, so hopefully we can get the best of both worlds! (Floating point is always tricky!)

[Mark-Simulacrum]

We should probably revert the initial patch so as to not regress behavior here until we come to a more principled solution. Could someone prepare that PR for beta?

[varkor]

@Mark-Simulacrum: see the discussion in #48622. I think this this is an intended effect (the accuracy has regressed for some inputs, but overall the accuracy across the range of floating points has increased).

I'll write up a more detailed response later.

[varkor]

My inclination is to close this as expected breakage. I don't think we really explicitly acknowledged the cases where the double-rounding of the previous implementation of f32::to_degrees caused some inputs to be assigned correct outputs, but where the single-rounding of the new implementation caused the input to be off by 1 ULP. It's also unfortunate that some of the inputs for which the function has regressed are on "well-known" mathematical identities.

That said, as #48622 (comment) describes, the overall accuracy for f32::to_degrees has increased as expected: the number of inputs for which the output is off by 1 ULP is roughly 4.5x lower than the previous implementation. I think behaviour over the entire range is more important than particular identities (though if the proportion was close to 1, I would be disposed to keeping the implementation preserving them).

There are ways to fix this completely (e.g. an arbitrary-precision algorithm, or possibly even just using f64, as was the approach in #48622), but they incur costs that are not desirable for the typical user, who just wants a fast radians-to-degrees conversion.

Using to_degrees as if it were correct to 0.5 ULP on all inputs previously was incorrect, and that assumption has not changed (although I would treat an error of more than 1 ULP a bug) — in particular, assuming any particular equality holds exactly is fallacious. (Tangentially, I can't even think where one might even need a maximally accurate to_degrees function, but that's by-the-by.)

Therefore, though I concede it's a little disappointing that those particular cases have regressed, I think the current implementation is the preferred one.

(@rkruppe has proved more prudent on this issue previously, though, so I'd like to hear what they think before any final decision is made.)

[hanna-kruppe]

I generally agree that improvements "across the board" are better than spot fixes on identities. I do have some questions about the data in #48622 (comment) which I now posted there (sorry, I meant to post these earlier but I've been busy).

Regardless what we decide here, I think this isn't a drastic enough regression to justify any quick reverts, since as @varkor said:

Using to_degrees as if it were correct to 0.5 ULP on all inputs previously was incorrect, and that assumption has not changed (although I would treat an error of more than 1 ULP a bug) — in particular, assuming any particular equality holds exactly is fallacious. (Tangentially, I can't even think where one might even need a maximally accurate to_degrees function, but that's by-the-by.)

If we decide to revert, we can do so after careful deliberation. This slipping in a stable release, even if we ultimately revert it, doesn't seem very bad.

[cuviper]

I tried a random probe of inputs from [0, PI) and compared them to f64 on the playground. I acknowledge that f64 can't represent a perfect truth either, but it's better than nothing. The results are consistently in favor of the current method using the constant 57.29... -- like these ULP diffs:

old: {-1: 650138, 0: 349862}
new: {0: 861547, 1: 138453}

[cuviper]

Surprisingly, a quick test of n / (consts::PI / 180.0f32) does even better:

old: {-1: 650259, 0: 349741}
new: {0: 862103, 1: 137897}
test: {0: 908434, 1: 91566}

[varkor]

@cuviper: I can confirm that my tests are saying the same thing: between 0x00004000 and 0x00008000, there are 273 inputs with erroneous output for the current method and only 179 for the division method :)

(It's also pleasant that f32::consts::PI / 180.0f32 rounds correctly without having to provide a more precise constant, although sadly it's still incorrect for FRAC_PI_6.)

[BubbaSheen]

What test is this refrenced too??

…

On Mar 5, 2018 1:39 PM, "varkor" ***@***.***> wrote: @cuviper <https://github.com/cuviper>: I can confirm that my tests are saying the same thing: between 0x00004000 and 0x00008000, there are 273 inputs with erroneous output for the current method and only 179 for the division method :) (It's also pleasant that f32::consts::PI / 180.0f32 rounds correctly without having to provide a more precise constant, although sadly it's still incorrect for FRAC_PI_6.) — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#48617 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AhXLZh-y6D-MrxaM9grXqFYAMLGDFRynks5tbYZngaJpZM4SXMRY> .

[varkor]

@BubbaSheen: I have an exhaustive test set for this issue, described here.

@cuviper: unfortunately, though the method seems to be more accurate on average, there are a significant number of inputs for which the error is more than 1 ULP (e.g. 0.000000000000000000000000000000000000000000003f32). Specifically, I suspect it's specifically the denormal numbers, though I haven't double-checked this yet.

[cuviper]

It's pretty compelling that the old way was correct only ~35% of the time, and the current way is ~86% correct. The division appears only slightly more accurate on average, ~91%, but it will be slower, and if it has worse ULP errors, that may not be an acceptable trade-off.

[hanna-kruppe]

The division [...] will be slower,

It can and will be constant folded.

That is not to say I am in favor of the division method. Being exactly right on a few more inputs (but not all) seems pretty useless, while having a stricter bound on the rounding error is pretty nice (at least when it takes you from 2 ULP to 1 ULP, as opposed to 300 ULP to 299 ULP).

[varkor]

Update: I was in a bit of a rush when I made that last comment, and hadn't spotted that I hadn't adjusted the code for the new test in the Mathematica code, which is why it was causing issues 😑

The division method looks correct for really small inputs too: from 0x0 to 0x8000, the current method makes 361 mistakes, whereas the division-based approach makes 243, which seems to match up with random testing. So perhaps the division is worth looking at...

[cuviper]

The division [...] will be slower,

It can and will be constant folded.

I'm sure the PI/180 will be, but what about this self / some_const when self is not constant? Will LLVM transform that to a multiplication, the same way it sometimes does for integer division?

[hanna-kruppe]

@cuviper Oh, duh, you're right. No, transforming a float division by a constant into a multiplication is not generally correct.

[varkor]

This seems to come down to accuracy vs speed again, then. I'm probably still inclined to stay away from the division, if only because the motivation all along has been to prioritise performance where possible. Really though, at this point, neither seems to have a significant drawback, and I think either could be reasonably chosen.

[nikomatsakis]

I'm trying to get up to speed here. It definitely seems to me like the consensus between @cuviper, @rkruppe, and @varkor is that we do not want to revert -- that the new method is far more accurate, even if some identities and particular cases no longer work.

(Point of clarification: @cuviper, @varkor -- when you say "the current method", you are referring to the behavior on nightly that landed in #47919?)

Do we know of code in the wild that is affected by this change? Do we expect that? (e.g., crater might be detecting tests that fail?)

[cuviper]

Yeah, I think what we have now is better in general.

(Point of clarification: @cuviper, @varkor -- when you say "the current method", you are referring to the behavior on nightly that landed in #47919?)

Yes.

Do we know of code in the wild that is affected by this change? Do we expect that? (e.g., crater might be detecting tests that fail?)

I discovered this in num-traits when I switched Float and FloatCore from implementing these methods on their own, to instead forward to the native methods in std. See here.

But this new error is only 1 ULP from the ideal, which makes me realize that the delta < 1e-6 test we had was already too narrow. At this magnitude, it was effectively demanding perfect accuracy, rather than allowing some floating point slop as intended.

[alexcrichton]

It sounds like the current plan of action here is to not revert, so should this be closed as "wontfix" effectively?

[cuviper]

@alexcrichton I think so -- I just didn't want to claim final authority. 😄

[alexcrichton]

Ok! In that case I'm going to close, but we can of course reopen if it keeps coming up!