[8/9] i386: Move hypot implementation to C

  The generic hypotf is slight slower, mostly due the tricks the assembly
does to optimize the isinf/isnan/issignaling.  Results on a Ryzen 5900X
with gcc 10.3.1:

master:
  "hypotf": {
   "workload-random": {
    "duration": 3.76493e+09,
    "iterations": 6.8e+07,
    "reciprocal-throughput": 38.4243,
    "latency": 72.309,
    "max-throughput": 2.60252e+07,
    "min-throughput": 1.38295e+07
   }
  }

patched:
  "hypotf": {
   "workload-random": {
    "duration": 3.78098e+09,
    "iterations": 6.8e+07,
    "reciprocal-throughput": 28.9206,
    "latency": 82.2848,
    "max-throughput": 3.45774e+07,
    "min-throughput": 1.21529e+07
   }
  }

The generic hypot is way slower, since the optimized implementation
uses the i386 default excessive precision to issue the operation
directly.  A similar implementation is provided instead of using
the generic implementation:

master:
  "hypot": {
   "workload-random": {
    "duration": 3.7452e+09,
    "iterations": 6.6e+07,
    "reciprocal-throughput": 40.3203,
    "latency": 73.1707,
    "max-throughput": 2.48014e+07,
    "min-throughput": 1.36667e+07
   }
  }

patched:
./benchtests/bench-hypot
  "hypot": {
   "workload-random": {
    "duration": 3.72606e+09,
    "iterations": 7.6e+07,
    "reciprocal-throughput": 25.9437,
    "latency": 72.1104,
    "max-throughput": 3.8545e+07,
    "min-throughput": 1.38676e+07
   }
  }

Checked on i686-linux-gnu.
---
 sysdeps/i386/fpu/e_hypot.S  | 75 -------------------------------------
 sysdeps/i386/fpu/e_hypot.c  | 42 +++++++++++++++++++++
 sysdeps/i386/fpu/e_hypotf.S | 64 -------------------------------
 3 files changed, 42 insertions(+), 139 deletions(-)
 delete mode 100644 sysdeps/i386/fpu/e_hypot.S
 create mode 100644 sysdeps/i386/fpu/e_hypot.c
 delete mode 100644 sysdeps/i386/fpu/e_hypotf.S
  

On Wed, 6 Oct 2021, Adhemerval Zanella via Libc-alpha wrote:

