patching file sysdeps/ieee754/dbl-64/e_hypot.c Hunk #1 FAILED at 37. Hunk #2 FAILED at 112. 2 out of 2 hunks FAILED -- saving rejects to file sysdeps/ieee754/dbl-64/e_hypot.c.rej Reject file sysdeps/ieee754/dbl-64/e_hypot.c.rej: --- sysdeps/ieee754/dbl-64/e_hypot.c +++ sysdeps/ieee754/dbl-64/e_hypot.c @@ -37,73 +37,40 @@ #include #include +#define FAST_FMINMAX 1 +//#undef __FP_FAST_FMA + +#define SCALE 0x1p-600 +#define LARGE_VAL 0x1p+511 +#define TINY_VAL 0x1p-459 +#define EPS 0x1p-54 + + static inline double handle_errno (double r) { + r = math_narrow_eval (r); if (isinf (r)) __set_errno (ERANGE); return r; } -/* sqrt (DBL_EPSILON / 2.0) */ -#define SQRT_EPS_DIV_2 0x1.6a09e667f3bcdp-27 -/* DBL_MIN / (sqrt (DBL_EPSILON / 2.0)) */ -#define DBL_MIN_THRESHOLD 0x1.6a09e667f3bcdp-996 -/* eps (double) * sqrt (DBL_MIN)) */ -#define SCALE 0x1p-563 -/* 1 / eps (sqrt (DBL_MIN) */ -#define INV_SCALE 0x1p+563 -/* sqrt (DBL_MAX) */ -#define SQRT_DBL_MAX 0x1.6a09e667f3bccp+511 -/* sqrt (DBL_MIN) */ -#define SQRT_DBL_MIN 0x1p-511 - -double -__hypot (double x, double y) +static inline double +kernel (double ax, double ay) { - if ((isinf (x) || isinf (y)) - && !issignaling (x) && !issignaling (y)) - return INFINITY; - if (isnan (x) || isnan (y)) - return x + y; - - double ax = fabs (x); - double ay = fabs (y); - if (ay > ax) - { - double tmp = ax; - ax = ay; - ay = tmp; - } - - /* Widely varying operands. The DBL_MIN_THRESHOLD check is used to avoid - a spurious underflow from the multiplication. */ - if (ax >= DBL_MIN_THRESHOLD && ay <= ax * SQRT_EPS_DIV_2) - return (ay == 0.0) - ? ax - : handle_errno (math_narrow_eval (ax + DBL_TRUE_MIN)); + double t1, t2; +#ifdef __FP_FAST_FMA + t1 = ay + ay; + t2 = ax - ay; - double scale = SCALE; - if (ax > SQRT_DBL_MAX) - { - ax *= scale; - ay *= scale; - scale = INV_SCALE; - } - else if (ay < SQRT_DBL_MIN) - { - ax /= scale; - ay /= scale; - } + if (t1 >= ax) + return sqrt (fma (t1, ax, t2 * t2)); else - scale = 1.0; - + return sqrt (fma (ax, ax, ay * ay)); +#else double h = sqrt (ax * ax + ay * ay); - double t1, t2; - if (h == 0.0) - return h; - else if (h <= 2.0 * ay) + if (h <= 2.0 * ay) { double delta = h - ay; t1 = ax * (2.0 * delta - ax); @@ -112,14 +79,57 @@ __hypot (double x, double y) else { double delta = h - ax; - t1 = 2.0 * delta * (ax - 2 * ay); + t1 = 2.0 * delta * (ax - 2.0 * ay); t2 = (4.0 * delta - ay) * ay + delta * delta; } h -= (t1 + t2) / (2.0 * h); - h = math_narrow_eval (h * scale); - math_check_force_underflow_nonneg (h); - return handle_errno (h); + return h; +#endif +} + + +double +__hypot (double x, double y) +{ + if (!isfinite (x) || !isfinite (y)) + { + if ((isinf (x) || isinf (y)) + && !issignaling_inline (x) && !issignaling_inline (y)) + return INFINITY; + return x + y; + } + + x = fabs (x); + y = fabs (y); + + double ax = FAST_FMINMAX ? fmax (x, y) : (x < y ? y : x); + double ay = FAST_FMINMAX ? fmin (x, y) : (x < y ? x : y); + + if (__glibc_unlikely (ax > LARGE_VAL)) + { + if (__glibc_unlikely (ay <= ax * EPS)) + return handle_errno (ax + ay); + + return handle_errno (kernel (ax * SCALE, ay * SCALE) / SCALE); + } + + if (__glibc_unlikely (ay < TINY_VAL)) + { + if (__glibc_unlikely (ax >= ay / EPS)) + return math_narrow_eval (ax + ay); + + ax = math_narrow_eval (kernel (ax / SCALE, ay / SCALE) * SCALE); + math_check_force_underflow_nonneg (ax); + return ax; + } + + /* Common case: ax is not huge and ay is not tiny. */ + if (__glibc_unlikely (ay <= ax * EPS)) + return math_narrow_eval (ax + ay); + + return math_narrow_eval (kernel (ax, ay)); } + strong_alias (__hypot, __ieee754_hypot) libm_alias_finite (__ieee754_hypot, __hypot) #if LIBM_SVID_COMPAT