[libc-commits] [PATCH] D118157: [libc] Improve hypotf performance with different algorithm correctly rounded to all rounding modes.

[Tue Ly via Phabricator via libc-commits]

lntue marked an inline comment as not done.
lntue added inline comments.


================
Comment at: libc/src/math/generic/hypotf.cpp:29-30
+  // Compute the rounding error with Dekker's algorithm.
+  double err = ((sumSq - xSq) - ySq) + ((sumSq - ySq) - xSq);
+
+  // Take sqrt in double precision.
----------------
vinc17 wrote:
> zimmermann6 wrote:
> > I hadn't seen that trick to compute the rounding error, do you have a reference?
> > By the way, I'm not sure the reference to Dekker is appropriate. For me, Dekker's algorithm splits
> > two floating-point numbers in two each, and computes their product (high + low part) using 4 multiplies.
> If I understand correctly, `err` should get the rounding error of the sum. The algorithm is known as TwoSum. It needs 6 operations, including the sum `sumSq`, and this is the same number of operations as you have. But with 6 add/sub operations, I proved that there is only one algorithm (up to obvious symmetries) that works. And the above one is different. Thus it will be sometimes incorrect. I think that if `xSq` and `ySq` are close to each other and their sum is not exact, then the above algorithm will give you twice the rounding error.
Thanks Paul and Vincent for finding the issue with this!  Actually I was trying to implement the Fast2Sum version since we are in radix-2.  Moreover, since both `xSq` and `ySq` are non-negative:

```
max(xSq, ySq) <= sumSq <= sqrt(2) max(xSq, ySq)
```
and so `sumSq - max(xSq, ySq)` is exact.  I was trying to implement it without branch and thought that `(sumSq - min(xSq, ySq)) - max(xSq, ySq) = 0`, which could be easily disproved by the following example:  Consider single precision with `xSq = 1 + 2^(-23)` (I know it's not a square) and `ySq = 2^(-24)` with default rounding mode, then `sumSq = xSq + ySq = 1 + 2^(-22)`.  Then:
```
sumSq - xSq = 2^(-23)  and  sumSq - ySq = 1 + 2^(-22) = sumSq
```
And hence:
```
(sumSq - xSq) - ySq = 2^-24
(sumSq - ySq) - xSq = 2^-23
((sumSq - xSq) - ySq) + ((sumSq - ySq) - xSq) != 2^(-24)  which is the rounding error.
```

I change it back to the normal Fast2Sum implementation with branching now, so the rounding error computation should be correct.


Repository:
  rG LLVM Github Monorepo

CHANGES SINCE LAST ACTION
  https://reviews.llvm.org/D118157/new/

https://reviews.llvm.org/D118157