diff options
29 files changed, 269 insertions, 296 deletions
diff --git a/SHARED-FILES b/SHARED-FILES index b9627afdfe..99f8554393 100644 --- a/SHARED-FILES +++ b/SHARED-FILES @@ -314,3 +314,7 @@ sysdeps/ieee754/flt-32/s_atanf.c: (src/binary32/atan/atanf.c in CORE-MATH) - The code was adapted to use glibc code style and internal functions to handle errno, overflow, and underflow. +sysdeps/ieee754/flt-32/e_atan2f.c: + (src/binary32/atan2/atan2f.c in CORE-MATH) + - The code was adapted to use glibc code style and internal + functions to handle errno, overflow, and underflow. diff --git a/sysdeps/aarch64/libm-test-ulps b/sysdeps/aarch64/libm-test-ulps index 33d282231a..931fac6bb0 100644 --- a/sysdeps/aarch64/libm-test-ulps +++ b/sysdeps/aarch64/libm-test-ulps @@ -142,7 +142,6 @@ double: 1 ldouble: 1 Function: "atan2": -float: 1 ldouble: 2 Function: "atan2_advsimd": @@ -151,7 +150,6 @@ float: 2 Function: "atan2_downward": double: 1 -float: 2 ldouble: 2 Function: "atan2_sve": @@ -160,12 +158,10 @@ float: 2 Function: "atan2_towardzero": double: 1 -float: 2 ldouble: 3 Function: "atan2_upward": double: 1 -float: 1 ldouble: 2 Function: "atan2pi": @@ -354,22 +350,18 @@ ldouble: 4 Function: "carg": double: 1 -float: 1 ldouble: 2 Function: "carg_downward": double: 1 -float: 2 ldouble: 2 Function: "carg_towardzero": double: 1 -float: 2 ldouble: 3 Function: "carg_upward": double: 1 -float: 1 ldouble: 2 Function: Real part of "casin": diff --git a/sysdeps/alpha/fpu/libm-test-ulps b/sysdeps/alpha/fpu/libm-test-ulps index f9c1cf7cf5..d098eaf820 100644 --- a/sysdeps/alpha/fpu/libm-test-ulps +++ b/sysdeps/alpha/fpu/libm-test-ulps @@ -70,22 +70,18 @@ double: 1 ldouble: 1 Function: "atan2": -float: 2 ldouble: 2 Function: "atan2_downward": double: 1 -float: 2 ldouble: 2 Function: "atan2_towardzero": double: 1 -float: 2 ldouble: 3 Function: "atan2_upward": double: 1 -float: 2 ldouble: 2 Function: "atan_downward": @@ -217,22 +213,18 @@ float: 2 ldouble: 4 Function: "carg": -float: 1 ldouble: 2 Function: "carg_downward": double: 1 -float: 2 ldouble: 2 Function: "carg_towardzero": double: 1 -float: 2 ldouble: 3 Function: "carg_upward": double: 1 -float: 2 ldouble: 2 Function: Real part of "casin": diff --git a/sysdeps/arc/fpu/libm-test-ulps b/sysdeps/arc/fpu/libm-test-ulps index 37b0efae66..5b89ee2696 100644 --- a/sysdeps/arc/fpu/libm-test-ulps +++ b/sysdeps/arc/fpu/libm-test-ulps @@ -54,19 +54,15 @@ double: 1 Function: "atan2": double: 7 -float: 2 Function: "atan2_downward": double: 5 -float: 2 Function: "atan2_towardzero": double: 5 -float: 2 Function: "atan2_upward": double: 8 -float: 2 Function: "atan_downward": double: 1 @@ -173,19 +169,15 @@ float: 3 Function: "carg": double: 7 -float: 2 Function: "carg_downward": double: 5 -float: 2 Function: "carg_towardzero": double: 5 -float: 2 Function: "carg_upward": double: 8 -float: 2 Function: Real part of "casin": double: 3 diff --git a/sysdeps/arc/nofpu/libm-test-ulps b/sysdeps/arc/nofpu/libm-test-ulps index 8d283f0627..ac43ab8dcc 100644 --- a/sysdeps/arc/nofpu/libm-test-ulps +++ b/sysdeps/arc/nofpu/libm-test-ulps @@ -17,7 +17,6 @@ Function: "atan": double: 1 Function: "atan2": -float: 2 Function: "atanh": double: 2 @@ -43,7 +42,6 @@ double: 1 float: 2 Function: "carg": -float: 1 Function: Real part of "casin": double: 1 diff --git a/sysdeps/arm/libm-test-ulps b/sysdeps/arm/libm-test-ulps index 7b1a541378..234609cae5 100644 --- a/sysdeps/arm/libm-test-ulps +++ b/sysdeps/arm/libm-test-ulps @@ -85,19 +85,15 @@ Function: "atan": double: 1 Function: "atan2": -float: 2 Function: "atan2_downward": double: 1 -float: 2 Function: "atan2_towardzero": double: 1 -float: 2 Function: "atan2_upward": double: 1 -float: 2 Function: "atan2pi": double: 1 @@ -233,19 +229,15 @@ double: 3 float: 2 Function: "carg": -float: 1 Function: "carg_downward": double: 1 -float: 2 Function: "carg_towardzero": double: 1 -float: 2 Function: "carg_upward": double: 1 -float: 2 Function: Real part of "casin": double: 1 diff --git a/sysdeps/csky/fpu/libm-test-ulps b/sysdeps/csky/fpu/libm-test-ulps index 9d3fcf693d..b8de7af6b9 100644 --- a/sysdeps/csky/fpu/libm-test-ulps +++ b/sysdeps/csky/fpu/libm-test-ulps @@ -50,19 +50,15 @@ double: 3 Function: "atan": Function: "atan2": -float: 1 Function: "atan2_downward": double: 1 -float: 2 Function: "atan2_towardzero": double: 1 -float: 2 Function: "atan2_upward": double: 1 -float: 2 Function: "atan_downward": double: 1 @@ -166,19 +162,15 @@ double: 3 float: 2 Function: "carg": -float: 1 Function: "carg_downward": double: 1 -float: 2 Function: "carg_towardzero": double: 1 -float: 2 Function: "carg_upward": double: 1 -float: 2 Function: Real part of "casin": double: 1 diff --git a/sysdeps/csky/nofpu/libm-test-ulps b/sysdeps/csky/nofpu/libm-test-ulps index 1bab8effc7..a7594b6b1c 100644 --- a/sysdeps/csky/nofpu/libm-test-ulps +++ b/sysdeps/csky/nofpu/libm-test-ulps @@ -50,19 +50,15 @@ double: 3 Function: "atan": Function: "atan2": -float: 1 Function: "atan2_downward": double: 1 -float: 2 Function: "atan2_towardzero": double: 1 -float: 2 Function: "atan2_upward": double: 1 -float: 2 Function: "atan_downward": double: 1 @@ -166,19 +162,15 @@ double: 2 float: 2 Function: "carg": -float: 1 Function: "carg_downward": double: 1 -float: 2 Function: "carg_towardzero": double: 1 -float: 2 Function: "carg_upward": double: 1 -float: 2 Function: Real part of "casin": double: 1 diff --git a/sysdeps/hppa/fpu/libm-test-ulps b/sysdeps/hppa/fpu/libm-test-ulps index 74936184b5..5d2074b098 100644 --- a/sysdeps/hppa/fpu/libm-test-ulps +++ b/sysdeps/hppa/fpu/libm-test-ulps @@ -85,19 +85,15 @@ Function: "atan": double: 1 Function: "atan2": -float: 2 Function: "atan2_downward": double: 1 -float: 2 Function: "atan2_towardzero": double: 1 -float: 2 Function: "atan2_upward": double: 1 -float: 2 Function: "atan2pi": double: 1 @@ -233,19 +229,15 @@ double: 3 float: 2 Function: "carg": -float: 1 Function: "carg_downward": double: 1 -float: 2 Function: "carg_towardzero": double: 1 -float: 2 Function: "carg_upward": double: 1 -float: 2 Function: Real part of "casin": double: 1 diff --git a/sysdeps/i386/fpu/e_atan2f.S b/sysdeps/i386/fpu/e_atan2f.S deleted file mode 100644 index 606d4a6f2d..0000000000 --- a/sysdeps/i386/fpu/e_atan2f.S +++ /dev/null @@ -1,30 +0,0 @@ -/* - * Public domain. - */ - -#include <machine/asm.h> -#include <i386-math-asm.h> -#include <libm-alias-finite.h> - -RCSID("$NetBSD: e_atan2f.S,v 1.1 1995/05/08 23:35:10 jtc Exp $") - -DEFINE_FLT_MIN - -#ifdef PIC -# define MO(op) op##@GOTOFF(%ecx) -#else -# define MO(op) op -#endif - - .text -ENTRY(__ieee754_atan2f) -#ifdef PIC - LOAD_PIC_REG (cx) -#endif - flds 4(%esp) - flds 8(%esp) - fpatan - FLT_CHECK_FORCE_UFLOW_NARROW - ret -END (__ieee754_atan2f) -libm_alias_finite (__ieee754_atan2f, __atan2f) diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps index 1ff3b5cf3a..3c16d8887b 100644 --- a/sysdeps/i386/fpu/libm-test-ulps +++ b/sysdeps/i386/fpu/libm-test-ulps @@ -139,19 +139,16 @@ ldouble: 1 Function: "atan2_downward": double: 1 -float: 1 float128: 2 ldouble: 1 Function: "atan2_towardzero": double: 1 -float: 1 float128: 3 ldouble: 1 Function: "atan2_upward": double: 1 -float: 1 float128: 2 ldouble: 1 @@ -163,19 +160,19 @@ ldouble: 1 Function: "atan2pi_downward": double: 1 -float: 1 +float: 2 float128: 2 ldouble: 2 Function: "atan2pi_towardzero": double: 1 -float: 1 +float: 2 float128: 2 ldouble: 2 Function: "atan2pi_upward": double: 1 -float: 1 +float: 2 float128: 2 ldouble: 2 @@ -364,19 +361,16 @@ ldouble: 1 Function: "carg_downward": double: 1 -float: 1 float128: 2 ldouble: 1 Function: "carg_towardzero": double: 1 -float: 1 float128: 3 ldouble: 1 Function: "carg_upward": double: 1 -float: 1 float128: 2 ldouble: 1 diff --git a/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps b/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps index c296fefd0f..302d5c63c8 100644 --- a/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps +++ b/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps @@ -139,19 +139,16 @@ ldouble: 1 Function: "atan2_downward": double: 1 -float: 1 float128: 2 ldouble: 1 Function: "atan2_towardzero": double: 1 -float: 1 float128: 3 ldouble: 1 Function: "atan2_upward": double: 1 -float: 1 float128: 2 ldouble: 1 @@ -163,19 +160,19 @@ ldouble: 1 Function: "atan2pi_downward": double: 1 -float: 1 +float: 2 float128: 2 ldouble: 2 Function: "atan2pi_towardzero": double: 1 -float: 1 +float: 2 float128: 2 ldouble: 2 Function: "atan2pi_upward": double: 1 -float: 1 +float: 2 float128: 2 ldouble: 2 @@ -364,19 +361,16 @@ ldouble: 1 Function: "carg_downward": double: 1 -float: 1 float128: 2 ldouble: 1 Function: "carg_towardzero": double: 1 -float: 1 float128: 3 ldouble: 1 Function: "carg_upward": double: 1 -float: 1 float128: 2 ldouble: 1 diff --git a/sysdeps/ieee754/flt-32/e_atan2f.c b/sysdeps/ieee754/flt-32/e_atan2f.c index 75174062e8..6c7bebc24b 100644 --- a/sysdeps/ieee754/flt-32/e_atan2f.c +++ b/sysdeps/ieee754/flt-32/e_atan2f.c @@ -1,94 +1,274 @@ -/* e_atan2f.c -- float version of e_atan2.c. - */ +/* Correctly-rounded arctangent function of two binary32 values. -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ +Copyright (c) 2022-2024 Alexei Sibidanov and Paul Zimmermann. + +The original version of this file was copied from the CORE-MATH +project (file src/binary32/atan2/atan2f.c, revision 7835c5d). + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. + +*/ #include <math.h> -#include <math_private.h> +#include <stdint.h> #include <libm-alias-finite.h> +#include "math_config.h" -static const float -tiny = 1.0e-30, -zero = 0.0, -pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */ -pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */ -pi = 3.1415927410e+00, /* 0x40490fdb */ -pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */ +static inline double +muldd (double xh, double xl, double ch, double cl, double *l) +{ + double ahlh = ch * xl; + double alhh = cl * xh; + double ahhh = ch * xh; + double ahhl = fma (ch, xh, -ahhh); + ahhl += alhh + ahlh; + ch = ahhh + ahhl; + *l = (ahhh - ch) + ahhl; + return ch; +} -float -__ieee754_atan2f (float y, float x) +static double +polydd (double xh, double xl, int n, const double c[][2], double *l) { - float z; - int32_t k,m,hx,hy,ix,iy; + int i = n - 1; + double ch = c[i][0]; + double cl = c[i][1]; + while (--i >= 0) + { + ch = muldd (xh, xl, ch, cl, &cl); + double th = ch + c[i][0]; + double tl = (c[i][0] - th) + ch; + ch = th; + cl += tl + c[i][1]; + } + *l = cl; + return ch; +} - GET_FLOAT_WORD(hx,x); - ix = hx&0x7fffffff; - GET_FLOAT_WORD(hy,y); - iy = hy&0x7fffffff; - if((ix>0x7f800000)|| - (iy>0x7f800000)) /* x or y is NaN */ - return x+y; - if(hx==0x3f800000) return __atanf(y); /* x=1.0 */ - m = ((hy>>31)&1)|((hx>>30)&2); /* 2*sign(x)+sign(y) */ +/* for y/x tiny, use Taylor approximation z - z^3/3 where z=y/x */ +static float +cr_atan2f_tiny (float y, float x) +{ + double dy = y; + double dx = x; + double z = dy / dx; + double e = fma (-z, x, y); + /* z * x + e = y thus y/x = z + e/x */ + static const double c = -0x1.5555555555555p-2; /* -1/3 rounded to nearest */ + double zz = z * z; + double cz = c * z; + e = e / x + cz * zz; + uint64_t t = asuint |