aboutsummaryrefslogtreecommitdiff
path: root/pl/math/atan2f_3u.c
diff options
context:
space:
mode:
Diffstat (limited to 'pl/math/atan2f_3u.c')
-rw-r--r--pl/math/atan2f_3u.c167
1 files changed, 0 insertions, 167 deletions
diff --git a/pl/math/atan2f_3u.c b/pl/math/atan2f_3u.c
deleted file mode 100644
index 38e1df59c102..000000000000
--- a/pl/math/atan2f_3u.c
+++ /dev/null
@@ -1,167 +0,0 @@
-/*
- * Single-precision scalar atan2(x) function.
- *
- * Copyright (c) 2021-2023, Arm Limited.
- * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
- */
-
-#include <stdbool.h>
-
-#include "atanf_common.h"
-#include "math_config.h"
-#include "pl_sig.h"
-#include "pl_test.h"
-
-#define Pi (0x1.921fb6p+1f)
-#define PiOver2 (0x1.921fb6p+0f)
-#define PiOver4 (0x1.921fb6p-1f)
-#define SignMask (0x80000000)
-
-/* We calculate atan2f by P(n/d), where n and d are similar to the input
- arguments, and P is a polynomial. The polynomial may underflow.
- POLY_UFLOW_BOUND is the lower bound of the difference in exponents of n and d
- for which P underflows, and is used to special-case such inputs. */
-#define POLY_UFLOW_BOUND 24
-
-static inline int32_t
-biased_exponent (float f)
-{
- uint32_t fi = asuint (f);
- int32_t ex = (int32_t) ((fi & 0x7f800000) >> 23);
- if (unlikely (ex == 0))
- {
- /* Subnormal case - we still need to get the exponent right for subnormal
- numbers as division may take us back inside the normal range. */
- return ex - __builtin_clz (fi << 9);
- }
- return ex;
-}
-
-/* Fast implementation of scalar atan2f. Largest observed error is
- 2.88ulps in [99.0, 101.0] x [99.0, 101.0]:
- atan2f(0x1.9332d8p+6, 0x1.8cb6c4p+6) got 0x1.964646p-1
- want 0x1.964640p-1. */
-float
-atan2f (float y, float x)
-{
- uint32_t ix = asuint (x);
- uint32_t iy = asuint (y);
-
- uint32_t sign_x = ix & SignMask;
- uint32_t sign_y = iy & SignMask;
-
- uint32_t iax = ix & ~SignMask;
- uint32_t iay = iy & ~SignMask;
-
- /* x or y is NaN. */
- if ((iax > 0x7f800000) || (iay > 0x7f800000))
- return x + y;
-
- /* m = 2 * sign(x) + sign(y). */
- uint32_t m = ((iy >> 31) & 1) | ((ix >> 30) & 2);
-
- /* The following follows glibc ieee754 implementation, except
- that we do not use +-tiny shifts (non-nearest rounding mode). */
-
- int32_t exp_diff = biased_exponent (x) - biased_exponent (y);
-
- /* Special case for (x, y) either on or very close to the x axis. Either y =
- 0, or y is tiny and x is huge (difference in exponents >=
- POLY_UFLOW_BOUND). In the second case, we only want to use this special
- case when x is negative (i.e. quadrants 2 or 3). */
- if (unlikely (iay == 0 || (exp_diff >= POLY_UFLOW_BOUND && m >= 2)))
- {
- switch (m)
- {
- case 0:
- case 1:
- return y; /* atan(+-0,+anything)=+-0. */
- case 2:
- return Pi; /* atan(+0,-anything) = pi. */
- case 3:
- return -Pi; /* atan(-0,-anything) =-pi. */
- }
- }
- /* Special case for (x, y) either on or very close to the y axis. Either x =
- 0, or x is tiny and y is huge (difference in exponents >=
- POLY_UFLOW_BOUND). */
- if (unlikely (iax == 0 || exp_diff <= -POLY_UFLOW_BOUND))
- return sign_y ? -PiOver2 : PiOver2;
-
- /* x is INF. */
- if (iax == 0x7f800000)
- {
- if (iay == 0x7f800000)
- {
- switch (m)
- {
- case 0:
- return PiOver4; /* atan(+INF,+INF). */
- case 1:
- return -PiOver4; /* atan(-INF,+INF). */
- case 2:
- return 3.0f * PiOver4; /* atan(+INF,-INF). */
- case 3:
- return -3.0f * PiOver4; /* atan(-INF,-INF). */
- }
- }
- else
- {
- switch (m)
- {
- case 0:
- return 0.0f; /* atan(+...,+INF). */
- case 1:
- return -0.0f; /* atan(-...,+INF). */
- case 2:
- return Pi; /* atan(+...,-INF). */
- case 3:
- return -Pi; /* atan(-...,-INF). */
- }
- }
- }
- /* y is INF. */
- if (iay == 0x7f800000)
- return sign_y ? -PiOver2 : PiOver2;
-
- uint32_t sign_xy = sign_x ^ sign_y;
-
- float ax = asfloat (iax);
- float ay = asfloat (iay);
-
- bool pred_aygtax = (ay > ax);
-
- /* Set up z for call to atanf. */
- float n = pred_aygtax ? -ax : ay;
- float d = pred_aygtax ? ay : ax;
- float z = n / d;
-
- float ret;
- if (unlikely (m < 2 && exp_diff >= POLY_UFLOW_BOUND))
- {
- /* If (x, y) is very close to x axis and x is positive, the polynomial
- will underflow and evaluate to z. */
- ret = z;
- }
- else
- {
- /* Work out the correct shift. */
- float shift = sign_x ? -2.0f : 0.0f;
- shift = pred_aygtax ? shift + 1.0f : shift;
- shift *= PiOver2;
-
- ret = eval_poly (z, z, shift);
- }
-
- /* Account for the sign of x and y. */
- return asfloat (asuint (ret) ^ sign_xy);
-}
-
-/* Arity of 2 means no mathbench entry emitted. See test/mathbench_funcs.h. */
-PL_SIG (S, F, 2, atan2)
-PL_TEST_ULP (atan2f, 2.4)
-PL_TEST_INTERVAL (atan2f, -10.0, 10.0, 50000)
-PL_TEST_INTERVAL (atan2f, -1.0, 1.0, 40000)
-PL_TEST_INTERVAL (atan2f, 0.0, 1.0, 40000)
-PL_TEST_INTERVAL (atan2f, 1.0, 100.0, 40000)
-PL_TEST_INTERVAL (atan2f, 1e6, 1e32, 40000)