]>
Commit | Line | Data |
---|---|---|
d614a753 | 1 | /* Copyright (C) 1997-2020 Free Software Foundation, Inc. |
57a52ec8 UD |
2 | This file is part of the GNU C Library. |
3 | ||
4 | The GNU C Library is free software; you can redistribute it and/or | |
3214b89b AJ |
5 | modify it under the terms of the GNU Lesser General Public |
6 | License as published by the Free Software Foundation; either | |
7 | version 2.1 of the License, or (at your option) any later version. | |
57a52ec8 UD |
8 | |
9 | The GNU C Library is distributed in the hope that it will be useful, | |
10 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
11 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
3214b89b | 12 | Lesser General Public License for more details. |
57a52ec8 | 13 | |
3214b89b | 14 | You should have received a copy of the GNU Lesser General Public |
ab84e3ff | 15 | License along with the GNU C Library. If not, see |
5a82c748 | 16 | <https://www.gnu.org/licenses/>. */ |
57a52ec8 | 17 | |
57a52ec8 | 18 | #include <math.h> |
e038d690 | 19 | #include <math_private.h> |
5da9d124 | 20 | #include "mathimpl.h" |
57a52ec8 UD |
21 | |
22 | #ifndef SUFF | |
23 | #define SUFF | |
24 | #endif | |
25 | #ifndef float_type | |
26 | #define float_type double | |
27 | #endif | |
28 | ||
62f075cd UD |
29 | #define CONCATX(a,b) __CONCAT(a,b) |
30 | #define s(name) CONCATX(name,SUFF) | |
31 | #define m81(func) __m81_u(s(func)) | |
57a52ec8 UD |
32 | |
33 | float_type | |
62f075cd | 34 | s(__ieee754_pow) (float_type x, float_type y) |
57a52ec8 UD |
35 | { |
36 | float_type z; | |
37 | float_type ax; | |
4291e757 | 38 | unsigned long x_cond, y_cond; |
57a52ec8 | 39 | |
4291e757 UD |
40 | y_cond = __m81_test (y); |
41 | if (y_cond & __M81_COND_ZERO) | |
57a52ec8 | 42 | return 1.0; |
5d0b1535 AS |
43 | if (y_cond & __M81_COND_NAN) |
44 | return x == 1.0 ? x : x + y; | |
4291e757 UD |
45 | |
46 | x_cond = __m81_test (x); | |
5d0b1535 | 47 | if (x_cond & __M81_COND_NAN) |
57a52ec8 UD |
48 | return x + y; |
49 | ||
4291e757 | 50 | if (y_cond & __M81_COND_INF) |
57a52ec8 | 51 | { |
62f075cd | 52 | ax = s(fabs) (x); |
5d0b1535 AS |
53 | if (ax == 1.0) |
54 | return ax; | |
55 | if (ax > 1.0) | |
4291e757 | 56 | return y_cond & __M81_COND_NEG ? 0 : y; |
57a52ec8 | 57 | else |
4291e757 | 58 | return y_cond & __M81_COND_NEG ? -y : 0; |
57a52ec8 UD |
59 | } |
60 | ||
5d0b1535 | 61 | if (s(fabs) (y) == 1.0) |
4291e757 | 62 | return y_cond & __M81_COND_NEG ? 1 / x : x; |
57a52ec8 UD |
63 | |
64 | if (y == 2) | |
65 | return x * x; | |
4291e757 | 66 | if (y == 0.5 && !(x_cond & __M81_COND_NEG)) |
f67a8147 | 67 | return m81(sqrt) (x); |
57a52ec8 UD |
68 | |
69 | if (x == 10.0) | |
70 | { | |
71 | __asm ("ftentox%.x %1, %0" : "=f" (z) : "f" (y)); | |
72 | return z; | |
73 | } | |
74 | if (x == 2.0) | |
75 | { | |
76 | __asm ("ftwotox%.x %1, %0" : "=f" (z) : "f" (y)); | |
77 | return z; | |
78 | } | |
79 | ||
62f075cd | 80 | ax = s(fabs) (x); |
5d0b1535 | 81 | if (x_cond & (__M81_COND_INF | __M81_COND_ZERO) || ax == 1.0) |
57a52ec8 UD |
82 | { |
83 | z = ax; | |
4291e757 | 84 | if (y_cond & __M81_COND_NEG) |
57a52ec8 | 85 | z = 1 / z; |
4291e757 | 86 | if (x_cond & __M81_COND_NEG) |
57a52ec8 | 87 | { |
334ca657 | 88 | if (y != m81(__rint) (y)) |
57a52ec8 UD |
89 | { |
90 | if (x == -1) | |
4291e757 | 91 | z = (z - z) / (z - z); |
57a52ec8 UD |
92 | } |
93 | else | |
334ca657 | 94 | goto maybe_negate; |
57a52ec8 UD |
95 | } |
96 | return z; | |
97 | } | |
98 | ||
4291e757 | 99 | if (x_cond & __M81_COND_NEG) |
57a52ec8 | 100 | { |
334ca657 | 101 | if (y == m81(__rint) (y)) |
57a52ec8 | 102 | { |
62f075cd | 103 | z = m81(__ieee754_exp) (y * m81(__ieee754_log) (-x)); |
334ca657 UD |
104 | maybe_negate: |
105 | /* We always use the long double format, since y is already in | |
106 | this format and rounding won't change the result. */ | |
107 | { | |
108 | int32_t exponent; | |
24ab7723 | 109 | uint32_t i0, i1; |
334ca657 UD |
110 | GET_LDOUBLE_WORDS (exponent, i0, i1, y); |
111 | exponent = (exponent & 0x7fff) - 0x3fff; | |
112 | if (exponent <= 31 | |
113 | ? i0 & (1 << (31 - exponent)) | |
114 | : (exponent <= 63 | |
115 | && i1 & (1 << (63 - exponent)))) | |
116 | z = -z; | |
117 | } | |
57a52ec8 UD |
118 | } |
119 | else | |
4291e757 | 120 | z = (y - y) / (y - y); |
57a52ec8 UD |
121 | } |
122 | else | |
62f075cd | 123 | z = m81(__ieee754_exp) (y * m81(__ieee754_log) (x)); |
57a52ec8 UD |
124 | return z; |
125 | } | |
d3d9bde5 | 126 | strong_alias (s(__ieee754_pow), CONCATX (s(__pow), _finite)) |