[libc-commits] [libc] [libc][math] Add float-only option for atan2f. (PR #122979)
Nick Desaulniers via libc-commits
libc-commits at lists.llvm.org
Tue Feb 11 11:27:56 PST 2025
================
@@ -0,0 +1,239 @@
+//===-- 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/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/double_double.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/config.h"
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/math/atan2f.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace {
+
+using FloatFloat = fputil::FloatFloat;
+
+// atan(i/64) with i = 0..16, generated by Sollya with:
+// > for i from 0 to 16 do {
+// a = round(atan(i/16), SG, RN);
+// b = round(atan(i/16) - a, SG, RN);
+// print("{", b, ",", a, "},");
+// };
+constexpr FloatFloat ATAN_I[17] = {
+ {0.0f, 0.0f},
+ {-0x1.1a6042p-30f, 0x1.ff55bcp-5f},
+ {-0x1.54f424p-30f, 0x1.fd5baap-4f},
+ {0x1.79cb6p-28f, 0x1.7b97b4p-3f},
+ {-0x1.b4dfc8p-29f, 0x1.f5b76p-3f},
+ {-0x1.1f0286p-27f, 0x1.362774p-2f},
+ {0x1.e4defp-30f, 0x1.6f6194p-2f},
+ {0x1.e611fep-29f, 0x1.a64eecp-2f},
+ {0x1.586ed4p-28f, 0x1.dac67p-2f},
+ {-0x1.6499e6p-26f, 0x1.0657eap-1f},
+ {0x1.7bdfd6p-26f, 0x1.1e00bap-1f},
+ {-0x1.98e422p-28f, 0x1.345f02p-1f},
+ {0x1.934f7p-28f, 0x1.4978fap-1f},
+ {0x1.c5a6c6p-27f, 0x1.5d5898p-1f},
+ {0x1.5e118cp-27f, 0x1.700a7cp-1f},
+ {-0x1.1d4eb6p-26f, 0x1.819d0cp-1f},
+ {-0x1.777a5cp-26f, 0x1.921fb6p-1f},
+};
+
+// Approximate atan(x) for |x| <= 2^-5.
+// Using degree-3 Taylor polynomial:
+// P = x - x^3/3
+// Then the absolute error is bounded by:
+// |atan(x) - P(x)| < |x|^5/5 < 2^(-5*5) / 5 < 2^-27.
+// And the relative error is bounded by:
+// |(atan(x) - P(x))/atan(x)| < |x|^4 / 4 < 2^-22.
+// For x = x_hi + x_lo, fully expand the polynomial and drop any terms less than
+// ulp(x_hi^3 / 3) gives us:
+// P(x) ~ x_hi - x_hi^3/3 + x_lo * (1 - x_hi^2)
+FloatFloat atan_eval(const FloatFloat &x) {
+ FloatFloat p;
+ p.hi = x.hi;
+ float x_hi_sq = x.hi * x.hi;
+ // c0 ~ - x_hi^2 / 3
+ float c0 = -0x1.555556p-2f * x_hi_sq;
+ // c1 ~ x_lo * (1 - x_hi^2)
+ float c1 = fputil::multiply_add(x_hi_sq, -x.lo, x.lo);
+ // p.lo ~ - x_hi^3 / 3 + x_lo * (1 - x_hi*2)
+ p.lo = fputil::multiply_add(x.hi, c0, c1);
+ return p;
+}
+
+} // 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)
+// with Q(x) = x - x^3/3 be the cubic Taylor polynomial of atan(x).
+// It's error in float-float precision is estimated in Sollya to be:
+// > P = x - x^3/3;
+// > dirtyinfnorm(atan(x) - P, [-2^-5, 2^-5]);
+// 0x1.995...p-28.
+
+LLVM_LIBC_FUNCTION(float, atan2f, (float y, float x)) {
+ using FPBits = typename fputil::FPBits<float>;
+ constexpr float IS_NEG[2] = {1.0f, -1.0f};
+ constexpr FloatFloat ZERO = {0.0f, 0.0f};
+ constexpr FloatFloat MZERO = {-0.0f, -0.0f};
+ constexpr FloatFloat PI = {-0x1.777a5cp-24f, 0x1.921fb6p1f};
+ constexpr FloatFloat MPI = {0x1.777a5cp-24f, -0x1.921fb6p1f};
+ constexpr FloatFloat PI_OVER_4 = {-0x1.777a5cp-26f, 0x1.921fb6p-1f};
+ constexpr FloatFloat PI_OVER_2 = {-0x1.777a5cp-25f, 0x1.921fb6p0f};
+ constexpr FloatFloat MPI_OVER_2 = {-0x1.777a5cp-25f, 0x1.921fb6p0f};
+ constexpr FloatFloat THREE_PI_OVER_4 = {-0x1.99bc5cp-28f, 0x1.2d97c8p1f};
+ // Adjustment for constant term:
+ // CONST_ADJ[x_sign][y_sign][recip]
+ constexpr FloatFloat CONST_ADJ[2][2][2] = {
+ {{ZERO, MPI_OVER_2}, {MZERO, MPI_OVER_2}},
+ {{MPI, PI_OVER_2}, {MPI, 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 = x_bits.abs();
+ y_bits = y_bits.abs();
+ uint32_t x_abs = x_bits.uintval();
+ uint32_t y_abs = y_bits.uintval();
+ bool recip = x_abs < y_abs;
+ uint32_t min_abs = recip ? x_abs : y_abs;
+ uint32_t max_abs = !recip ? x_abs : y_abs;
+ unsigned min_exp = static_cast<unsigned>(min_abs >> FPBits::FRACTION_LEN);
+ unsigned max_exp = static_cast<unsigned>(max_abs >> FPBits::FRACTION_LEN);
----------------
nickdesaulniers wrote:
A common convention used throughout llvm on when to use `auto` is when the resulting type would be specified on the right hand side of an assignment due to a cast.
https://llvm.org/docs/CodingStandards.html#use-auto-type-deduction-to-make-code-more-readable
```suggestion
auto min_exp = static_cast<unsigned>(min_abs >> FPBits::FRACTION_LEN);
auto max_exp = static_cast<unsigned>(max_abs >> FPBits::FRACTION_LEN);
```
https://github.com/llvm/llvm-project/pull/122979
More information about the libc-commits
mailing list