llvm-project

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

[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.