math: Use an improved algorithm for hypot (dbl-64)
authorWilco Dijkstra <Wilco.Dijkstra@arm.com>
Mon, 8 Mar 2021 20:07:39 +0000 (17:07 -0300)
committerAdhemerval Zanella <adhemerval.zanella@linaro.org>
Mon, 13 Dec 2021 12:02:34 +0000 (09:02 -0300)
commit6c848d70383e1dbe932ef41723ac0abfdeec7ca8
tree366dd2d23fb72a4a31b7e44a81d355f46627355a
parent7fe0ace3e289c88cab5014cef94e946fd695221f
math: Use an improved algorithm for hypot (dbl-64)

This implementation is based on the 'An Improved Algorithm for
hypot(a,b)' by Carlos F. Borges [1] using the MyHypot3 with the
following changes:

 - Handle qNaN and sNaN.
 - Tune the 'widely varying operands' to avoid spurious underflow
   due the multiplication and fix the return value for upwards
   rounding mode.
 - Handle required underflow exception for denormal results.

The main advantage of the new algorithm is its precision: with a
random 1e9 input pairs in the range of [DBL_MIN, DBL_MAX], glibc
current implementation shows around 0.34% results with an error of
1 ulp (3424869 results) while the new implementation only shows
0.002% of total (18851).

The performance result are also only slight worse than current
implementation.  On x86_64 (Ryzen 5900X) with gcc 12:

Before:

  "hypot": {
   "workload-random": {
    "duration": 3.73319e+09,
    "iterations": 1.12e+08,
    "reciprocal-throughput": 22.8737,
    "latency": 43.7904,
    "max-throughput": 4.37184e+07,
    "min-throughput": 2.28361e+07
   }
  }

After:

  "hypot": {
   "workload-random": {
    "duration": 3.7597e+09,
    "iterations": 9.8e+07,
    "reciprocal-throughput": 23.7547,
    "latency": 52.9739,
    "max-throughput": 4.2097e+07,
    "min-throughput": 1.88772e+07
   }
  }

Co-Authored-By: Adhemerval Zanella <adhemerval.zanella@linaro.org>
Checked on x86_64-linux-gnu and aarch64-linux-gnu.

[1] https://arxiv.org/pdf/1904.09481.pdf
sysdeps/ieee754/dbl-64/e_hypot.c