Merge "Sync with upstream FreeBSD libm." into main

Author: enh-google

libm/upstream-freebsd/lib/msun/src/e_acosf.c

@@ -22,11 +22,17 @@ pi =  3.1415925026e+00, /* 0x40490fda */
 pio2_hi =  1.5707962513e+00; /* 0x3fc90fda */
 static volatile float
 pio2_lo =  7.5497894159e-08; /* 0x33a22168 */
+
+/*
+ * The coefficients for the rational approximation were generated over
+ *  0x1p-12f <= x <= 0.5f.  The maximum error satisfies log2(e) < -30.084.
+ */
 static const float
-pS0 =  1.6666586697e-01,
-pS1 = -4.2743422091e-02,
-pS2 = -8.6563630030e-03,
-qS1 = -7.0662963390e-01;
+pS0 =  1.66666672e-01f, /* 0x3e2aaaab */
+pS1 = -1.19510300e-01f, /* 0xbdf4c1d1 */
+pS2 =  5.47002675e-03f, /* 0x3bb33de9 */
+qS1 = -1.16706085e+00f, /* 0xbf956240 */
+qS2 =  2.90115148e-01f; /* 0x3e9489f9 */
 
 float
 acosf(float x)
@@ -46,13 +52,13 @@ acosf(float x)
 	    if(ix<=0x32800000) return pio2_hi+pio2_lo;/*if|x|<2**-26*/
 	    z = x*x;
 	    p = z*(pS0+z*(pS1+z*pS2));
-	    q = one+z*qS1;
+	    q = one+z*(qS1+z*qS2);
 	    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));
-	    q = one+z*qS1;
+	    q = one+z*(qS1+z*qS2);
 	    s = sqrtf(z);
 	    r = p/q;
 	    w = r*s-pio2_lo;
@@ -66,7 +72,7 @@ acosf(float x)
 	    SET_FLOAT_WORD(df,idf&0xfffff000);
 	    c  = (z-df*df)/(s+df);
 	    p = z*(pS0+z*(pS1+z*pS2));
-	    q = one+z*qS1;
+	    q = one+z*(qS1+z*qS2);
 	    r = p/q;
 	    w = r*s+c;
 	    return (float)2.0*(df+w);

libm/upstream-freebsd/lib/msun/src/e_asinf.c

@@ -18,12 +18,18 @@
 
 static const float
 one =  1.0000000000e+00, /* 0x3F800000 */
-huge =  1.000e+30,
-	/* coefficient for R(x^2) */
-pS0 =  1.6666586697e-01,
-pS1 = -4.2743422091e-02,
-pS2 = -8.6563630030e-03,
-qS1 = -7.0662963390e-01;
+huge =  1.000e+30;
+
+/*
+ * The coefficients for the rational approximation were generated over
+ *  0x1p-12f <= x <= 0.5f.  The maximum error satisfies log2(e) < -30.084.
+ */
+static const float
+pS0 =  1.66666672e-01f, /* 0x3e2aaaab */
+pS1 = -1.19510300e-01f, /* 0xbdf4c1d1 */
+pS2 =  5.47002675e-03f, /* 0x3bb33de9 */
+qS1 = -1.16706085e+00f, /* 0xbf956240 */
+qS2 =  2.90115148e-01f; /* 0x3e9489f9 */
 
 static const double
 pio2 =  1.570796326794896558e+00;
@@ -46,15 +52,15 @@ asinf(float x)
 	    }
 	    t = x*x;
 	    p = t*(pS0+t*(pS1+t*pS2));
-	    q = one+t*qS1;
+	    q = one+t*(qS1+t*qS2);
 	    w = p/q;
 	    return x+x*w;
 	}
 	/* 1> |x|>= 0.5 */
 	w = one-fabsf(x);
 	t = w*(float)0.5;
 	p = t*(pS0+t*(pS1+t*pS2));
-	q = one+t*qS1;
+	q = one+t*(qS1+t*qS2);
 	s = sqrt(t);
 	w = p/q;
 	t = pio2-2.0*(s+s*w);

libm/upstream-freebsd/lib/msun/src/math_private.h

@@ -405,7 +405,7 @@ do {								\
  * any extra precision into the type of 'a' -- 'a' should have type float_t,
  * double_t or long double.  b's type should be no larger than 'a's type.
  * Callers should use these types with scopes as large as possible, to
- * reduce their own extra-precision and efficiciency problems.  In
+ * reduce their own extra-precision and efficiency problems.  In
  * particular, they shouldn't convert back and forth just to call here.
  */
 #ifdef DEBUG

libm/upstream-freebsd/lib/msun/src/s_fma.c

@@ -260,14 +260,14 @@ fma(double x, double y, double z)
 
 	spread = ex + ey;
 
-	if (r.hi == 0.0) {
+	if (r.hi == 0.0 && xy.lo == 0) {
 		/*
 		 * When the addends cancel to 0, ensure that the result has
 		 * the correct sign.
 		 */
 		fesetround(oround);
 		volatile double vzs = zs; /* XXX gcc CSE bug workaround */
-		return (xy.hi + vzs + ldexp(xy.lo, spread));
+		return (xy.hi + vzs);
 	}
 
 	if (oround != FE_TONEAREST) {

libm/upstream-freebsd/lib/msun/src/s_fmal.c

@@ -241,14 +241,14 @@ fmal(long double x, long double y, long double z)
 
 	spread = ex + ey;
 
-	if (r.hi == 0.0) {
+	if (r.hi == 0.0 && xy.lo == 0) {
 		/*
 		 * When the addends cancel to 0, ensure that the result has
 		 * the correct sign.
 		 */
 		fesetround(oround);
 		volatile long double vzs = zs; /* XXX gcc CSE bug workaround */
-		return (xy.hi + vzs + ldexpl(xy.lo, spread));
+		return (xy.hi + vzs);
 	}
 
 	if (oround != FE_TONEAREST) {