]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/ldbl-96/e_hypotl.c
1 /* Euclidean distance function. Long Double/Binary96 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 /* 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:
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
26 - Handle required underflow exception for subnormal results.
28 [1] https://arxiv.org/pdf/1904.09481.pdf */
31 #include <math_private.h>
32 #include <math-underflow.h>
33 #include <libm-alias-finite.h>
35 #define SCALE 0x8p-8257L
36 #define LARGE_VAL 0xb.504f333f9de6484p+8188L
37 #define TINY_VAL 0x8p-8194L
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 long double
43 kernel (long double ax
, long double ay
)
46 long double h
= sqrtl (ax
* ax
+ ay
* ay
);
49 long double delta
= h
- ay
;
50 t1
= ax
* (2.0L * delta
- ax
);
51 t2
= (delta
- 2.0L * (ax
- ay
)) * delta
;
55 long double delta
= h
- ax
;
56 t1
= 2.0L * delta
* (ax
- 2.0L * ay
);
57 t2
= (4.0L * delta
- ay
) * ay
+ delta
* delta
;
60 h
-= (t1
+ t2
) / (2.0L * h
);
65 __ieee754_hypotl (long double x
, long double y
)
67 if (!isfinite(x
) || !isfinite(y
))
69 if ((isinf (x
) || isinf (y
))
70 && !issignaling (x
) && !issignaling (y
))
78 long double ax
= x
< y
? y
: x
;
79 long double ay
= x
< y
? x
: y
;
81 /* If ax is huge, scale both inputs down. */
82 if (__glibc_unlikely (ax
> LARGE_VAL
))
84 if (__glibc_unlikely (ay
<= ax
* EPS
))
87 return kernel (ax
* SCALE
, ay
* SCALE
) / SCALE
;
90 /* If ay is tiny, scale both inputs up. */
91 if (__glibc_unlikely (ay
< TINY_VAL
))
93 if (__glibc_unlikely (ax
>= ay
/ EPS
))
96 ax
= kernel (ax
/ SCALE
, ay
/ SCALE
) * SCALE
;
97 math_check_force_underflow_nonneg (ax
);
101 /* Common case: ax is not huge and ay is not tiny. */
102 if (__glibc_unlikely (ay
<= ax
* EPS
))
105 return kernel (ax
, ay
);
107 libm_alias_finite (__ieee754_hypotl
, __hypotl
)