1/*
2 * Single-precision scalar atan2(x) function.
3 *
4 * Copyright (c) 2021-2023, Arm Limited.
5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6 */
7
8#include <stdbool.h>
9
10#include "atanf_common.h"
11#include "math_config.h"
12#include "pl_sig.h"
13#include "pl_test.h"
14
15#define Pi (0x1.921fb6p+1f)
16#define PiOver2 (0x1.921fb6p+0f)
17#define PiOver4 (0x1.921fb6p-1f)
18#define SignMask (0x80000000)
19
20/* We calculate atan2f by P(n/d), where n and d are similar to the input
21   arguments, and P is a polynomial. The polynomial may underflow.
22   POLY_UFLOW_BOUND is the lower bound of the difference in exponents of n and d
23   for which P underflows, and is used to special-case such inputs.  */
24#define POLY_UFLOW_BOUND 24
25
26static inline int32_t
27biased_exponent (float f)
28{
29  uint32_t fi = asuint (f);
30  int32_t ex = (int32_t) ((fi & 0x7f800000) >> 23);
31  if (unlikely (ex == 0))
32    {
33      /* Subnormal case - we still need to get the exponent right for subnormal
34	 numbers as division may take us back inside the normal range.  */
35      return ex - __builtin_clz (fi << 9);
36    }
37  return ex;
38}
39
40/* Fast implementation of scalar atan2f. Largest observed error is
41   2.88ulps in [99.0, 101.0] x [99.0, 101.0]:
42   atan2f(0x1.9332d8p+6, 0x1.8cb6c4p+6) got 0x1.964646p-1
43				       want 0x1.964640p-1.  */
44float
45atan2f (float y, float x)
46{
47  uint32_t ix = asuint (x);
48  uint32_t iy = asuint (y);
49
50  uint32_t sign_x = ix & SignMask;
51  uint32_t sign_y = iy & SignMask;
52
53  uint32_t iax = ix & ~SignMask;
54  uint32_t iay = iy & ~SignMask;
55
56  /* x or y is NaN.  */
57  if ((iax > 0x7f800000) || (iay > 0x7f800000))
58    return x + y;
59
60  /* m = 2 * sign(x) + sign(y).  */
61  uint32_t m = ((iy >> 31) & 1) | ((ix >> 30) & 2);
62
63  /* The following follows glibc ieee754 implementation, except
64     that we do not use +-tiny shifts (non-nearest rounding mode).  */
65
66  int32_t exp_diff = biased_exponent (x) - biased_exponent (y);
67
68  /* Special case for (x, y) either on or very close to the x axis. Either y =
69     0, or y is tiny and x is huge (difference in exponents >=
70     POLY_UFLOW_BOUND). In the second case, we only want to use this special
71     case when x is negative (i.e. quadrants 2 or 3).  */
72  if (unlikely (iay == 0 || (exp_diff >= POLY_UFLOW_BOUND && m >= 2)))
73    {
74      switch (m)
75	{
76	case 0:
77	case 1:
78	  return y; /* atan(+-0,+anything)=+-0.  */
79	case 2:
80	  return Pi; /* atan(+0,-anything) = pi.  */
81	case 3:
82	  return -Pi; /* atan(-0,-anything) =-pi.  */
83	}
84    }
85  /* Special case for (x, y) either on or very close to the y axis. Either x =
86     0, or x is tiny and y is huge (difference in exponents >=
87     POLY_UFLOW_BOUND).  */
88  if (unlikely (iax == 0 || exp_diff <= -POLY_UFLOW_BOUND))
89    return sign_y ? -PiOver2 : PiOver2;
90
91  /* x is INF.  */
92  if (iax == 0x7f800000)
93    {
94      if (iay == 0x7f800000)
95	{
96	  switch (m)
97	    {
98	    case 0:
99	      return PiOver4; /* atan(+INF,+INF).  */
100	    case 1:
101	      return -PiOver4; /* atan(-INF,+INF).  */
102	    case 2:
103	      return 3.0f * PiOver4; /* atan(+INF,-INF).  */
104	    case 3:
105	      return -3.0f * PiOver4; /* atan(-INF,-INF).  */
106	    }
107	}
108      else
109	{
110	  switch (m)
111	    {
112	    case 0:
113	      return 0.0f; /* atan(+...,+INF).  */
114	    case 1:
115	      return -0.0f; /* atan(-...,+INF).  */
116	    case 2:
117	      return Pi; /* atan(+...,-INF).  */
118	    case 3:
119	      return -Pi; /* atan(-...,-INF).  */
120	    }
121	}
122    }
123  /* y is INF.  */
124  if (iay == 0x7f800000)
125    return sign_y ? -PiOver2 : PiOver2;
126
127  uint32_t sign_xy = sign_x ^ sign_y;
128
129  float ax = asfloat (iax);
130  float ay = asfloat (iay);
131
132  bool pred_aygtax = (ay > ax);
133
134  /* Set up z for call to atanf.  */
135  float n = pred_aygtax ? -ax : ay;
136  float d = pred_aygtax ? ay : ax;
137  float z = n / d;
138
139  float ret;
140  if (unlikely (m < 2 && exp_diff >= POLY_UFLOW_BOUND))
141    {
142      /* If (x, y) is very close to x axis and x is positive, the polynomial
143	 will underflow and evaluate to z.  */
144      ret = z;
145    }
146  else
147    {
148      /* Work out the correct shift.  */
149      float shift = sign_x ? -2.0f : 0.0f;
150      shift = pred_aygtax ? shift + 1.0f : shift;
151      shift *= PiOver2;
152
153      ret = eval_poly (z, z, shift);
154    }
155
156  /* Account for the sign of x and y.  */
157  return asfloat (asuint (ret) ^ sign_xy);
158}
159
160/* Arity of 2 means no mathbench entry emitted. See test/mathbench_funcs.h.  */
161PL_SIG (S, F, 2, atan2)
162PL_TEST_ULP (atan2f, 2.4)
163PL_TEST_INTERVAL (atan2f, -10.0, 10.0, 50000)
164PL_TEST_INTERVAL (atan2f, -1.0, 1.0, 40000)
165PL_TEST_INTERVAL (atan2f, 0.0, 1.0, 40000)
166PL_TEST_INTERVAL (atan2f, 1.0, 100.0, 40000)
167PL_TEST_INTERVAL (atan2f, 1e6, 1e32, 40000)
168