Complete loss of precision for trig reduction in non-default rounding modes

[stephentyrone]

The trig-reduction implementations assume that the rounding mode is round-to-nearest, and do not function correctly if it isn't. Example: sinf(-0x1.4665d2p+27). The mathematically precise result is approximately -0x1.927bcc77af475p-25, but in round-toward-zero or round-down, openlibm computes a result of -0x1.a2c7eep-17.

This behavior is especially bad for the "medium argument" path in rem_pio2f. When the argument is so large that you punt to rem_pio2, double-precision has enough extra bits to still deliver relatively good results.

[simonbyrne]

I'm far from an expert on these things, having only really discovered this by a similar accident, but as I understand it rounding modes typically are only correct for basic arithmetic operation (and sqrt). There is CRlibm, which claims to be much more accurate for these cases, but even then at some point you still suffer from the table-maker's dilemma for any operation involving transcendental numbers.

This should be documented though.

[simonbyrne]

Or alternatively, we could save and restore the rounding mode, doing the actual computation in round-to-nearest, which I understand is what glibc does.

[StefanKarpinski]

I think that saving and restoring is probably the best approach. The main use case for changing rounding modes is verifying algorithm stability, which doesn't really work if our libm functions return garbage – testing their algorithm stability is not the point – they are one of the few things in the system that we actually know are really stable across the full range of values. The real trick is going back to default rounding in libm without a big performance hit.

[stephentyrone]

@simonbyrne I wouldn't expect openlibm to deliver correctly rounded results (we have CRlibm for that), but a high-quality math library should at the very least satisfy a small (in ULPs) error bound regardless of the prevailing rounding mode -- ideally less than an ULP.

@StefanKarpinski That's undoubtedly the simplest fix. In the long run, if someone gets ambitious they may want to investigate doing argument reduction using integer operations (on 64-bit platforms this can actually be quite efficient, and they are of course not effected by rounding mode).