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 26 static inline int32_t 27 biased_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. */ 44 float 45 atan2f (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. */ 161 PL_SIG (S, F, 2, atan2) 162 PL_TEST_ULP (atan2f, 2.4) 163 PL_TEST_INTERVAL (atan2f, -10.0, 10.0, 50000) 164 PL_TEST_INTERVAL (atan2f, -1.0, 1.0, 40000) 165 PL_TEST_INTERVAL (atan2f, 0.0, 1.0, 40000) 166 PL_TEST_INTERVAL (atan2f, 1.0, 100.0, 40000) 167 PL_TEST_INTERVAL (atan2f, 1e6, 1e32, 40000) 168