[0/9] Improve hypot()

  This patchset uses a different algorithm for hypot()  based on the
'An Improved Algorithm for hypot(a,b)' by Carlos F. Borges [1].  This
method is also used by Julia language.  The motivation for this change
are:

  1. The new algorithm is more precise without large performance
     differences.
  2. It allows to consolidate the implementations, since it favor
     floating-point operations over integer ones (as done for powerpc).
     This might be a boost for some architectures as well.

The current hypot() implementation seems already to be bounded to a
maximum of 1 ulp of error, however the new proposed algorithm shows an
slight precision improvement by showing more correctly rounded results.

 With a random 1e9 inputs for different float format I see:

  - An improvement from 3427362 to 18457 results with 1 ulp of
    error for Binary64.
  - An improvement from 233442 to 1274 results with 1 ulp of
    error for Binary96 (x86_64).
  - An improvement from 453045 to 1294 results with 1 ulp of
    error for Binary96 (x86_64).

Along with the newer implementation, the last patch also optimizes the
call with similar work done for other math functions by removing the
POSIX error handling (which is a big performance boost in some cases
and architectures).

I have adapted the dbl-64, ldbl-96, and ldbl-128, the flt-32 is not
required since it calls the dbl-64 one.  I have not adapated ldbl-128ibm
since the format has a lot of caveats and IBM aims to move to ldbl-128.

[1] https://arxiv.org/pdf/1904.09481.pdf
[2] https://github.com/JuliaLang/julia/commit/4a046009a3362ab5e17d369641dbbc9657eb680c

Adhemerval Zanella (9):
  benchtests: Make hypot input random
  benchtests: Add hypotf
  math: Simplify hypotf implementation
  math: Use an improved algorithm for hypot (dbl-64)
  math: Use an improved algorithm for hypotl (ldbl-96)
  math: Use an improved algorithm for hypotl (ldbl-128)
  math: Remove powerpc e_hypot
  i386: Move hypot implementation to C
  math: Remove the error handling wrapper from hypot and hypotf

 benchtests/Makefile                           |    2 +-
 benchtests/hypot-inputs                       | 1014 ++++++++++++++++-
 benchtests/hypotf-inputs                      | 1006 ++++++++++++++++
 math/Versions                                 |    2 +
 math/w_hypot.c                                |    8 +
 math/w_hypot_compat.c                         |   13 +-
 math/w_hypotf.c                               |    8 +
 math/w_hypotf_compat.c                        |    6 +-
 sysdeps/i386/fpu/e_hypot.S                    |   75 --
 sysdeps/i386/fpu/e_hypot.c                    |   54 +
 sysdeps/i386/fpu/e_hypotf.S                   |   64 --
 sysdeps/ieee754/dbl-64/e_hypot.c              |  247 ++--
 sysdeps/ieee754/dbl-64/w_hypot.c              |    1 +
 sysdeps/ieee754/flt-32/e_hypotf.c             |   72 +-
 sysdeps/ieee754/flt-32/w_hypotf.c             |    1 +
 sysdeps/ieee754/ldbl-128/e_hypotl.c           |  222 ++--
 sysdeps/ieee754/ldbl-96/e_hypotl.c            |  227 ++--
 sysdeps/mach/hurd/i386/libm.abilist           |    2 +
 sysdeps/mach/hurd/libhurduser.abilist         |    0
 sysdeps/mach/libmachuser.abilist              |    0
 sysdeps/powerpc/fpu/e_hypot.c                 |   87 --
 sysdeps/powerpc/fpu/e_hypotf.c                |   78 --
 .../powerpc32/power4/fpu/multiarch/Makefile   |    5 +-
 .../power4/fpu/multiarch/e_hypot-power7.c     |   23 -
 .../power4/fpu/multiarch/e_hypot-ppc32.c      |   23 -
 .../powerpc32/power4/fpu/multiarch/e_hypot.c  |   33 -
 .../power4/fpu/multiarch/e_hypotf-power7.c    |   23 -
 .../power4/fpu/multiarch/e_hypotf-ppc32.c     |   23 -
 .../powerpc32/power4/fpu/multiarch/e_hypotf.c |   33 -
 sysdeps/unix/sysv/linux/aarch64/libm.abilist  |    2 +
 sysdeps/unix/sysv/linux/alpha/libm.abilist    |    2 +
 sysdeps/unix/sysv/linux/arm/be/libm.abilist   |    2 +
 sysdeps/unix/sysv/linux/arm/le/libm.abilist   |    2 +
 sysdeps/unix/sysv/linux/hppa/libm.abilist     |    2 +
 sysdeps/unix/sysv/linux/i386/libm.abilist     |    2 +
 .../sysv/linux/m68k/coldfire/libm.abilist     |    2 +
 .../unix/sysv/linux/m68k/m680x0/libm.abilist  |    2 +
 .../sysv/linux/microblaze/be/libm.abilist     |    2 +
 .../sysv/linux/microblaze/le/libm.abilist     |    2 +
 .../unix/sysv/linux/mips/mips32/libm.abilist  |    2 +
 .../unix/sysv/linux/mips/mips64/libm.abilist  |    2 +
 sysdeps/unix/sysv/linux/nios2/libm.abilist    |    2 +
 .../linux/powerpc/powerpc32/fpu/libm.abilist  |    2 +
 .../powerpc/powerpc32/nofpu/libm.abilist      |    2 +
 .../linux/powerpc/powerpc64/be/libm.abilist   |    2 +
 .../linux/powerpc/powerpc64/le/libm.abilist   |    2 +
 .../unix/sysv/linux/s390/s390-32/libm.abilist |    2 +
 .../unix/sysv/linux/s390/s390-64/libm.abilist |    2 +
 sysdeps/unix/sysv/linux/sh/be/libm.abilist    |    2 +
 sysdeps/unix/sysv/linux/sh/le/libm.abilist    |    2 +
 .../sysv/linux/sparc/sparc32/libm.abilist     |    2 +
 .../sysv/linux/sparc/sparc64/libm.abilist     |    2 +
 .../unix/sysv/linux/x86_64/64/libm.abilist    |    2 +
 .../unix/sysv/linux/x86_64/x32/libm.abilist   |    2 +
 54 files changed, 2483 insertions(+), 919 deletions(-)
 create mode 100644 benchtests/hypotf-inputs
 create mode 100644 math/w_hypot.c
 create mode 100644 math/w_hypotf.c
 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
 create mode 100644 sysdeps/ieee754/dbl-64/w_hypot.c
 create mode 100644 sysdeps/ieee754/flt-32/w_hypotf.c
 delete mode 100644 sysdeps/mach/hurd/libhurduser.abilist
 delete mode 100644 sysdeps/mach/libmachuser.abilist
 delete mode 100644 sysdeps/powerpc/fpu/e_hypot.c
 delete mode 100644 sysdeps/powerpc/fpu/e_hypotf.c
 delete mode 100644 sysdeps/powerpc/powerpc32/power4/fpu/multiarch/e_hypot-power7.c
 delete mode 100644 sysdeps/powerpc/powerpc32/power4/fpu/multiarch/e_hypot-ppc32.c
 delete mode 100644 sysdeps/powerpc/powerpc32/power4/fpu/multiarch/e_hypot.c
 delete mode 100644 sysdeps/powerpc/powerpc32/power4/fpu/multiarch/e_hypotf-power7.c
 delete mode 100644 sysdeps/powerpc/powerpc32/power4/fpu/multiarch/e_hypotf-ppc32.c
 delete mode 100644 sysdeps/powerpc/powerpc32/power4/fpu/multiarch/e_hypotf.c
  

