1/* Euclidean distance function. Long Double/Binary128 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/* 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. */
42static inline _Float128
43kernel (_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}
107libm_alias_finite (__ieee754_hypotl, __hypotl)
108