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.
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.
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.
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/>. */
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
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
27 - Handle required underflow exception for subnormal results.
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)
35 [1] https://arxiv.org/pdf/1904.09481.pdf */
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"
48 #define SCALE 0x1p-600
49 #define LARGE_VAL 0x1p+511
50 #define TINY_VAL 0x1p-459
54 handle_errno (double r
)
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. */
64 kernel (double ax
, double ay
)
72 return sqrt (fma (t1
, ax
, t2
* t2
));
74 return sqrt (fma (ax
, ax
, ay
* ay
));
77 double h
= sqrt (ax
* ax
+ ay
* ay
);
80 double delta
= h
- ay
;
81 t1
= ax
* (2.0 * delta
- ax
);
82 t2
= (delta
- 2.0 * (ax
- ay
)) * delta
;
86 double delta
= h
- ax
;
87 t1
= 2.0 * delta
* (ax
- 2.0 * ay
);
88 t2
= (4.0 * delta
- ay
) * ay
+ delta
* delta
;
91 h
-= (t1
+ t2
) / (2.0 * h
);
97 __hypot (double x
, double y
)
99 if (!isfinite(x
) || !isfinite(y
))
101 if ((isinf (x
) || isinf (y
))
102 && !issignaling_inline (x
) && !issignaling_inline (y
))
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
);
113 /* If ax is huge, scale both inputs down. */
114 if (__glibc_unlikely (ax
> LARGE_VAL
))
116 if (__glibc_unlikely (ay
<= ax
* EPS
))
117 return handle_errno (math_narrow_eval (ax
+ ay
));
119 return handle_errno (math_narrow_eval (kernel (ax
* SCALE
, ay
* SCALE
)
123 /* If ay is tiny, scale both inputs up. */
124 if (__glibc_unlikely (ay
< TINY_VAL
))
126 if (__glibc_unlikely (ax
>= ay
/ EPS
))
127 return math_narrow_eval (ax
+ ay
);
129 ax
= math_narrow_eval (kernel (ax
/ SCALE
, ay
/ SCALE
) * SCALE
);
130 math_check_force_underflow_nonneg (ax
);
134 /* Common case: ax is not huge and ay is not tiny. */
135 if (__glibc_unlikely (ay
<= ax
* EPS
))
138 return kernel (ax
, ay
);
140 strong_alias (__hypot
, __ieee754_hypot
)
141 libm_alias_finite (__ieee754_hypot
, __hypot
)
143 versioned_symbol (libm
, __hypot
, hypot
, GLIBC_2_35
);
144 libm_alias_double_other (__hypot
, hypot
)
146 libm_alias_double (__hypot
, hypot
)