xref: /freebsd/contrib/arm-optimized-routines/pl/math/atanf_common.h (revision 5a02ffc32e777041dd2dad4e651ed2a0865a0a5d)
1072a4ba8SAndrew Turner /*
2*5a02ffc3SAndrew Turner  * Single-precision polynomial evaluation function for scalar
3072a4ba8SAndrew Turner  * atan(x) and atan2(y,x).
4072a4ba8SAndrew Turner  *
5072a4ba8SAndrew Turner  * Copyright (c) 2021-2023, Arm Limited.
6072a4ba8SAndrew Turner  * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
7072a4ba8SAndrew Turner  */
8072a4ba8SAndrew Turner 
9072a4ba8SAndrew Turner #ifndef PL_MATH_ATANF_COMMON_H
10072a4ba8SAndrew Turner #define PL_MATH_ATANF_COMMON_H
11072a4ba8SAndrew Turner 
12072a4ba8SAndrew Turner #include "math_config.h"
13*5a02ffc3SAndrew Turner #include "poly_scalar_f32.h"
14072a4ba8SAndrew Turner 
15072a4ba8SAndrew Turner /* Polynomial used in fast atanf(x) and atan2f(y,x) implementations
16072a4ba8SAndrew Turner    The order 7 polynomial P approximates (atan(sqrt(x))-sqrt(x))/x^(3/2).  */
17*5a02ffc3SAndrew Turner static inline float
eval_poly(float z,float az,float shift)18*5a02ffc3SAndrew Turner eval_poly (float z, float az, float shift)
19072a4ba8SAndrew Turner {
20072a4ba8SAndrew Turner   /* Use 2-level Estrin scheme for P(z^2) with deg(P)=7. However,
21072a4ba8SAndrew Turner      a standard implementation using z8 creates spurious underflow
22072a4ba8SAndrew Turner      in the very last fma (when z^8 is small enough).
23072a4ba8SAndrew Turner      Therefore, we split the last fma into a mul and and an fma.
24072a4ba8SAndrew Turner      Horner and single-level Estrin have higher errors that exceed
25072a4ba8SAndrew Turner      threshold.  */
26*5a02ffc3SAndrew Turner   float z2 = z * z;
27*5a02ffc3SAndrew Turner   float z4 = z2 * z2;
28072a4ba8SAndrew Turner 
29072a4ba8SAndrew Turner   /* Then assemble polynomial.  */
30*5a02ffc3SAndrew Turner   float y = fmaf (
31*5a02ffc3SAndrew Turner       z4, z4 * pairwise_poly_3_f32 (z2, z4, __atanf_poly_data.poly + 4),
32*5a02ffc3SAndrew Turner       pairwise_poly_3_f32 (z2, z4, __atanf_poly_data.poly));
33072a4ba8SAndrew Turner   /* Finalize:
34072a4ba8SAndrew Turner      y = shift + z * P(z^2).  */
35*5a02ffc3SAndrew Turner   return fmaf (y, z2 * az, az) + shift;
36072a4ba8SAndrew Turner }
37072a4ba8SAndrew Turner 
38072a4ba8SAndrew Turner #endif // PL_MATH_ATANF_COMMON_H
39