-/* e_acosf.c -- float version of e_acos.c.
- */
+/* Correctly-rounded arc-cosine function for binary32 value.
-/*
- * ====================================================
- * 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) 2023-2024 Alexei Sibidanov.
+The original version of this file was copied from the CORE-MATH
+project (file src/binary32/acos/acosf.c, revision 61d7bef).
+
+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 <errno.h>
#include <math.h>
#include <math_private.h>
#include <libm-alias-finite.h>
+#include "math_config.h"
-static const float
-one = 1.0000000000e+00, /* 0x3F800000 */
-pi = 3.1415925026e+00, /* 0x40490fda */
-pio2_hi = 1.5707962513e+00, /* 0x3fc90fda */
-pio2_lo = 7.5497894159e-08, /* 0x33a22168 */
-pS0 = 1.6666667163e-01, /* 0x3e2aaaab */
-pS1 = -3.2556581497e-01, /* 0xbea6b090 */
-pS2 = 2.0121252537e-01, /* 0x3e4e0aa8 */
-pS3 = -4.0055535734e-02, /* 0xbd241146 */
-pS4 = 7.9153501429e-04, /* 0x3a4f7f04 */
-pS5 = 3.4793309169e-05, /* 0x3811ef08 */
-qS1 = -2.4033949375e+00, /* 0xc019d139 */
-qS2 = 2.0209457874e+00, /* 0x4001572d */
-qS3 = -6.8828397989e-01, /* 0xbf303361 */
-qS4 = 7.7038154006e-02; /* 0x3d9dc62e */
+static __attribute__ ((noinline)) float
+as_special (float x)
+{
+ const float pih = 0x1.921fb6p+1;
+ const float pil = -0x1p-24f;
+ uint32_t t = asuint (x);
+ if (t == (0x7fu << 23))
+ return 0.0f; /* x=1 */
+ if (t == (0x17fu << 23))
+ return pih + pil; /* x=-1 */
+ uint32_t ax = t << 1;
+ if (ax > (0xffu << 24))
+ return x + x; /* nan */
+ return __math_invalidf (0.0);
+}
+
+static inline double
+poly12 (double z, const double *c)
+{
+ double z2 = z * z, z4 = z2 * z2;
+ double c0 = c[0] + z * c[1];
+ double c2 = c[2] + z * c[3];
+ double c4 = c[4] + z * c[5];
+ double c6 = c[6] + z * c[7];
+ double c8 = c[8] + z * c[9];
+ double c10 = c[10] + z * c[11];
+ c0 += c2 * z2;
+ c4 += c6 * z2;
+ c8 += z2 * c10;
+ c0 += z4 * (c4 + z4 * c8);
+ return c0;
+}
float
-__ieee754_acosf(float x)
+__ieee754_acosf (float x)
{
- float z,p,q,r,w,s,c,df;
- int32_t hx,ix;
- GET_FLOAT_WORD(hx,x);
- ix = hx&0x7fffffff;
- if(ix==0x3f800000) { /* |x|==1 */
- if(hx>0) return 0.0; /* acos(1) = 0 */
- else return pi+(float)2.0*pio2_lo; /* acos(-1)= pi */
- } else if(ix>0x3f800000) { /* |x| >= 1 */
- return (x-x)/(x-x); /* acos(|x|>1) is NaN */
- }
- if(ix<0x3f000000) { /* |x| < 0.5 */
- if(ix<=0x32800000) return pio2_hi+pio2_lo;/*if|x|<=2**-26*/
- z = x*x;
- p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
- q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
- r = p/q;
- return pio2_hi - (x - (pio2_lo-x*r));
- } else if (hx<0) { /* x < -0.5 */
- z = (one+x)*(float)0.5;
- p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
- q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
- s = sqrtf(z);
- r = p/q;
- w = r*s-pio2_lo;
- return pi - (float)2.0*(s+w);
- } else { /* x > 0.5 */
- int32_t idf;
- z = (one-x)*(float)0.5;
- s = sqrtf(z);
- df = s;
- GET_FLOAT_WORD(idf,df);
- SET_FLOAT_WORD(df,idf&0xfffff000);
- c = (z-df*df)/(s+df);
- p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
- q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
- r = p/q;
- w = r*s+c;
- return (float)2.0*(df+w);
- }
+ const double pi2 = 0x1.921fb54442d18p+0;
+ static const double o[] = { 0, 0x1.921fb54442d18p+1 };
+ double xs = x;
+ double r;
+ uint32_t t = asuint (x);
+ uint32_t ax = t << 1;
+ if (__glibc_unlikely (ax >= 0x7f<<24))
+ return as_special (x);
+ if (__glibc_likely (ax < 0x7ec2a1dcu)) /* |x| < 0x1.c2a1dcp-1 */
+ {
+ static const double b[] =
+ {
+ 0x1.fffffffd9ccb8p-1, 0x1.5555c94838007p-3, 0x1.32ded4b7c20fap-4,
+ 0x1.8566df703309ep-5, -0x1.980c959bec9a3p-6, 0x1.56fbb04998344p-1,
+ -0x1.403d8e4c49f52p+2, 0x1.b06c3e9f311eap+4, -0x1.9ea97c4e2c21fp+6,
+ 0x1.200b8261cc61bp+8, -0x1.2274c2799a5c7p+9, 0x1.a558a59cc19d3p+9,
+ -0x1.aca4b6a529ffp+9, 0x1.228744703f813p+9, -0x1.d7dbb0b322228p+7,
+ 0x1.5c2018c0c0105p+5
+ };
+ /* Avoid spurious underflow exception. */
+ if (__glibc_unlikely (ax <= 0x40000000u)) /* |x| < 2^-63 */
+ return (float) pi2;
+ double z = xs;
+ double z2 = z * z;
+ double z4 = z2 * z2;
+ double z8 = z4 * z4;
+ double z16 = z8 * z8;
+ r = z * ((((b[0] + z2 * b[1]) + z4 * (b[2] + z2 * b[3]))
+ + z8 * ((b[4] + z2 * b[5]) + z4 * (b[6] + z2 * b[7])))
+ + z16 * (((b[8] + z2 * b[9]) + z4 * (b[10] + z2 * b[11]))
+ + z8
+ * ((b[12] + z2 * b[13])+ z4 * (b[14] + z2 * b[15]))));
+ float ub = 0x1.921fb54574191p+0 - r;
+ float lb = 0x1.921fb543118ap+0 - r;
+ if (ub == lb)
+ return ub;
+ }
+ /* accurate path */
+ if (ax < (0x7eu << 24))
+ {
+ static const double c[] =
+ {
+ 0x1.555555555529cp-3, 0x1.333333337e0ddp-4, 0x1.6db6db3b4465ep-5,
+ 0x1.f1c72e13ac306p-6, 0x1.6e89cebe06bc4p-6, 0x1.1c6dcf5289094p-6,
+ 0x1.c6dbbcc7c6315p-7, 0x1.8f8dc2615e996p-7, 0x1.a5833b7bf15e8p-8,
+ 0x1.43f44ace1665cp-6, -0x1.0fb17df881c73p-6, 0x1.07520c026b2d6p-5
+ };
+ if (t == 0x328885a3u)
+ return 0x1.921fb6p+0f + 0x1p-25;
+ if (t == 0x39826222u)
+ return 0x1.920f6ap+0f + 0x1p-25;
+ double x2 = xs * xs;
+ r = (pi2 - xs) - (xs * x2) * poly12 (x2, c);
+ }
+ else
+ {
+ static const double c[] =
+ {
+ 0x1.6a09e667f3bcbp+0, 0x1.e2b7dddff2db9p-4, 0x1.b27247ab42dbcp-6,
+ 0x1.02995cc4e0744p-7, 0x1.5ffb0276ec8eap-9, 0x1.033885a928decp-10,
+ 0x1.911f2be23f8c7p-12, 0x1.4c3c55d2437fdp-13, 0x1.af477e1d7b461p-15,
+ 0x1.abd6bdff67dcbp-15, -0x1.1717e86d0fa28p-16, 0x1.6ff526de46023p-16
+ };
+ double bx = fabs (xs);
+ double z = 1.0 - bx;
+ double s = copysign (sqrt (z), xs);
+ r = o[t >> 31] + s * poly12 (z, c);
+ }
+ return r;
}
libm_alias_finite (__ieee754_acosf, __acosf)
projects / thirdparty / glibc.git / commitdiff
