1  /* Euclidean distance function.  Double/Binary64 version.
       2     Copyright (C) 2021-2023 Free Software Foundation, Inc.
       3     This file is part of the GNU C Library.
       4  
       5     The GNU C Library is free software; you can redistribute it and/or
       6     modify it under the terms of the GNU Lesser General Public
       7     License as published by the Free Software Foundation; either
       8     version 2.1 of the License, or (at your option) any later version.
       9  
      10     The GNU C Library is distributed in the hope that it will be useful,
      11     but WITHOUT ANY WARRANTY; without even the implied warranty of
      12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      13     Lesser General Public License for more details.
      14  
      15     You should have received a copy of the GNU Lesser General Public
      16     License along with the GNU C Library; if not, see
      17     <https://www.gnu.org/licenses/>.  */
      18  
      19  /* The implementation uses a correction based on 'An Improved Algorithm for
      20     hypot(a,b)' by Carlos F. Borges [1] usingthe MyHypot3 with the following
      21     changes:
      22  
      23     - Handle qNaN and sNaN.
      24     - Tune the 'widely varying operands' to avoid spurious underflow
      25       due the multiplication and fix the return value for upwards
      26       rounding mode.
      27     - Handle required underflow exception for subnormal results.
      28  
      29     The expected ULP is ~0.792 or ~0.948 if FMA is used.  For FMA, the
      30     correction is not used and the error of sqrt (x^2 + y^2) is below 1 ULP
      31     if x^2 + y^2 is computed with less than 0.707 ULP error.  If |x| >= |2y|,
      32     fma (x, x, y^2) has ~0.625 ULP.  If |x| < |2y|, fma (|2x|, |y|, (x - y)^2)
      33     has ~0.625 ULP.
      34  
      35     [1] https://arxiv.org/pdf/1904.09481.pdf  */
      36  
      37  #include <errno.h>
      38  #include <math.h>
      39  #include <math_private.h>
      40  #include <math-underflow.h>
      41  #include <math-narrow-eval.h>
      42  #include <math-use-builtins.h>
      43  #include <math-svid-compat.h>
      44  #include <libm-alias-finite.h>
      45  #include <libm-alias-double.h>
      46  #include "math_config.h"
      47  
      48  #define SCALE     0x1p-600
      49  #define LARGE_VAL 0x1p+511
      50  #define TINY_VAL  0x1p-459
      51  #define EPS       0x1p-54
      52  
      53  static inline double
      54  handle_errno (double r)
      55  {
      56    if (isinf (r))
      57      __set_errno (ERANGE);
      58    return r;
      59  }
      60  
      61  /* Hypot kernel. The inputs must be adjusted so that ax >= ay >= 0
      62     and squaring ax, ay and (ax - ay) does not overflow or underflow.  */
      63  static inline double
      64  kernel (double ax, double ay)
      65  {
      66    double t1, t2;
      67  #ifdef __FP_FAST_FMA
      68    t1 = ay + ay;
      69    t2 = ax - ay;
      70  
      71    if (t1 >= ax)
      72      return sqrt (fma (t1, ax, t2 * t2));
      73    else
      74      return sqrt (fma (ax, ax, ay * ay));
      75  
      76  #else
      77    double h = sqrt (ax * ax + ay * ay);
      78    if (h <= 2.0 * ay)
      79      {
      80        double delta = h - ay;
      81        t1 = ax * (2.0 * delta - ax);
      82        t2 = (delta - 2.0 * (ax - ay)) * delta;
      83      }
      84    else
      85      {
      86        double delta = h - ax;
      87        t1 = 2.0 * delta * (ax - 2.0 * ay);
      88        t2 = (4.0 * delta - ay) * ay + delta * delta;
      89      }
      90  
      91    h -= (t1 + t2) / (2.0 * h);
      92    return h;
      93  #endif
      94  }
      95  
      96  double
      97  __hypot (double x, double y)
      98  {
      99    if (!isfinite(x) || !isfinite(y))
     100      {
     101        if ((isinf (x) || isinf (y))
     102  	  && !issignaling_inline (x) && !issignaling_inline (y))
     103  	return INFINITY;
     104        return x + y;
     105      }
     106  
     107    x = fabs (x);
     108    y = fabs (y);
     109  
     110    double ax = USE_FMAX_BUILTIN ? fmax (x, y) : (x < y ? y : x);
     111    double ay = USE_FMIN_BUILTIN ? fmin (x, y) : (x < y ? x : y);
     112  
     113    /* If ax is huge, scale both inputs down.  */
     114    if (__glibc_unlikely (ax > LARGE_VAL))
     115      {
     116        if (__glibc_unlikely (ay <= ax * EPS))
     117  	return handle_errno (math_narrow_eval (ax + ay));
     118  
     119        return handle_errno (math_narrow_eval (kernel (ax * SCALE, ay * SCALE)
     120  					     / SCALE));
     121      }
     122  
     123    /* If ay is tiny, scale both inputs up.  */
     124    if (__glibc_unlikely (ay < TINY_VAL))
     125      {
     126        if (__glibc_unlikely (ax >= ay / EPS))
     127  	return math_narrow_eval (ax + ay);
     128  
     129        ax = math_narrow_eval (kernel (ax / SCALE, ay / SCALE) * SCALE);
     130        math_check_force_underflow_nonneg (ax);
     131        return ax;
     132      }
     133  
     134    /* Common case: ax is not huge and ay is not tiny.  */
     135    if (__glibc_unlikely (ay <= ax * EPS))
     136      return ax + ay;
     137  
     138    return kernel (ax, ay);
     139  }
     140  strong_alias (__hypot, __ieee754_hypot)
     141  libm_alias_finite (__ieee754_hypot, __hypot)
     142  #if LIBM_SVID_COMPAT
     143  versioned_symbol (libm, __hypot, hypot, GLIBC_2_35);
     144  libm_alias_double_other (__hypot, hypot)
     145  #else
     146  libm_alias_double (__hypot, hypot)
     147  #endif