sysdeps/ieee754/ldbl-128/e_hypotl.c

   1 /* Euclidean distance function.  Long Double/Binary128 version.
   2    Copyright (C) 2021 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 /* This implementation is based on 'An Improved Algorithm for hypot(a,b)' by
  20    Carlos F. Borges [1] using the MyHypot3 with the following changes:
  21
  22    - Handle qNaN and sNaN.
  23    - Tune the 'widely varying operands' to avoid spurious underflow
  24      due the multiplication and fix the return value for upwards
  25      rounding mode.
  26    - Handle required underflow exception for subnormal results.
  27
  28    [1] https://arxiv.org/pdf/1904.09481.pdf  */
  29
  30 #include <math.h>
  31 #include <math_private.h>
  32 #include <math-underflow.h>
  33 #include <libm-alias-finite.h>
  34
  35 #define SCALE      L(0x1p-8303)
  36 #define LARGE_VAL  L(0x1.6a09e667f3bcc908b2fb1366ea95p+8191)
  37 #define TINY_VAL   L(0x1p-8191)
  38 #define EPS        L(0x1p-114)
  39
  40 /* Hypot kernel. The inputs must be adjusted so that ax >= ay >= 0
  41    and squaring ax, ay and (ax - ay) does not overflow or underflow.  */
  42 static inline _Float128
  43 kernel (_Float128 ax, _Float128 ay)
  44 {
  45   _Float128 t1, t2;
  46   _Float128 h = sqrtl (ax * ax + ay * ay);
  47   if (h <= L(2.0) * ay)
  48     {
  49       _Float128 delta = h - ay;
  50       t1 = ax * (L(2.0) * delta - ax);
  51       t2 = (delta - L(2.0) * (ax - ay)) * delta;
  52     }
  53   else
  54     {
  55       _Float128 delta = h - ax;
  56       t1 = L(2.0) * delta * (ax - L(2.0) * ay);
  57       t2 = (L(4.0) * delta - ay) * ay + delta * delta;
  58     }
  59
  60   h -= (t1 + t2) / (L(2.0) * h);
  61   return h;
  62 }
  63
  64 _Float128
  65 __ieee754_hypotl(_Float128 x, _Float128 y)
  66 {
  67   if (!isfinite(x) || !isfinite(y))
  68     {
  69       if ((isinf (x) || isinf (y))
  70           && !issignaling (x) && !issignaling (y))
  71         return INFINITY;
  72       return x + y;
  73     }
  74
  75   x = fabsl (x);
  76   y = fabsl (y);
  77
  78   _Float128 ax = x < y ? y : x;
  79   _Float128 ay = x < y ? x : y;
  80
  81   /* If ax is huge, scale both inputs down.  */
  82   if (__glibc_unlikely (ax > LARGE_VAL))
  83     {
  84       if (__glibc_unlikely (ay <= ax * EPS))
  85         return ax + ay;
  86
  87       return kernel (ax * SCALE, ay * SCALE) / SCALE;
  88     }
  89
  90   /* If ay is tiny, scale both inputs up.  */
  91   if (__glibc_unlikely (ay < TINY_VAL))
  92     {
  93       if (__glibc_unlikely (ax >= ay / EPS))
  94         return ax + ay;
  95
  96       ax = kernel (ax / SCALE, ay / SCALE) * SCALE;
  97       math_check_force_underflow_nonneg (ax);
  98       return ax;
  99     }
 100
 101   /* Common case: ax is not huge and ay is not tiny.  */
 102   if (__glibc_unlikely (ay <= ax * EPS))
 103     return ax + ay;
 104
 105   return kernel (ax, ay);
 106 }
 107 libm_alias_finite (__ieee754_hypotl, __hypotl)