1/* Euclidean distance function. Double/Binary64 version.
2 Copyright (C) 2021-2022 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
53static inline double
54handle_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. */
63static inline double
64kernel (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
96double
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}
140strong_alias (__hypot, __ieee754_hypot)
141libm_alias_finite (__ieee754_hypot, __hypot)
142#if LIBM_SVID_COMPAT
143versioned_symbol (libm, __hypot, hypot, GLIBC_2_35);
144libm_alias_double_other (__hypot, hypot)
145#else
146libm_alias_double (__hypot, hypot)
147#endif
148