[all-commits] [llvm/llvm-project] 7fc9fb: [libc][math] Implement double precision cbrt corre...

[lntue via All-commits]

  Branch: refs/heads/main
  Home:   https://github.com/llvm/llvm-project
  Commit: 7fc9fb9f3f671ee4b1ccfefaf03ed18cc0c3e3c3
      https://github.com/llvm/llvm-project/commit/7fc9fb9f3f671ee4b1ccfefaf03ed18cc0c3e3c3
  Author: lntue <35648136+lntue at users.noreply.github.com>
  Date:   2024-07-17 (Wed, 17 Jul 2024)

  Changed paths:
    M libc/config/darwin/arm/entrypoints.txt
    M libc/config/gpu/entrypoints.txt
    M libc/config/linux/aarch64/entrypoints.txt
    M libc/config/linux/arm/entrypoints.txt
    M libc/config/linux/riscv/entrypoints.txt
    M libc/config/linux/x86_64/entrypoints.txt
    M libc/config/windows/entrypoints.txt
    M libc/docs/math/index.rst
    M libc/spec/stdc.td
    M libc/src/math/CMakeLists.txt
    A libc/src/math/cbrt.h
    M libc/src/math/generic/CMakeLists.txt
    A libc/src/math/generic/cbrt.cpp
    M libc/test/src/math/CMakeLists.txt
    A libc/test/src/math/cbrt_test.cpp
    M libc/test/src/math/smoke/CMakeLists.txt
    A libc/test/src/math/smoke/cbrt_test.cpp

  Log Message:
  -----------
  [libc][math] Implement double precision cbrt correctly rounded to all rounding modes. (#99262)

Division-less Newton iterations algorithm for cube roots.

1. **Range reduction**

For `x = (-1)^s * 2^e * (1.m)`, we get 2 reduced arguments `x_r` and `a`
as:
```
  x_r = 1.m
  a   = (-1)^s * 2^(e % 3) * (1.m)
```
Then `cbrt(x) = x^(1/3)` can be computed as:
```
  x^(1/3) = 2^(e / 3) * a^(1/3).
```

In order to avoid division, we compute `a^(-2/3)` using Newton method
and then
multiply the results by a:
```
  a^(1/3) = a * a^(-2/3).
```

2. **First approximation to a^(-2/3)**

First, we use a degree-7 minimax polynomial generated by Sollya to
approximate `x_r^(-2/3)` for `1 <= x_r < 2`.
```
  p = P(x_r) ~ x_r^(-2/3),
```
with relative errors bounded by:
```
  | p / x_r^(-2/3) - 1 | < 1.16 * 2^-21.
```

Then we multiply with `2^(e % 3)` from a small lookup table to get:
```
  x_0 = 2^(-2*(e % 3)/3) * p
      ~ 2^(-2*(e % 3)/3) * x_r^(-2/3)
      = a^(-2/3)
```
with relative errors:
```
  | x_0 / a^(-2/3) - 1 | < 1.16 * 2^-21.
```
This step is done in double precision.

3. **First Newton iteration**

We follow the method described in:
Sibidanov, A. and Zimmermann, P., "Correctly rounded cubic root
evaluation
in double precision", https://core-math.gitlabpages.inria.fr/cbrt64.pdf
to derive multiplicative Newton iterations as below:
Let `x_n` be the nth approximation to `a^(-2/3)`. Define the n^th error
as:
```
  h_n = x_n^3 * a^2 - 1
```
Then:
```
  a^(-2/3) = x_n / (1 + h_n)^(1/3)
           = x_n * (1 - (1/3) * h_n + (2/9) * h_n^2 - (14/81) * h_n^3 + ...)
```
using the Taylor series expansion of `(1 + h_n)^(-1/3)`.

Apply to `x_0` above:
```
  h_0 = x_0^3 * a^2 - 1
      = a^2 * (x_0 - a^(-2/3)) * (x_0^2 + x_0 * a^(-2/3) + a^(-4/3)),
```
it's bounded by:
```
  |h_0| < 4 * 3 * 1.16 * 2^-21 * 4 < 2^-17.
```
So in the first iteration step, we use:
```
  x_1 = x_0 * (1 - (1/3) * h_n + (2/9) * h_n^2 - (14/81) * h_n^3)
```
Its relative error is bounded by:
```
  | x_1 / a^(-2/3) - 1 | < 35/242 * |h_0|^4 < 2^-70.
```
Then we perform Ziv's rounding test and check if the answer is exact.
This step is done in double-double precision.

4. **Second Newton iteration**

If the Ziv's rounding test from the previous step fails, we define the
error
term:
```
  h_1 = x_1^3 * a^2 - 1,
```
And perform another iteration:
```
  x_2 = x_1 * (1 - h_1 / 3)
```
with the relative errors exceed the precision of double-double.
We then check the Ziv's accuracy test with relative errors < 2^-102 to
compensate for rounding errors.

5. **Final iteration**
 
If the Ziv's accuracy test from the previous step fails, we perform
another
iteration in 128-bit precision and check for exact outputs.



To unsubscribe from these emails, change your notification settings at https://github.com/llvm/llvm-project/settings/notifications