[PATCH] Ensure calculations happen with desired rounding mode in y1lf128

[PATCH] Ensure calculations happen with desired rounding mode in y1lf128

[Wilco Dijkstra]

Hi Michael,

> I don't know if this patch should be committed as is, although it does
> fix an observed failure for me.

The patch looks good to me, but this is a more general issue...

> My read of the discussion in the gcc bug I filed about this
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=106574 is that every single
> use of a SET_RESTORE_ROUND macro is vulnerable to this. I can't think of
> a generic way to fix this -- I guess you could have a macro to call an
> always_inline function with a rounding mode set? But that would uglify
> the control flow in the code quite a bit.

Indeed, there is no general way to block optimizations across a boundary in
GCC if they do not involve memory accesses. As a result it is never safe to
change the rounding mode inside a function...

All math functions using the SET_RESTORE_ROUND macros will need similar
barriers. Note that it is feasible to remove these macros altogether and fix
any issues (a slightly larger ULP is acceptable for non-nearest rounding).
Given rounding mode changes are generally expensive, this also improves
performance (though that may not be important for 128-bit floats).

Cheers,
Wilco

Re: [PATCH] Ensure calculations happen with desired rounding mode in y1lf128

[Joseph Myers]

On Fri, 12 Aug 2022, Wilco Dijkstra via Libc-alpha wrote:

> All math functions using the SET_RESTORE_ROUND macros will need similar
> barriers. Note that it is feasible to remove these macros altogether and fix
> any issues (a slightly larger ULP is acceptable for non-nearest rounding).
> Given rounding mode changes are generally expensive, this also improves
> performance (though that may not be important for 128-bit floats).

This one is a case where SET_RESTORE_ROUND is used to reduce error 
accumulation to keep the errors within the bounds accepted by the 
testsuite (see bug 16824).  In such cases, it may indeed be possible to 
change the algorithm to one that has less total error accumulation 
possible in any rounding mode so the results are sufficiently accurate 
independent of rounding mode without needing SET_RESTORE_ROUND.

In other cases, the manipulation of the floating-point environment is 
needed for correctness, e.g. to avoid spurious exceptions or to implement 
round-to-odd for functions that must be correctly rounding, or it's 
because algorithms for higher internal precision are used (Dekker / Knuth) 
that are only correct in round-to-nearest most and much larger errors 
might occur if those are used in the wrong rounding mode.

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: [PATCH] Ensure calculations happen with desired rounding mode in y1lf128

[Michael Hudson-Doyle]

On Tue, 16 Aug 2022 at 09:00, Joseph Myers <joseph@codesourcery.com> wrote:

> On Fri, 12 Aug 2022, Wilco Dijkstra via Libc-alpha wrote:
>
> > All math functions using the SET_RESTORE_ROUND macros will need similar
> > barriers.


Sounds fun. Really, wrapping mode changes around a function call sounds
saner but also tedious in a different way.


> > Note that it is feasible to remove these macros altogether and fix
> > any issues (a slightly larger ULP is acceptable for non-nearest
> rounding).
>

Removing the rounding mode change entirely makes the errors quite a lot
worse, as the bug Joseph links to shows.


> > Given rounding mode changes are generally expensive, this also improves
> > performance (though that may not be important for 128-bit floats).


> This one is a case where SET_RESTORE_ROUND is used to reduce error
> accumulation to keep the errors within the bounds accepted by the
> testsuite (see bug 16824).  In such cases, it may indeed be possible to
> change the algorithm to one that has less total error accumulation
> possible in any rounding mode so the results are sufficiently accurate
> independent of rounding mode without needing SET_RESTORE_ROUND.
>
> In other cases, the manipulation of the floating-point environment is
> needed for correctness, e.g. to avoid spurious exceptions or to implement
> round-to-odd for functions that must be correctly rounding, or it's
> because algorithms for higher internal precision are used (Dekker / Knuth)
> that are only correct in round-to-nearest most and much larger errors
> might occur if those are used in the wrong rounding mode.
>

It does seem likely that an algorithm that does not require setting a
rounding mode would in general be better than one that does but also that
this is not very realistic. But I guess my point is that SET_RESTORE_ROUND
without barriers is a footgun. I guess I should commit my patch and perhaps
see about writing some more for other uses of the macro?

Cheers,
mwh

