]>
Commit | Line | Data |
---|---|---|
77a2a8b4 AZ |
1 | /* Pythagorean addition using floats |
2 | Copyright (C) 2011 Free Software Foundation, Inc. | |
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 | |
17 | License along with the GNU C Library; see the file COPYING.LIB. If not, | |
18 | write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, | |
19 | Boston, MA 02111-1307, USA. */ | |
20 | ||
21 | #include "math.h" | |
22 | #include "math_private.h" | |
23 | ||
24 | ||
25 | static const float two30 = 1.0737418e09; | |
26 | static const float two50 = 1.1259000e15; | |
27 | static const float two60 = 1.1529221e18; | |
28 | static const float two126 = 8.5070592e+37; | |
29 | static const float twoM50 = 8.8817842e-16; | |
30 | static const float twoM60 = 6.7762644e-21; | |
31 | static const float pdnum = 1.1754939e-38; | |
32 | ||
33 | ||
34 | /* __ieee754_hypotf(x,y) | |
35 | * | |
36 | * This a FP only version without any FP->INT conversion. | |
37 | * It is similar to default C version, making appropriates | |
38 | * overflow and underflows checks as well scaling when it | |
39 | * is needed. | |
40 | */ | |
41 | ||
42 | #ifdef _ARCH_PWR7 | |
43 | /* POWER7 isinf and isnan optimizations are fast. */ | |
44 | # define TEST_INF_NAN(x, y) \ | |
45 | if (isinff(x) || isinff(y)) \ | |
46 | return INFINITY; \ | |
47 | if (isnanf(x) || isnanf(y)) \ | |
48 | return NAN; | |
49 | # else | |
50 | /* For POWER6 and below isinf/isnan triggers LHS and PLT calls are | |
51 | * costly (especially for POWER6). */ | |
52 | # define GET_TWO_FLOAT_WORD(f1,f2,i1,i2) \ | |
53 | do { \ | |
54 | ieee_float_shape_type gf_u1; \ | |
55 | ieee_float_shape_type gf_u2; \ | |
56 | gf_u1.value = (f1); \ | |
57 | gf_u2.value = (f2); \ | |
58 | (i1) = gf_u1.word; \ | |
59 | (i2) = gf_u2.word; \ | |
60 | } while (0) | |
61 | ||
62 | # define TEST_INF_NAN(x, y) \ | |
63 | do { \ | |
64 | int32_t hx, hy; \ | |
65 | GET_TWO_FLOAT_WORD(x, y, hx, hy); \ | |
66 | if (hy > hx) { \ | |
67 | uint32_t ht = hx; hx = hy; hy = ht; \ | |
68 | } \ | |
69 | if (hx >= 0x7f800000) { \ | |
70 | if (hx == 0x7f800000 || hy == 0x7f800000) \ | |
71 | return INFINITY; \ | |
72 | return NAN; \ | |
73 | } \ | |
74 | } while (0) | |
75 | #endif | |
76 | ||
77 | ||
78 | float | |
79 | __ieee754_hypotf (float x, float y) | |
80 | { | |
81 | x = fabsf (x); | |
82 | y = fabsf (y); | |
83 | ||
84 | TEST_INF_NAN (x, y); | |
85 | ||
86 | if (y > x) | |
87 | { | |
88 | float t = y; | |
89 | y = x; | |
90 | x = t; | |
91 | } | |
92 | if (y == 0.0 || (x / y) > two30) | |
93 | { | |
94 | return x + y; | |
95 | } | |
96 | if (x > two50) | |
97 | { | |
98 | x *= twoM60; | |
99 | y *= twoM60; | |
edc121be | 100 | return __ieee754_sqrtf (x * x + y * y) / twoM60; |
77a2a8b4 AZ |
101 | } |
102 | if (y < twoM50) | |
103 | { | |
104 | if (y <= pdnum) | |
105 | { | |
106 | x *= two126; | |
107 | y *= two126; | |
edc121be | 108 | return __ieee754_sqrtf (x * x + y * y) / two126; |
77a2a8b4 AZ |
109 | } |
110 | else | |
111 | { | |
112 | x *= two60; | |
113 | y *= two60; | |
edc121be | 114 | return __ieee754_sqrtf (x * x + y * y) / two60; |
77a2a8b4 AZ |
115 | } |
116 | } | |
edc121be | 117 | return __ieee754_sqrtf (x * x + y * y); |
77a2a8b4 | 118 | } |