1 /* 2 * Double-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 "atan_common.h" 11 #include "math_config.h" 12 #include "pl_sig.h" 13 #include "pl_test.h" 14 15 #define Pi (0x1.921fb54442d18p+1) 16 #define PiOver2 (0x1.921fb54442d18p+0) 17 #define PiOver4 (0x1.921fb54442d18p-1) 18 #define SignMask (0x8000000000000000) 19 #define ExpMask (0x7ff0000000000000) 20 21 /* We calculate atan2 by P(n/d), where n and d are similar to the input 22 arguments, and P is a polynomial. Evaluating P(x) requires calculating x^8, 23 which may underflow if n and d have very different magnitude. 24 POW8_EXP_UFLOW_BOUND is the lower bound of the difference in exponents of n 25 and d for which P underflows, and is used to special-case such inputs. */ 26 #define POW8_EXP_UFLOW_BOUND 62 27 28 static inline int64_t 29 biased_exponent (double f) 30 { 31 uint64_t fi = asuint64 (f); 32 return (fi & ExpMask) >> 52; 33 } 34 35 /* Fast implementation of scalar atan2. Largest errors are when y and x are 36 close together. The greatest observed error is 2.28 ULP: 37 atan2(-0x1.5915b1498e82fp+732, 0x1.54d11ef838826p+732) 38 got -0x1.954f42f1fa841p-1 want -0x1.954f42f1fa843p-1. */ 39 double 40 atan2 (double y, double x) 41 { 42 uint64_t ix = asuint64 (x); 43 uint64_t iy = asuint64 (y); 44 45 uint64_t sign_x = ix & SignMask; 46 uint64_t sign_y = iy & SignMask; 47 48 uint64_t iax = ix & ~SignMask; 49 uint64_t iay = iy & ~SignMask; 50 51 bool xisnan = isnan (x); 52 if (unlikely (isnan (y) && !xisnan)) 53 return __math_invalid (y); 54 if (unlikely (xisnan)) 55 return __math_invalid (x); 56 57 /* m = 2 * sign(x) + sign(y). */ 58 uint32_t m = ((iy >> 63) & 1) | ((ix >> 62) & 2); 59 60 int64_t exp_diff = biased_exponent (x) - biased_exponent (y); 61 62 /* y = 0. */ 63 if (iay == 0) 64 { 65 switch (m) 66 { 67 case 0: 68 case 1: 69 return y; /* atan(+-0,+anything)=+-0. */ 70 case 2: 71 return Pi; /* atan(+0,-anything) = pi. */ 72 case 3: 73 return -Pi; /* atan(-0,-anything) =-pi. */ 74 } 75 } 76 /* Special case for (x, y) either on or very close to the y axis. Either x = 77 0, or y is much larger than x (difference in exponents >= 78 POW8_EXP_UFLOW_BOUND). */ 79 if (unlikely (iax == 0 || exp_diff <= -POW8_EXP_UFLOW_BOUND)) 80 return sign_y ? -PiOver2 : PiOver2; 81 82 /* Special case for either x is INF or (x, y) is very close to x axis and x is 83 negative. */ 84 if (unlikely (iax == 0x7ff0000000000000 85 || (exp_diff >= POW8_EXP_UFLOW_BOUND && m >= 2))) 86 { 87 if (iay == 0x7ff0000000000000) 88 { 89 switch (m) 90 { 91 case 0: 92 return PiOver4; /* atan(+INF,+INF). */ 93 case 1: 94 return -PiOver4; /* atan(-INF,+INF). */ 95 case 2: 96 return 3.0 * PiOver4; /* atan(+INF,-INF). */ 97 case 3: 98 return -3.0 * PiOver4; /* atan(-INF,-INF). */ 99 } 100 } 101 else 102 { 103 switch (m) 104 { 105 case 0: 106 return 0.0; /* atan(+...,+INF). */ 107 case 1: 108 return -0.0; /* atan(-...,+INF). */ 109 case 2: 110 return Pi; /* atan(+...,-INF). */ 111 case 3: 112 return -Pi; /* atan(-...,-INF). */ 113 } 114 } 115 } 116 /* y is INF. */ 117 if (iay == 0x7ff0000000000000) 118 return sign_y ? -PiOver2 : PiOver2; 119 120 uint64_t sign_xy = sign_x ^ sign_y; 121 122 double ax = asdouble (iax); 123 double ay = asdouble (iay); 124 uint64_t pred_aygtax = (ay > ax); 125 126 /* Set up z for call to atan. */ 127 double n = pred_aygtax ? -ax : ay; 128 double d = pred_aygtax ? ay : ax; 129 double z = n / d; 130 131 double ret; 132 if (unlikely (m < 2 && exp_diff >= POW8_EXP_UFLOW_BOUND)) 133 { 134 /* If (x, y) is very close to x axis and x is positive, the polynomial 135 will underflow and evaluate to z. */ 136 ret = z; 137 } 138 else 139 { 140 /* Work out the correct shift. */ 141 double shift = sign_x ? -2.0 : 0.0; 142 shift = pred_aygtax ? shift + 1.0 : shift; 143 shift *= PiOver2; 144 145 ret = eval_poly (z, z, shift); 146 } 147 148 /* Account for the sign of x and y. */ 149 return asdouble (asuint64 (ret) ^ sign_xy); 150 } 151 152 /* Arity of 2 means no mathbench entry emitted. See test/mathbench_funcs.h. */ 153 PL_SIG (S, D, 2, atan2) 154 PL_TEST_ULP (atan2, 1.78) 155 PL_TEST_INTERVAL (atan2, -10.0, 10.0, 50000) 156 PL_TEST_INTERVAL (atan2, -1.0, 1.0, 40000) 157 PL_TEST_INTERVAL (atan2, 0.0, 1.0, 40000) 158 PL_TEST_INTERVAL (atan2, 1.0, 100.0, 40000) 159 PL_TEST_INTERVAL (atan2, 1e6, 1e32, 40000) 160