Re: [PATCH] Ensure calculations happen with desired rounding mode in y1lf128

[Wilco Dijkstra]

Hi Michael,

> It does seem likely that an algorithm that does not require setting a rounding mode
> would in general be better than one that does but also that this is not very realistic.

Rewriting this function would not be trivial, but the math functions that are already
rewritten (like exp, pow, log) prove you can get fast and accurate results without
ever needing rounding mode changes.

When you design a polynomial to be very accurate, it actually works in all rounding
modes without extra effort. For this function it seems like whoever wrote it didn't
understand how to accurately evaluate polynomials - the input range of neval/deval
is (0.0, 4.0] so the repeated multiplies actually multiply the rounding errors...

> But I guess my point is that SET_RESTORE_ROUND without barriers is a footgun. 
> I guess I should commit my patch and perhaps see about writing some more for other
> uses of the macro?

The patch LGTM. Yes more patches would be welcome. It's a good idea to check
what happens if you remove the SET_RESTORE_ROUND - there will be cases where
the ULP is good enough.

Cheers,
Wilco

Re: [PATCH] Ensure calculations happen with desired rounding mode in y1lf128

[Joseph Myers]

On Wed, 17 Aug 2022, Wilco Dijkstra via Libc-alpha wrote:

> When you design a polynomial to be very accurate, it actually works in 
> all rounding modes without extra effort. For this function it seems like 

In the Bessel function case, accurate results near zeros in the 
intermediate range (the range where multiple terms of a series are needed 
before applying range reduction, see Harrison's paper 
<https://www.cl.cam.ac.uk/~jrh13/papers/bessel.pdf>) may require extra 
internal precision when computing the phase.  And some precision extension 
techniques only work in round-to-nearest (we only have that Bessel 
function approach implemented in glibc for float, where using double 
internally avoids such rounding mode problems for precision extension).

-- 
Joseph S. Myers
joseph@codesourcery.com

Re: [PATCH] Ensure calculations happen with desired rounding mode in y1lf128

[Paul Zimmermann]

> Rewriting this function would not be trivial, but the math functions that are already
> rewritten (like exp, pow, log) prove you can get fast and accurate results without
> ever needing rounding mode changes.

indeed, you can even get correct rounding (for all rounding modes) faster than
the current glibc implementations. See https://hal.inria.fr/hal-03721525
for example, or the current LLVM libc implementations (cf bottom table of
https://core-math.gitlabpages.inria.fr/).

Paul Zimmermann

[PATCH] Ensure calculations happen with desired rounding mode in y1lf128

[Michael Hudson-Doyle]

math/test-float128-y1 fails on x86_64 and ppc64el with gcc 12 and -O3,
because code inside a block guarded by SET_RESTORE_ROUNDL is being moved
after the rounding mode has been restored. Use math_force_eval to
prevent this (and insert some math_opt_barrier calls to prevent code
from being moved before the rounding mode is set).

Fixes #29463
---
I don't know if this patch should be committed as is, although it does
fix an observed failure for me.

My read of the discussion in the gcc bug I filed about this
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=106574 is that every single
use of a SET_RESTORE_ROUND macro is vulnerable to this. I can't think of
a generic way to fix this -- I guess you could have a macro to call an
always_inline function with a rounding mode set? But that would uglify
the control flow in the code quite a bit.
---
 sysdeps/ieee754/ldbl-128/e_j1l.c | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/sysdeps/ieee754/ldbl-128/e_j1l.c b/sysdeps/ieee754/ldbl-128/e_j1l.c
index 54c457681a..9a9c5c6f00 100644
--- a/sysdeps/ieee754/ldbl-128/e_j1l.c
+++ b/sysdeps/ieee754/ldbl-128/e_j1l.c
@@ -869,10 +869,13 @@ __ieee754_y1l (_Float128 x)
     {
       /* 0 <= x <= 2 */
       SET_RESTORE_ROUNDL (FE_TONEAREST);
+      xx = math_opt_barrier (xx);
+      x = math_opt_barrier (x);
       z = xx * xx;
       p = xx * neval (z, Y0_2N, NY0_2N) / deval (z, Y0_2D, NY0_2D);
       p = -TWOOPI / xx + p;
       p = TWOOPI * __ieee754_logl (x) * __ieee754_j1l (x) + p;
+      math_force_eval (p);
       return p;
     }
 
-- 
2.34.1