msvcrt: Import acosf from musl.

dlls/msvcrt/math.c

@@ -219,18 +219,63 @@ INT CDECL MSVCRT__isnanf( float num )
 
 /*********************************************************************
  *      MSVCRT_acosf (MSVCRT.@)
+ *
+ * Copied from musl: src/math/acosf.c
  */
+static float acosf_R(float z)
+{
+    static const float pS0 = 1.6666586697e-01,
+                 pS1 = -4.2743422091e-02,
+                 pS2 = -8.6563630030e-03,
+                 qS1 = -7.0662963390e-01;
+
+    float p, q;
+    p = z * (pS0 + z * (pS1 + z * pS2));
+    q = 1.0f + z * qS1;
+    return p / q;
+}
+
 float CDECL MSVCRT_acosf( float x )
 {
-  /* glibc implements acos() as the FPU equivalent of atan2(sqrt(1 - x ^ 2), x).
-   * asin() uses a similar construction. This is bad because as x gets nearer to
-   * 1 the error in the expression "1 - x^2" can get relatively large due to
-   * cancellation. The sqrt() makes things worse. A safer way to calculate
-   * acos() is to use atan2(sqrt((1 - x) * (1 + x)), x). */
-  float ret = atan2f(sqrtf((1 - x) * (1 + x)), x);
-  if (x < -1.0 || x > 1.0 || !finitef(x))
-    return math_error(_DOMAIN, "acosf", x, 0, ret);
-  return ret;
+    static const float pio2_hi = 1.5707962513e+00,
+                 pio2_lo = 7.5497894159e-08;
+
+    float z, w, s, c, df;
+    unsigned int hx, ix;
+
+    hx = *(unsigned int*)&x;
+    ix = hx & 0x7fffffff;
+    /* |x| >= 1 or nan */
+    if (ix >= 0x3f800000) {
+        if (ix == 0x3f800000) {
+            if (hx >> 31)
+                return 2 * pio2_hi + 7.5231638453e-37;
+            return 0;
+        }
+        if (MSVCRT__isnanf(x)) return x;
+        return math_error(_DOMAIN, "acosf", x, 0, 0 / (x - x));
+    }
+    /* |x| < 0.5 */
+    if (ix < 0x3f000000) {
+        if (ix <= 0x32800000) /* |x| < 2**-26 */
+            return pio2_hi + 7.5231638453e-37;
+        return pio2_hi - (x - (pio2_lo - x * acosf_R(x * x)));
+    }
+    /* x < -0.5 */
+    if (hx >> 31) {
+        z = (1 + x) * 0.5f;
+        s = sqrtf(z);
+        w = acosf_R(z) * s - pio2_lo;
+        return 2 * (pio2_hi - (s + w));
+    }
+    /* x > 0.5 */
+    z = (1 - x) * 0.5f;
+    s = sqrtf(z);
+    hx = *(unsigned int*)&s & 0xfffff000;
+    df = *(float*)&hx;
+    c = (z - df * df) / (s + df);
+    w = acosf_R(z) * s + c;
+    return 2 * (df + w);
 }
 
 /*********************************************************************