> +/* The i386 allows ot use the default excess of precision to optimize the

s/ot/to/

> +   hypot implementation, since internal multiplication and sqrt is carried
> +   with 80-bit FP type.  */
> +double
> +__ieee754_hypot (double x, double y)
> +{
> +  if ((isinf (x) || isinf (y))
> +      && !issignaling (x) && !issignaling (y))
> +    return INFINITY;
> +  if (isnan (x) || isnan (y))
> +    return x + y;
> +
> +  double r = math_narrow_eval (sqrt (x * x + y * y));

There is no guarantee of when excess precision might or might not be used 
for intermediate computations, or whether it might not be used (so 
resulting in spurious overflows).  I'd expect explicit conversions of x 
and y to long double here, use of sqrtl and conversion of the result to 
double.
  

On 06/10/2021 15:37, Joseph Myers wrote:
> On Wed, 6 Oct 2021, Adhemerval Zanella via Libc-alpha wrote:
> 
>> +/* The i386 allows ot use the default excess of precision to optimize the
> 
> s/ot/to/

Ack.

> 
>> +   hypot implementation, since internal multiplication and sqrt is carried
>> +   with 80-bit FP type.  */
>> +double
>> +__ieee754_hypot (double x, double y)
>> +{
>> +  if ((isinf (x) || isinf (y))
>> +      && !issignaling (x) && !issignaling (y))
>> +    return INFINITY;
>> +  if (isnan (x) || isnan (y))
>> +    return x + y;
>> +
>> +  double r = math_narrow_eval (sqrt (x * x + y * y));
> 
> There is no guarantee of when excess precision might or might not be used 
> for intermediate computations, or whether it might not be used (so 
> resulting in spurious overflows).  I'd expect explicit conversions of x 
> and y to long double here, use of sqrtl and conversion of the result to 
> double.
> 

Something like the below?

  long double lx = x;
  long double ly = y;
  double r = sqrtl (lx * lx + ly * ly);
  return r;
  

On 06/10/2021 16:19, Adhemerval Zanella wrote:
> 
> 
> On 06/10/2021 15:37, Joseph Myers wrote:
>> On Wed, 6 Oct 2021, Adhemerval Zanella via Libc-alpha wrote:
>>
>>> +/* The i386 allows ot use the default excess of precision to optimize the
>>
>> s/ot/to/
> 
> Ack.
> 
>>
>>> +   hypot implementation, since internal multiplication and sqrt is carried
>>> +   with 80-bit FP type.  */
>>> +double
>>> +__ieee754_hypot (double x, double y)
>>> +{
>>> +  if ((isinf (x) || isinf (y))
>>> +      && !issignaling (x) && !issignaling (y))
>>> +    return INFINITY;
>>> +  if (isnan (x) || isnan (y))
>>> +    return x + y;
>>> +
>>> +  double r = math_narrow_eval (sqrt (x * x + y * y));
>>
>> There is no guarantee of when excess precision might or might not be used 
>> for intermediate computations, or whether it might not be used (so 
>> resulting in spurious overflows).  I'd expect explicit conversions of x 
>> and y to long double here, use of sqrtl and conversion of the result to 
>> double.
>>
> 
> Something like the below?
> 
>   long double lx = x;
>   long double ly = y;
>   double r = sqrtl (lx * lx + ly * ly);
>   return r;
> 

In fact we need to handle underflow:

  long double lx = x;
  long double ly = y;
  double r = sqrtl (lx * lx + ly * ly);
  math_check_force_underflow_nonneg (r);
  return r;
  

On Wed, 6 Oct 2021, Adhemerval Zanella via Libc-alpha wrote:

> In fact we need to handle underflow:
> 
>   long double lx = x;
>   long double ly = y;
>   double r = sqrtl (lx * lx + ly * ly);

You probably also need a math_narrow_eval here to ensure the value checked 
for overflow and returned has definitely been properly narrowed to double 
first.

>   math_check_force_underflow_nonneg (r);
>   return r;
  

> > In fact we need to handle underflow:
> > 
> >   long double lx = x;
> >   long double ly = y;
> >   double r = sqrtl (lx * lx + ly * ly);
> 
> You probably also need a math_narrow_eval here to ensure the value checked 
> for overflow and returned has definitely been properly narrowed to double 
> first.

right, If x = y = 0x1.fffffffffffffp1023, then lx * lx + ly * ly will not
overflow as long double, nor sqrtl(.), but the conversion to double will.
This should probably be added as a test case (if not already).

Paul
  

On Thu, 7 Oct 2021, Paul Zimmermann wrote:

> > > In fact we need to handle underflow:
> > > 
> > >   long double lx = x;
> > >   long double ly = y;
> > >   double r = sqrtl (lx * lx + ly * ly);
> > 
> > You probably also need a math_narrow_eval here to ensure the value checked 
> > for overflow and returned has definitely been properly narrowed to double 
> > first.
> 
> right, If x = y = 0x1.fffffffffffffp1023, then lx * lx + ly * ly will not
> overflow as long double, nor sqrtl(.), but the conversion to double will.
> This should probably be added as a test case (if not already).

This is the sort of test that could be represented simply by a line

hypot max max

in auto-libm-test-in (which would then include all pairs of the maximum 
values of supported floating-point formats), but is problematic to add in 
practice because of limitations of the libm test machinery.

For example, that line would generate tests including hypot (DBL_MAX, 
FLT_MAX) for double.  Logically, by the rules glibc follows about 
acceptable errors (that most libm functions should act like they compute 
an infinite-precision result close to the mathematical value, then round 
it to the relevant type):

* If, for round-to-nearest, this overflows (generates Inf as the result 
with FE_OVERFLOW raised), that's a small (1ulp) error and so OK.

* If, for round-upward, this does not overflow (generates DBL_MAX as the 
result without FE_OVERFLOW raised), again that's a 1ulp error and so OK.

But the libm-test machinery isn't able to represent the information that a 
small error might result in overflowing when a correctly rounded result 
does not, or vice versa (a similar issue can also arise for underflow).  
And if the mathematical result in an overflow case is more than a few ulp 
above the overflow threshold, such a variation becomes a large error which 
is no longer OK.

We have a few cases close to the overflow threshold that are XFAILed in 
auto-libm-test-in (for lgamma, for example).  But XFAILing like that is 
problematic for "hypot max max", because that line generates lots of 
different tests, and only a few of them are appropriate to XFAIL.