1 | /* Copyright (C) 1997-2024 Free Software Foundation, Inc. |
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 |
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. |
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 |
12 | Lesser General Public License for more details. |
13 | |
14 | You should have received a copy of the GNU Lesser General Public |
15 | License along with the GNU C Library. If not, see |
16 | <https://www.gnu.org/licenses/>. */ |
17 | |
18 | #include <math.h> |
19 | #include <math_private.h> |
20 | #include "mathimpl.h" |
21 | #include <libm-alias-finite.h> |
22 | |
23 | #ifndef SUFF |
24 | #define SUFF |
25 | #endif |
26 | #ifndef float_type |
27 | #define float_type double |
28 | #endif |
29 | |
30 | #define CONCATX(a,b) __CONCAT(a,b) |
31 | #define s(name) CONCATX(name,SUFF) |
32 | #define m81(func) __m81_u(s(func)) |
33 | |
34 | float_type |
35 | s(__ieee754_pow) (float_type x, float_type y) |
36 | { |
37 | float_type z; |
38 | float_type ax; |
39 | unsigned long x_cond, y_cond; |
40 | |
41 | y_cond = __m81_test (val: y); |
42 | if (y_cond & __M81_COND_ZERO) |
43 | return 1.0; |
44 | if (y_cond & __M81_COND_NAN) |
45 | return x == 1.0 ? x : x + y; |
46 | |
47 | x_cond = __m81_test (val: x); |
48 | if (x_cond & __M81_COND_NAN) |
49 | return x + y; |
50 | |
51 | if (y_cond & __M81_COND_INF) |
52 | { |
53 | ax = s(fabs) (x: x); |
54 | if (ax == 1.0) |
55 | return ax; |
56 | if (ax > 1.0) |
57 | return y_cond & __M81_COND_NEG ? 0 : y; |
58 | else |
59 | return y_cond & __M81_COND_NEG ? -y : 0; |
60 | } |
61 | |
62 | if (s(fabs) (x: y) == 1.0) |
63 | return y_cond & __M81_COND_NEG ? 1 / x : x; |
64 | |
65 | if (y == 2) |
66 | return x * x; |
67 | if (y == 0.5 && !(x_cond & __M81_COND_NEG)) |
68 | return m81(sqrt) (x: x); |
69 | |
70 | if (x == 10.0) |
71 | { |
72 | __asm ("ftentox%.x %1, %0" : "=f" (z) : "f" (y)); |
73 | return z; |
74 | } |
75 | if (x == 2.0) |
76 | { |
77 | __asm ("ftwotox%.x %1, %0" : "=f" (z) : "f" (y)); |
78 | return z; |
79 | } |
80 | |
81 | ax = s(fabs) (x: x); |
82 | if (x_cond & (__M81_COND_INF | __M81_COND_ZERO) || ax == 1.0) |
83 | { |
84 | z = ax; |
85 | if (y_cond & __M81_COND_NEG) |
86 | z = 1 / z; |
87 | if (x_cond & __M81_COND_NEG) |
88 | { |
89 | if (y != m81(__rint) (mathop_x: y)) |
90 | { |
91 | if (x == -1) |
92 | z = (z - z) / (z - z); |
93 | } |
94 | else |
95 | goto maybe_negate; |
96 | } |
97 | return z; |
98 | } |
99 | |
100 | if (x_cond & __M81_COND_NEG) |
101 | { |
102 | if (y == m81(__rint) (mathop_x: y)) |
103 | { |
104 | z = m81(__ieee754_exp) (mathop_x: y * m81(__ieee754_log) (mathop_x: -x)); |
105 | maybe_negate: |
106 | /* We always use the long double format, since y is already in |
107 | this format and rounding won't change the result. */ |
108 | { |
109 | int32_t exponent; |
110 | uint32_t i0, i1; |
111 | GET_LDOUBLE_WORDS (exponent, i0, i1, y); |
112 | exponent = (exponent & 0x7fff) - 0x3fff; |
113 | if (exponent <= 31 |
114 | ? i0 & (1 << (31 - exponent)) |
115 | : (exponent <= 63 |
116 | && i1 & (1 << (63 - exponent)))) |
117 | z = -z; |
118 | } |
119 | } |
120 | else |
121 | z = (y - y) / (y - y); |
122 | } |
123 | else |
124 | z = m81(__ieee754_exp) (mathop_x: y * m81(__ieee754_log) (mathop_x: x)); |
125 | return z; |
126 | } |
127 | libm_alias_finite (s(__ieee754_pow), s (__pow)) |
128 | |