>  With a random 1e9 inputs for different float format I see:
> 
>   - An improvement from 3427362 to 18457 results with 1 ulp of
>     error for Binary64.
>   - An improvement from 233442 to 1274 results with 1 ulp of
>     error for Binary96 (x86_64).
>   - An improvement from 453045 to 1294 results with 1 ulp of
>     error for Binary96 (x86_64).

I guess the last one should read Binary128.

Another comparison is to check for the maximal known error.
For glibc-2.34 the maximal known errors are as follows (in ulps, with
corresponding inputs), determined with [1]:

binary32  0.500 0x1.3ac98p+67,-0x1.ba5ec2p+77
binary64  0.987 -0x0.5a934b7eac967p-1022,-0x0.b5265a7e06b82p-1022
binary96  0.981 0x1.73f339f61eda21dp-16384l,0x2.e45f9f9500877e2p-16384l
binary128 0.985 -0x2.d8311789103b76133ea1d5bc38c4p-16384,-0x1.6d85492006d7dcc6cc52938684p-16384

With the new code I get (with threshold 1000000, except for binary128 where I
used 200000):

binary32  0.500 0x1.3ac98p+67,-0x1.ba5ec2p+77 [same]
binary64  0.792 0x0.603e52daf0bfdp-1022,-0x0.a622d0a9a433bp-1022
binary96  0.584 -0x2.97b86706043d619p+7240l,0x1.8256bdd12d2e163ep+7240l
binary128 0.749 0x2.2d5faf4036d6e68566f01054612p-8192,0x3.5738e8e2505f5d1fc2973716f05p-8192

Paul

[1] https://gitlab.inria.fr/zimmerma/math_accuracy/-/blob/master/binary64/check_sample2.c
  

On 07/10/2021 08:03, Paul Zimmermann wrote:
>>  With a random 1e9 inputs for different float format I see:
>>
>>   - An improvement from 3427362 to 18457 results with 1 ulp of
>>     error for Binary64.
>>   - An improvement from 233442 to 1274 results with 1 ulp of
>>     error for Binary96 (x86_64).
>>   - An improvement from 453045 to 1294 results with 1 ulp of
>>     error for Binary96 (x86_64).
> 
> I guess the last one should read Binary128.

Yeap.

> 
> Another comparison is to check for the maximal known error.
> For glibc-2.34 the maximal known errors are as follows (in ulps, with
> corresponding inputs), determined with [1]:
> 
> binary32  0.500 0x1.3ac98p+67,-0x1.ba5ec2p+77
> binary64  0.987 -0x0.5a934b7eac967p-1022,-0x0.b5265a7e06b82p-1022
> binary96  0.981 0x1.73f339f61eda21dp-16384l,0x2.e45f9f9500877e2p-16384l
> binary128 0.985 -0x2.d8311789103b76133ea1d5bc38c4p-16384,-0x1.6d85492006d7dcc6cc52938684p-16384
> 
> With the new code I get (with threshold 1000000, except for binary128 where I
> used 200000):
> 
> binary32  0.500 0x1.3ac98p+67,-0x1.ba5ec2p+77 [same]
> binary64  0.792 0x0.603e52daf0bfdp-1022,-0x0.a622d0a9a433bp-1022
> binary96  0.584 -0x2.97b86706043d619p+7240l,0x1.8256bdd12d2e163ep+7240l
> binary128 0.749 0x2.2d5faf4036d6e68566f01054612p-8192,0x3.5738e8e2505f5d1fc2973716f05p-8192

Thanks to double check that precision does improve in all cases.

> 
> Paul
> 
> [1] https://gitlab.inria.fr/zimmerma/math_accuracy/-/blob/master/binary64/check_sample2.c
>