]>
Commit | Line | Data |
---|---|---|
77a2a8b4 | 1 | /* Pythagorean addition using doubles |
2b778ceb | 2 | Copyright (C) 2011-2021 Free Software Foundation, Inc. |
77a2a8b4 AZ |
3 | This file is part of the GNU C Library |
4 | Contributed by Adhemerval Zanella <azanella@br.ibm.com>, 2011 | |
5 | ||
6 | The GNU C Library is free software; you can redistribute it and/or | |
7 | modify it under the terms of the GNU Library General Public License as | |
8 | published by the Free Software Foundation; either version 2 of the | |
9 | License, or (at your option) any later version. | |
10 | ||
11 | The GNU C Library is distributed in the hope that it will be useful, | |
12 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
13 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
14 | Library General Public License for more details. | |
15 | ||
16 | You should have received a copy of the GNU Library General Public | |
59ba27a6 | 17 | License along with the GNU C Library; see the file COPYING.LIB. If |
5a82c748 | 18 | not, see <https://www.gnu.org/licenses/>. */ |
77a2a8b4 | 19 | |
1ed0291c RH |
20 | #include <math.h> |
21 | #include <math_private.h> | |
8f5b00d3 | 22 | #include <math-underflow.h> |
e054f494 | 23 | #include <stdint.h> |
220622dd | 24 | #include <libm-alias-finite.h> |
77a2a8b4 | 25 | |
77a2a8b4 AZ |
26 | /* __ieee754_hypot(x,y) |
27 | * | |
28 | * This a FP only version without any FP->INT conversion. | |
29 | * It is similar to default C version, making appropriates | |
30 | * overflow and underflows checks as well scaling when it | |
31 | * is needed. | |
32 | */ | |
33 | ||
77a2a8b4 AZ |
34 | double |
35 | __ieee754_hypot (double x, double y) | |
36 | { | |
69461d98 AZ |
37 | if ((isinf (x) || isinf (y)) |
38 | && !issignaling (x) && !issignaling (y)) | |
39 | return INFINITY; | |
40 | if (isnan (x) || isnan (y)) | |
41 | return x + y; | |
42 | ||
77a2a8b4 AZ |
43 | x = fabs (x); |
44 | y = fabs (y); | |
45 | ||
77a2a8b4 AZ |
46 | if (y > x) |
47 | { | |
48 | double t = x; | |
49 | x = y; | |
50 | y = t; | |
51 | } | |
16e616a7 AZ |
52 | if (y == 0.0) |
53 | return x; | |
69461d98 | 54 | |
16e616a7 AZ |
55 | /* if y is higher enough, y * 2^60 might overflow. The tests if |
56 | y >= 1.7976931348623157e+308/2^60 (two60factor) and uses the | |
57 | appropriate check to avoid the overflow exception generation. */ | |
69461d98 AZ |
58 | if (y <= 0x1.fffffffffffffp+963 && x > (y * 0x1p+60)) |
59 | return x + y; | |
60 | ||
61 | if (x > 0x1p+500) | |
77a2a8b4 | 62 | { |
69461d98 AZ |
63 | x *= 0x1p-600; |
64 | y *= 0x1p-600; | |
65 | return sqrt (x * x + y * y) / 0x1p-600; | |
77a2a8b4 | 66 | } |
69461d98 | 67 | if (y < 0x1p-500) |
77a2a8b4 | 68 | { |
69461d98 | 69 | if (y <= 0x0.fffffffffffffp-1022) |
77a2a8b4 | 70 | { |
69461d98 AZ |
71 | x *= 0x1p+1022; |
72 | y *= 0x1p+1022; | |
73 | double ret = sqrt (x * x + y * y) / 0x1p+1022; | |
f6987f5a JM |
74 | math_check_force_underflow_nonneg (ret); |
75 | return ret; | |
77a2a8b4 AZ |
76 | } |
77 | else | |
78 | { | |
69461d98 AZ |
79 | x *= 0x1p+600; |
80 | y *= 0x1p+600; | |
81 | return sqrt (x * x + y * y) / 0x1p+600; | |
77a2a8b4 AZ |
82 | } |
83 | } | |
f67a8147 | 84 | return sqrt (x * x + y * y); |
77a2a8b4 | 85 | } |
220622dd WD |
86 | #ifndef __ieee754_hypot |
87 | libm_alias_finite (__ieee754_hypot, __hypot) | |
88 | #endif |