[PATCH 3/4] Use libc_fe* macros in ldbl-128/e_expl.c.

[Stefan Liebler]

Unfortunately, this patch is responsible for testfails on x86_64:

math/test-float128-exp.out:
Failure: exp (-0x1p-10000): Exception "Underflow" set
Failure: exp (-0x2p-16384): Exception "Underflow" set
...

math/test-float128-cexp.out:
Failure: Real part of: cexp (0x2p-16384 - 0x4p-1076 i): Exception 
"Underflow" set
Failure: Real part of: cexp (0x2p-16384 - 0x8p-152 i): Exception 
"Underflow" set
...

I've stepped through "expf128 (0x1p-10000)" in 
sysdeps/ieee754/float128/../ldbl-128/e_expl.c:
151: libc_feholdexcept_setroundl (&oldenv, FE_TONEAREST);    // before 
this patch: feholdexcept (&oldenv); fesetround (FE_TONEAREST);
199: x22 = x + x*x*(P1+x*(P2+x*(P3+x*(P4+x*(P5+x*P6)))));
203: libc_fesetenvl (&oldenv);    // before this patch: fesetenv (&oldenv);

During the evaluation of x22 the underflow exception occures while:
<__ieee754_expf128+920>  callq  0x7ffff7f2a7c0 <__multf3>
which calls __sfp_handle_exceptions():
ae8be:       40 f6 c7 10             test   $0x10,%dil
ae8c2:       74 0f                   je     ae8d3 
<__sfp_handle_exceptions+0x73>
ae8c4:       d9 74 24 d8             fnstenv -0x28(%rsp)
ae8c8:       66 83 4c 24 dc 10       orw    $0x10,-0x24(%rsp)
ae8ce:       d9 64 24 d8             fldenv -0x28(%rsp)
ae8d2:       9b                      fwait

According to sysdeps/x86/fpu/fenv_private.h:
#ifdef __x86_64__
# define libc_feholdexcept_setroundf128 libc_feholdexcept_setround_sse
# define libc_fesetenvf128	libc_fesetenv_sse
#else
# define libc_feholdexcept_setroundf128	default_libc_feholdexcept_setround
# define libc_fesetenvf128	default_libc_fesetenv
#endif

// On my machine:
# define STMXCSR "stmxcsr"
# define LDMXCSR "ldmxcsr"

static __always_inline void
libc_feholdexcept_setround_sse (fenv_t *e, int r)
{
   unsigned int mxcsr;
   asm (STMXCSR " %0" : "=m" (*&mxcsr));
   e->__mxcsr = mxcsr;
   mxcsr = ((mxcsr | 0x1f80) & ~0x603f) | (r << 3);
   asm volatile (LDMXCSR " %0" : : "m" (*&mxcsr));
}

Whereas the feholdexcept() function is using in
./sysdeps/x86_64/fpu/feholdexcpt.c:
int
__feholdexcept (fenv_t *envp)
{
   unsigned int mxcsr;

   /* Store the environment.  Recall that fnstenv has a side effect of
      masking all exceptions.  Then clear all exceptions.  */
   __asm__ ("fnstenv %0\n\t"
	   "stmxcsr %1\n\t"
	   "fnclex"
	   : "=m" (*envp), "=m" (envp->__mxcsr));

   /* Set the SSE MXCSR register.  */
   mxcsr = (envp->__mxcsr | 0x1f80) & ~0x3f;
   __asm__ ("ldmxcsr %0" : : "m" (*&mxcsr));

   return 0;
}

I assume that the underflow exception keeps active as the pair of 
fnstenv / fldenv is missing if libc_feholdexcept_setroundf128 / 
libc_fesetenvf128 is used instead of feholdexcept, fesetround and fesetenv.
As I'm not familiar with float128 on x86_64, can anybody please help?

Bye,
Stefan

On 3/25/20 11:06 AM, Stefan Liebler wrote:
> The calls to feholdexcept, fesetround and fesetenv are replaced
> by the libc_fe* macros.
> ---
>   sysdeps/ieee754/ldbl-128/e_expl.c | 8 +++++---
>   1 file changed, 5 insertions(+), 3 deletions(-)
> 
> diff --git a/sysdeps/ieee754/ldbl-128/e_expl.c b/sysdeps/ieee754/ldbl-128/e_expl.c
> index 37c1538c08..104ace1690 100644
> --- a/sysdeps/ieee754/ldbl-128/e_expl.c
> +++ b/sysdeps/ieee754/ldbl-128/e_expl.c
> @@ -66,6 +66,7 @@
>   #include <inttypes.h>
>   #include <math-barriers.h>
>   #include <math_private.h>
> +#include <fenv_private.h>
>   #include <math-underflow.h>
>   #include <stdlib.h>
>   #include "t_expl.h"
> @@ -146,9 +147,10 @@ __ieee754_expl (_Float128 x)
>         union ieee854_long_double ex2_u, scale_u;
>         fenv_t oldenv;
>   
> -      feholdexcept (&oldenv);
>   #ifdef FE_TONEAREST
> -      fesetround (FE_TONEAREST);
> +      libc_feholdexcept_setroundl (&oldenv, FE_TONEAREST);
> +#else
> +      libc_feholdexceptl (&oldenv);
>   #endif
>   
>         /* Calculate n.  */
> @@ -198,7 +200,7 @@ __ieee754_expl (_Float128 x)
>         math_force_eval (x22);
>   
>         /* Return result.  */
> -      fesetenv (&oldenv);
> +      libc_fesetenvl (&oldenv);
>   
>         result = x22 * ex2_u.d + ex2_u.d;
>   
>