1 /* 2 * Double-precision atan(x) function. 3 * 4 * Copyright (c) 2022-2023, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6 */ 7 8 #include "pl_sig.h" 9 #include "pl_test.h" 10 #include "atan_common.h" 11 12 #define AbsMask 0x7fffffffffffffff 13 #define PiOver2 0x1.921fb54442d18p+0 14 #define TinyBound 0x3e1 /* top12(asuint64(0x1p-30)). */ 15 #define BigBound 0x434 /* top12(asuint64(0x1p53)). */ 16 #define OneTop 0x3ff 17 18 /* Fast implementation of double-precision atan. 19 Based on atan(x) ~ shift + z + z^3 * P(z^2) with reduction to [0,1] using 20 z=1/x and shift = pi/2. Maximum observed error is 2.27 ulps: 21 atan(0x1.0005af27c23e9p+0) got 0x1.9225645bdd7c1p-1 22 want 0x1.9225645bdd7c3p-1. */ 23 double 24 atan (double x) 25 { 26 uint64_t ix = asuint64 (x); 27 uint64_t sign = ix & ~AbsMask; 28 uint64_t ia = ix & AbsMask; 29 uint32_t ia12 = ia >> 52; 30 31 if (unlikely (ia12 >= BigBound || ia12 < TinyBound)) 32 { 33 if (ia12 < TinyBound) 34 /* Avoid underflow by returning x. */ 35 return x; 36 if (ia > 0x7ff0000000000000) 37 /* Propagate NaN. */ 38 return __math_invalid (x); 39 /* atan(x) rounds to PiOver2 for large x. */ 40 return asdouble (asuint64 (PiOver2) ^ sign); 41 } 42 43 double z, az, shift; 44 if (ia12 >= OneTop) 45 { 46 /* For x > 1, use atan(x) = pi / 2 + atan(-1 / x). */ 47 z = -1.0 / x; 48 shift = PiOver2; 49 /* Use absolute value only when needed (odd powers of z). */ 50 az = -fabs (z); 51 } 52 else 53 { 54 /* For x < 1, approximate atan(x) directly. */ 55 z = x; 56 shift = 0; 57 az = asdouble (ia); 58 } 59 60 /* Calculate polynomial, shift + z + z^3 * P(z^2). */ 61 double y = eval_poly (z, az, shift); 62 /* Copy sign. */ 63 return asdouble (asuint64 (y) ^ sign); 64 } 65 66 PL_SIG (S, D, 1, atan, -10.0, 10.0) 67 PL_TEST_ULP (atan, 1.78) 68 PL_TEST_INTERVAL (atan, 0, 0x1p-30, 10000) 69 PL_TEST_INTERVAL (atan, -0, -0x1p-30, 1000) 70 PL_TEST_INTERVAL (atan, 0x1p-30, 0x1p53, 900000) 71 PL_TEST_INTERVAL (atan, -0x1p-30, -0x1p53, 90000) 72 PL_TEST_INTERVAL (atan, 0x1p53, inf, 10000) 73 PL_TEST_INTERVAL (atan, -0x1p53, -inf, 1000) 74