[libc-commits] [libc] [libc][math] Implement atan2f correctly rounded to all rounding modes. (PR #86716)
Nick Desaulniers via libc-commits
libc-commits at lists.llvm.org
Tue Mar 26 13:08:23 PDT 2024
================
@@ -0,0 +1,297 @@
+//===-- Single-precision atan2f function ----------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/atan2f.h"
+#include "inv_trigf_utils.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/except_value_utils.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/FPUtil/rounding_mode.h"
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+
+namespace LIBC_NAMESPACE {
+
+namespace {
+
+// Look up tables for accurate pass:
+
+// atan(i/16) with i = 0..16, generated by Sollya with:
+// > for i from 0 to 16 do {
+// a = round(atan(i/16), D, RN);
+// b = round(atan(i/16) - a, D, RN);
+// print("{", b, ",", a, "},");
+// };
+constexpr fputil::DoubleDouble ATAN_I[17] = {
+ {0.0, 0.0},
+ {-0x1.c934d86d23f1dp-60, 0x1.ff55bb72cfdeap-5},
+ {-0x1.cd37686760c17p-59, 0x1.fd5ba9aac2f6ep-4},
+ {0x1.347b0b4f881cap-58, 0x1.7b97b4bce5b02p-3},
+ {0x1.8ab6e3cf7afbdp-57, 0x1.f5b75f92c80ddp-3},
+ {-0x1.963a544b672d8p-57, 0x1.362773707ebccp-2},
+ {-0x1.c63aae6f6e918p-56, 0x1.6f61941e4def1p-2},
+ {-0x1.24dec1b50b7ffp-56, 0x1.a64eec3cc23fdp-2},
+ {0x1.a2b7f222f65e2p-56, 0x1.dac670561bb4fp-2},
+ {-0x1.d5b495f6349e6p-56, 0x1.0657e94db30dp-1},
+ {-0x1.928df287a668fp-58, 0x1.1e00babdefeb4p-1},
+ {0x1.1021137c71102p-55, 0x1.345f01cce37bbp-1},
+ {0x1.2419a87f2a458p-56, 0x1.4978fa3269ee1p-1},
+ {0x1.0028e4bc5e7cap-57, 0x1.5d58987169b18p-1},
+ {-0x1.8c34d25aadef6p-56, 0x1.700a7c5784634p-1},
+ {-0x1.bf76229d3b917p-56, 0x1.819d0b7158a4dp-1},
+ {0x1.1a62633145c07p-55, 0x1.921fb54442d18p-1},
+};
+
+// Taylor polynomial, generated by Sollya with:
+// > for i from 0 to 8 do {
+// j = (-1)^(i + 1)/(2*i + 1);
+// a = round(j, D, RN);
+// b = round(j - a, D, RN);
+// print("{", b, ",", a, "},");
+// };
+constexpr fputil::DoubleDouble COEFFS[9] = {
+ {0.0, 1.0}, // 1
+ {-0x1.5555555555555p-56, -0x1.5555555555555p-2}, // -1/3
+ {-0x1.999999999999ap-57, 0x1.999999999999ap-3}, // 1/5
+ {-0x1.2492492492492p-57, -0x1.2492492492492p-3}, // -1/7
+ {0x1.c71c71c71c71cp-58, 0x1.c71c71c71c71cp-4}, // 1/9
+ {0x1.745d1745d1746p-59, -0x1.745d1745d1746p-4}, // -1/11
+ {-0x1.3b13b13b13b14p-58, 0x1.3b13b13b13b14p-4}, // 1/13
+ {-0x1.1111111111111p-60, -0x1.1111111111111p-4}, // -1/15
+ {0x1.e1e1e1e1e1e1ep-61, 0x1.e1e1e1e1e1e1ep-5}, // 1/17
+};
+
+// Veltkamp's splitting of a double precision into hi + lo, where the hi part is
+// slightly smaller than an even split, so that the product of
+// hi * s * k is exact,
+// where:
+// s is single precsion,
+// 0 < k < 16 is an integer.
+// This is used when FMA instruction is not available.
+[[maybe_unused]] LIBC_INLINE constexpr fputil::DoubleDouble split_d(double a) {
+ fputil::DoubleDouble r{0.0, 0.0};
+ constexpr double C = 0x1.0p33 + 1.0;
+ double t1 = C * a;
+ double t2 = a - t1;
+ r.hi = t1 + t2;
+ r.lo = a - r.hi;
+ return r;
+}
+
+} // anonymous namespace
+
+// There are several range reduction steps we can take for atan2(y, x) as
+// follow:
+
+// * Range reduction 1: signness
+// atan2(y, x) will return a number between -PI and PI representing the angle
+// forming by the 0x axis and the vector (x, y) on the 0xy-plane.
+// In particular, we have that:
+// atan2(y, x) = atan( y/x ) if x >= 0 and y >= 0 (I-quadrant)
+// = pi + atan( y/x ) if x < 0 and y >= 0 (II-quadrant)
+// = -pi + atan( y/x ) if x < 0 and y < 0 (III-quadrant)
+// = atan( y/x ) if x >= 0 and y < 0 (IV-quadrant)
+// Since atan function is odd, we can use the formula:
+// atan(-u) = -atan(u)
+// to adjust the above conditions a bit further:
+// atan2(y, x) = atan( |y|/|x| ) if x >= 0 and y >= 0 (I-quadrant)
+// = pi - atan( |y|/|x| ) if x < 0 and y >= 0 (II-quadrant)
+// = -pi + atan( |y|/|x| ) if x < 0 and y < 0 (III-quadrant)
+// = -atan( |y|/|x| ) if x >= 0 and y < 0 (IV-quadrant)
+// Which can be simplified to:
+// atan2(y, x) = sign(y) * atan( |y|/|x| ) if x >= 0
+// = sign(y) * (pi - atan( |y|/|x| )) if x < 0
+
+// * Range reduction 2: reciprocal
+// Now that the argument inside atan is positive, we can use the formula:
+// atan(1/x) = pi/2 - atan(x)
+// to make the argument inside atan <= 1 as follow:
+// atan2(y, x) = sign(y) * atan( |y|/|x|) if 0 <= |y| <= x
+// = sign(y) * (pi/2 - atan( |x|/|y| ) if 0 <= x < |y|
+// = sign(y) * (pi - atan( |y|/|x| )) if 0 <= |y| <= -x
+// = sign(y) * (pi/2 + atan( |x|/|y| )) if 0 <= -x < |y|
+
+// * Range reduction 3: look up table.
+// After the previous two range reduction steps, we reduce the problem to
+// compute atan(u) with 0 <= u <= 1, or to be precise:
+// atan( n / d ) where n = min(|x|, |y|) and d = max(|x|, |y|).
+// An accurate polynomial approximation for the whole [0, 1] input range will
+// require a very large degree. To make it more efficient, we reduce the input
+// range further by finding an integer idx such that:
+// | n/d - idx/16 | <= 1/32.
+// In particular,
+// idx := 2^-4 * round(2^4 * n/d)
+// Then for the fast pass, we find a polynomial approximation for:
+// atan( n/d ) ~ atan( idx/16 ) + (n/d - idx/16) * Q(n/d - idx/16)
+// For the accurate pass, we use the addition formula:
+// atan( n/d ) - atan( idx/16 ) = atan( (n/d - idx/16)/(1 + (n*idx)/(16*d)) )
+// = atan( (n - d * idx/16)/(d + n * idx/16) )
+// And finally we use Taylor polynomial to compute the RHS in the accurate pass:
+// atan(u) ~ P(u) = u - u^3/3 + u^5/5 - u^7/7 + u^9/9 - u^11/11 + u^13/13 -
+// - u^15/15 + u^17/17
+// It's error in double-double precision is estimated in Sollya to be:
+// > P = x - x^3/3 + x^5/5 -x^7/7 + x^9/9 - x^11/11 + x^13/13 - x^15/15
+// + x^17/17;
+// > dirtyinfnorm(atan(x) - P, [-2^-5, 2^-5]);
+// 0x1.aec6f...p-100
+// which is about rounding errors of double-double (2^-104).
+
+LLVM_LIBC_FUNCTION(float, atan2f, (float y, float x)) {
+ using FPBits = typename fputil::FPBits<float>;
+ constexpr double IS_NEG[2] = {1.0, -1.0};
+ constexpr double PI = 0x1.921fb54442d18p1;
+ constexpr double PI_LO = 0x1.1a62633145c07p-53;
+ constexpr double PI_OVER_4 = 0x1.921fb54442d18p-1;
+ constexpr double PI_OVER_2 = 0x1.921fb54442d18p0;
+ constexpr double THREE_PI_OVER_4 = 0x1.2d97c7f3321d2p+1;
+ // Adjustment for constant term:
+ // CONST_ADJ[x_sign][y_sign][recip]
+ constexpr fputil::DoubleDouble CONST_ADJ[2][2][2] = {
+ {{{0.0, 0.0}, {-PI_LO / 2, -PI_OVER_2}},
+ {{-0.0, -0.0}, {-PI_LO / 2, -PI_OVER_2}}},
+ {{{-PI_LO, -PI}, {PI_LO / 2, PI_OVER_2}},
+ {{-PI_LO, -PI}, {PI_LO / 2, PI_OVER_2}}}};
+
+ FPBits x_bits(x), y_bits(y);
+ bool x_sign = x_bits.sign().is_neg();
+ bool y_sign = y_bits.sign().is_neg();
+ x_bits.set_sign(Sign::POS);
+ y_bits.set_sign(Sign::POS);
+ uint32_t x_abs = x_bits.uintval();
+ uint32_t y_abs = y_bits.uintval();
+ uint32_t max_abs = x_abs > y_abs ? x_abs : y_abs;
+ uint32_t min_abs = x_abs <= y_abs ? x_abs : y_abs;
+ bool recip = x_abs < y_abs;
----------------
nickdesaulniers wrote:
`recip` is not used until later, after the (unlikely) conditional that may return early without needing it. Consider sinking this closer to the use. Perhaps some of the expressions that the value of `recip` depends on can be sunk too, if they're also unused on the if statement starting on L173.
https://github.com/llvm/llvm-project/pull/86716
More information about the libc-commits
mailing list