1 /* 2 * Single-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 "atanf_common.h" 9 #include "pl_sig.h" 10 #include "pl_test.h" 11 12 #define PiOver2 0x1.921fb6p+0f 13 #define AbsMask 0x7fffffff 14 #define TinyBound 0x30800000 /* asuint(0x1p-30). */ 15 #define BigBound 0x4e800000 /* asuint(0x1p30). */ 16 #define One 0x3f800000 17 18 /* Approximation of single-precision atan(x) based on 19 atan(x) ~ shift + z + z^3 * P(z^2) with reduction to [0,1] 20 using z=-1/x and shift = pi/2. 21 Maximum error is 2.88 ulps: 22 atanf(0x1.0565ccp+0) got 0x1.97771p-1 23 want 0x1.97770ap-1. */ 24 float 25 atanf (float x) 26 { 27 uint32_t ix = asuint (x); 28 uint32_t sign = ix & ~AbsMask; 29 uint32_t ia = ix & AbsMask; 30 31 if (unlikely (ia < TinyBound)) 32 /* Avoid underflow by returning x. */ 33 return x; 34 35 if (unlikely (ia > BigBound)) 36 { 37 if (ia > 0x7f800000) 38 /* Propagate NaN. */ 39 return __math_invalidf (x); 40 /* atan(x) rounds to PiOver2 for large x. */ 41 return asfloat (asuint (PiOver2) ^ sign); 42 } 43 44 float z, az, shift; 45 if (ia > One) 46 { 47 /* For x > 1, use atan(x) = pi / 2 + atan(-1 / x). */ 48 z = -1.0f / x; 49 shift = PiOver2; 50 /* Use absolute value only when needed (odd powers of z). */ 51 az = -fabsf (z); 52 } 53 else 54 { 55 /* For x < 1, approximate atan(x) directly. */ 56 z = x; 57 az = asfloat (ia); 58 shift = 0; 59 } 60 61 /* Calculate polynomial, shift + z + z^3 * P(z^2). */ 62 float y = eval_poly (z, az, shift); 63 /* Copy sign. */ 64 return asfloat (asuint (y) ^ sign); 65 } 66 67 PL_SIG (S, F, 1, atan, -10.0, 10.0) 68 PL_TEST_ULP (atanf, 2.38) 69 PL_TEST_INTERVAL (atanf, 0, 0x1p-30, 5000) 70 PL_TEST_INTERVAL (atanf, -0, -0x1p-30, 5000) 71 PL_TEST_INTERVAL (atanf, 0x1p-30, 1, 40000) 72 PL_TEST_INTERVAL (atanf, -0x1p-30, -1, 40000) 73 PL_TEST_INTERVAL (atanf, 1, 0x1p30, 40000) 74 PL_TEST_INTERVAL (atanf, -1, -0x1p30, 40000) 75 PL_TEST_INTERVAL (atanf, 0x1p30, inf, 1000) 76 PL_TEST_INTERVAL (atanf, -0x1p30, -inf, 1000) 77