[Bug math/23253] tgamma test suite failures on i686 with -march=x86-64 -mtune=generic -mfpmath=sse

["joseph at codesourcery dot com" <sourceware-bugzilla at sourceware dot org>]

https://sourceware.org/bugzilla/show_bug.cgi?id=23253

--- Comment #5 from joseph at codesourcery dot com <joseph at codesourcery dot com> ---
On Fri, 8 Jun 2018, fweimer at redhat dot com wrote:

> --- Comment #4 from Florian Weimer <fweimer at redhat dot com> ---
> (In reply to joseph@codesourcery.com from comment #1)
> > On Fri, 1 Jun 2018, fweimer at redhat dot com wrote:
> > 
> > > Failure: tgamma_upward (-0xb.4ffffp+4): errno set to 0, expected 34 (ERANGE)
> > 
> > I advise checking what __ieee754_gamma_r returns (presumably not the 
> > expected zero),
> 
> GDB says it's 4.9406564584124654e-324, assuming that I got the right invocation
> in the test.  It's from the final invocation right before the Failure: printf.

The mathematical tgamma value here is (as per auto-libm-test-out-tgamma) 
about -0x5.fa2b7cb1342f0619b84e4c550dp-1088 (which thus rounds upward to 
zero for double).  The peculiar __ieee754_gamma_r interface (separate 
signgam only really makes sense for lgamma not tgamma...) thus means that 
what's meant to happen is: the initial computation of ret (the absolute 
value of the result), in round-to-nearest mode, produces 0 after 
math_narrow_eval.  The "else if (ret == 0)" logic should then come into 
play, compute math_narrow_eval (-__copysign (DBL_MIN, ret) * DBL_MIN) in 
the original rounding mode (which in this case means rounding upward to 
-0), and then negate it to +0.

> 100           double ret = (__ieee754_pow (x_adj_mant, x_adj)
> 101                         * __ieee754_exp2 (x_adj_log2 * x_adj_frac)
> 102                         * __ieee754_exp (-x_adj)
> 103                         * sqrt (2 * M_PI / x_adj)
> 104                         / prod);
> 
> And the hand-written assembler implementations use the FPU, not SSE2.  So in
> order to call them, the FPU would have to be switched to round-to-nearest as
> well.  Or we drop those assembler implementations.

For those implementations, I don't think the rounding mode should 
particularly matter; in the case where underflow isn't just using DBL_MIN 
* DBL_MIN (in the range above -184 where whether the result underflows 
depends on the fractional part of the argument), the underflow should 
actually come from the "ret = __scalbn (tret, -exp2_adj);".  So the issue 
is probably the i386 s_scalbn.S.

So, yes, in this particular case either the x87 rounding mode needs to be 
set to round-to-nearest, or you need to use a different scalbn 
implementation for building with -mfpmath=sse (but it can't just be using 
a generic one as-is because of the i386-specific symbol versioning for 
scalbln in libc).  There may well be a case that if building with 
-mfpmath=sse you want to avoid most of the float/double implementations 
using x87 assembly (subject to symbol versioning issues, of course), which 
I suppose would need different sysdeps directories depending on the 
-mfpmath= choice.

If you want just to set the x87 rounding mode in this particular case you 
need to make fenv_private.h do something different here, hopefully without 
also needing to set both rounding modes in lots of other places that only 
need one set.  Maybe a macro fenv_private.h tests and an i386-specific 
wrapper round e_gamma_r.c that defines that macro, for example?

-- 
You are receiving this mail because:
You are on the CC list for the bug.