Piotr Caban : msvcrt: Import _hypot implementation from musl.

[Alexandre Julliard]

Module: wine
Branch: master
Commit: e55d341c177befba78fc4d72e1959cf9b057d248
URL:    https://source.winehq.org/git/wine.git/?a=commit;h=e55d341c177befba78fc4d72e1959cf9b057d248

Author: Piotr Caban <piotr at codeweavers.com>
Date:   Thu Jun  3 16:28:58 2021 +0200

msvcrt: Import _hypot implementation from musl.

Signed-off-by: Piotr Caban <piotr at codeweavers.com>
Signed-off-by: Alexandre Julliard <julliard at winehq.org>

---

 dlls/msvcrt/math.c    | 58 +++++++++++++++++++++++++++++++++++++++++++++++++--
 dlls/msvcrt/unixlib.c |  9 --------
 dlls/msvcrt/unixlib.h |  1 -
 3 files changed, 56 insertions(+), 12 deletions(-)

diff --git a/dlls/msvcrt/math.c b/dlls/msvcrt/math.c
index edc5f6363d3..4d53f51e2f7 100644
--- a/dlls/msvcrt/math.c
+++ b/dlls/msvcrt/math.c
@@ -3198,13 +3198,67 @@ double CDECL _logb(double x)
     return __ilogb(x);
 }
 
+static void sq(double *hi, double *lo, double x)
+{
+    double xh, xl, xc;
+
+    xc = x * (0x1p27 + 1);
+    xh = x - xc + xc;
+    xl = x - xh;
+    *hi = x * x;
+    *lo = xh * xh - *hi + 2 * xh * xl + xl * xl;
+}
+
 /*********************************************************************
  *		_hypot (MSVCRT.@)
+ *
+ * Copied from musl: src/math/hypot.c
  */
 double CDECL _hypot(double x, double y)
 {
-  /* FIXME: errno handling */
-  return unix_funcs->hypot( x, y );
+    UINT64 ux = *(UINT64*)&x, uy = *(UINT64*)&y, ut;
+    double hx, lx, hy, ly, z;
+    int ex, ey;
+
+    /* arrange |x| >= |y| */
+    ux &= -1ULL >> 1;
+    uy &= -1ULL >> 1;
+    if (ux < uy) {
+        ut = ux;
+        ux = uy;
+        uy = ut;
+    }
+
+    /* special cases */
+    ex = ux >> 52;
+    ey = uy >> 52;
+    x = *(double*)&ux;
+    y = *(double*)&uy;
+    /* note: hypot(inf,nan) == inf */
+    if (ey == 0x7ff)
+        return y;
+    if (ex == 0x7ff || uy == 0)
+        return x;
+    /* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */
+    /* 64 difference is enough for ld80 double_t */
+    if (ex - ey > 64)
+        return x + y;
+
+    /* precise sqrt argument in nearest rounding mode without overflow */
+    /* xh*xh must not overflow and xl*xl must not underflow in sq */
+    z = 1;
+    if (ex > 0x3ff + 510) {
+        z = 0x1p700;
+        x *= 0x1p-700;
+        y *= 0x1p-700;
+    } else if (ey < 0x3ff - 450) {
+        z = 0x1p-700;
+        x *= 0x1p700;
+        y *= 0x1p700;
+    }
+    sq(&hx, &lx, x);
+    sq(&hy, &ly, y);
+    return z * sqrt(ly + lx + hy + hx);
 }
 
 /*********************************************************************
diff --git a/dlls/msvcrt/unixlib.c b/dlls/msvcrt/unixlib.c
index 742244bb287..dff51ee8fd4 100644
--- a/dlls/msvcrt/unixlib.c
+++ b/dlls/msvcrt/unixlib.c
@@ -94,14 +94,6 @@ static float CDECL unix_fmaf( float x, float y, float z )
 #endif
 }
 
-/*********************************************************************
- *      hypot
- */
-static double CDECL unix_hypot(double x, double y)
-{
-    return hypot( x, y );
-}
-
 /*********************************************************************
  *      hypotf
  */
@@ -265,7 +257,6 @@ static const struct unix_funcs funcs =
     unix_exp2,
     unix_exp2f,
     unix_fmaf,
-    unix_hypot,
     unix_hypotf,
     unix_lgamma,
     unix_lgammaf,
diff --git a/dlls/msvcrt/unixlib.h b/dlls/msvcrt/unixlib.h
index 9e44c49b691..2429bd0556d 100644
--- a/dlls/msvcrt/unixlib.h
+++ b/dlls/msvcrt/unixlib.h
@@ -28,7 +28,6 @@ struct unix_funcs
     double          (CDECL *exp2)(double x);
     float           (CDECL *exp2f)(float x);
     float           (CDECL *fmaf)(float x, float y, float z);
-    double          (CDECL *hypot)(double x, double y);
     float           (CDECL *hypotf)(float x, float y);
     double          (CDECL *lgamma)(double x);
     float           (CDECL *lgammaf)(float x);