|
14 | 14 | namespace boost { |
15 | 15 | namespace decimal { |
16 | 16 |
|
17 | | -/* |
18 | | -constexpr auto fmad32(decimal32 x, decimal32 y, decimal32 z) noexcept -> decimal32 |
19 | | -{ |
20 | | - // First calculate x * y without rounding |
21 | | - constexpr decimal32 zero {0, 0}; |
22 | | -
|
23 | | - const auto res {detail::check_non_finite(x, y)}; |
24 | | - if (res != zero) |
25 | | - { |
26 | | - return res; |
27 | | - } |
28 | | -
|
29 | | - auto sig_lhs {x.full_significand()}; |
30 | | - auto exp_lhs {x.biased_exponent()}; |
31 | | - detail::normalize(sig_lhs, exp_lhs); |
32 | | -
|
33 | | - auto sig_rhs {y.full_significand()}; |
34 | | - auto exp_rhs {y.biased_exponent()}; |
35 | | - detail::normalize(sig_rhs, exp_rhs); |
36 | | -
|
37 | | - auto mul_result {detail::mul_impl<detail::decimal32_components>(sig_lhs, exp_lhs, x.isneg(), sig_rhs, exp_rhs, y.isneg())}; |
38 | | - const decimal32 dec_result {mul_result.sig, mul_result.exp, mul_result.sign}; |
39 | | -
|
40 | | - const auto res_add {detail::check_non_finite(dec_result, z)}; |
41 | | - if (res_add != zero) |
42 | | - { |
43 | | - return res_add; |
44 | | - } |
45 | | -
|
46 | | - bool lhs_bigger {dec_result > z}; |
47 | | - if (dec_result.isneg() && z.isneg()) |
48 | | - { |
49 | | - lhs_bigger = !lhs_bigger; |
50 | | - } |
51 | | - bool abs_lhs_bigger {abs(dec_result) > abs(z)}; |
52 | | -
|
53 | | - // To avoid the rounding step we promote the constituent pieces to the next higher type |
54 | | - detail::decimal64_components promoted_mul_result {static_cast<std::uint64_t>(mul_result.sig), |
55 | | - mul_result.exp, mul_result.sign}; |
56 | | -
|
57 | | - detail::normalize<decimal64>(promoted_mul_result.sig, promoted_mul_result.exp); |
58 | | -
|
59 | | - auto sig_z {static_cast<std::uint64_t>(z.full_significand())}; |
60 | | - auto exp_z {z.biased_exponent()}; |
61 | | - detail::normalize<decimal64>(sig_z, exp_z); |
62 | | - detail::decimal64_components z_components {sig_z, exp_z, z.isneg()}; |
63 | | -
|
64 | | - if (!lhs_bigger) |
65 | | - { |
66 | | - detail::swap(promoted_mul_result, z_components); |
67 | | - abs_lhs_bigger = !abs_lhs_bigger; |
68 | | - } |
69 | | -
|
70 | | - detail::decimal64_components result {}; |
71 | | -
|
72 | | - if (!promoted_mul_result.sign && z_components.sign) |
73 | | - { |
74 | | - result = detail::d64_sub_impl<detail::decimal64_components>(promoted_mul_result.sig, promoted_mul_result.exp, promoted_mul_result.sign, |
75 | | - z_components.sig, z_components.exp, z_components.sign, |
76 | | - abs_lhs_bigger); |
77 | | - } |
78 | | - else |
79 | | - { |
80 | | - result = detail::d64_add_impl<detail::decimal64_components>(promoted_mul_result.sig, promoted_mul_result.exp, promoted_mul_result.sign, |
81 | | - z_components.sig, z_components.exp, z_components.sign); |
82 | | - } |
83 | | -
|
84 | | - return {result.sig, result.exp, result.sign}; |
85 | | -} |
86 | | -
|
87 | | -constexpr auto fmad64(decimal64 x, decimal64 y, decimal64 z) noexcept -> decimal64 |
88 | | -{ |
89 | | - // First calculate x * y without rounding |
90 | | - constexpr decimal64 zero {0, 0}; |
91 | | -
|
92 | | - const auto res {detail::check_non_finite(x, y)}; |
93 | | - if (res != zero) |
94 | | - { |
95 | | - return res; |
96 | | - } |
97 | | -
|
98 | | - auto sig_lhs {x.full_significand()}; |
99 | | - auto exp_lhs {x.biased_exponent()}; |
100 | | - detail::normalize<decimal64>(sig_lhs, exp_lhs); |
101 | | -
|
102 | | - auto sig_rhs {y.full_significand()}; |
103 | | - auto exp_rhs {y.biased_exponent()}; |
104 | | - detail::normalize<decimal64>(sig_rhs, exp_rhs); |
105 | | -
|
106 | | - auto mul_result {detail::d64_mul_impl<detail::decimal64_components>(sig_lhs, exp_lhs, x.isneg(), sig_rhs, exp_rhs, y.isneg())}; |
107 | | - const decimal64 dec_result {mul_result.sig, mul_result.exp, mul_result.sign}; |
108 | | -
|
109 | | - const auto res_add {detail::check_non_finite(dec_result, z)}; |
110 | | - if (res_add != zero) |
111 | | - { |
112 | | - return res_add; |
113 | | - } |
114 | | -
|
115 | | - bool lhs_bigger {dec_result > z}; |
116 | | - if (dec_result.isneg() && z.isneg()) |
117 | | - { |
118 | | - lhs_bigger = !lhs_bigger; |
119 | | - } |
120 | | - bool abs_lhs_bigger {abs(dec_result) > abs(z)}; |
121 | | -
|
122 | | - // To avoid the rounding step we promote the constituent pieces to the next higher type |
123 | | - detail::decimal128_components promoted_mul_result {static_cast<detail::uint128>(mul_result.sig), |
124 | | - mul_result.exp, mul_result.sign}; |
125 | | -
|
126 | | - detail::normalize<decimal128>(promoted_mul_result.sig, promoted_mul_result.exp); |
127 | | -
|
128 | | - auto sig_z {static_cast<detail::uint128>(z.full_significand())}; |
129 | | - auto exp_z {z.biased_exponent()}; |
130 | | - detail::normalize<decimal128>(sig_z, exp_z); |
131 | | - detail::decimal128_components z_components {sig_z, exp_z, z.isneg()}; |
132 | | -
|
133 | | - if (!lhs_bigger) |
134 | | - { |
135 | | - detail::swap(promoted_mul_result, z_components); |
136 | | - abs_lhs_bigger = !abs_lhs_bigger; |
137 | | - } |
138 | | -
|
139 | | - detail::decimal128_components result {}; |
140 | | -
|
141 | | - if (!promoted_mul_result.sign && z_components.sign) |
142 | | - { |
143 | | - result = detail::d128_sub_impl<detail::decimal128_components>( |
144 | | - promoted_mul_result.sig, promoted_mul_result.exp, promoted_mul_result.sign, |
145 | | - z_components.sig, z_components.exp, z_components.sign, |
146 | | - abs_lhs_bigger); |
147 | | - } |
148 | | - else |
149 | | - { |
150 | | - result = detail::d128_add_impl<detail::decimal128_components>( |
151 | | - promoted_mul_result.sig, promoted_mul_result.exp, promoted_mul_result.sign, |
152 | | - z_components.sig, z_components.exp, z_components.sign); |
153 | | - } |
154 | | -
|
155 | | - return {result.sig, result.exp, result.sign}; |
156 | | -} |
157 | | -
|
158 | | -constexpr auto fmad128(decimal128 x, decimal128 y, decimal128 z) noexcept -> decimal128 |
159 | | -{ |
160 | | - return x * y + z; |
161 | | -} |
162 | | -
|
163 | | -// TODO(mborland): promote to decimal64_fast instead of regular decimal64 once it is available |
164 | | -constexpr auto fmad32f(decimal32_fast x, decimal32_fast y, decimal32_fast z) noexcept -> decimal32_fast |
165 | | -{ |
166 | | - // First calculate x * y without rounding |
167 | | - constexpr decimal32_fast zero {0, 0}; |
168 | | -
|
169 | | - const auto res {detail::check_non_finite(x, y)}; |
170 | | - if (res != zero) |
171 | | - { |
172 | | - return res; |
173 | | - } |
174 | | -
|
175 | | - auto sig_lhs {x.full_significand()}; |
176 | | - auto exp_lhs {x.biased_exponent()}; |
177 | | - detail::normalize(sig_lhs, exp_lhs); |
178 | | -
|
179 | | - auto sig_rhs {y.full_significand()}; |
180 | | - auto exp_rhs {y.biased_exponent()}; |
181 | | - detail::normalize(sig_rhs, exp_rhs); |
182 | | -
|
183 | | - auto mul_result {detail::mul_impl<detail::decimal32_fast_components>(sig_lhs, exp_lhs, x.isneg(), sig_rhs, exp_rhs, y.isneg())}; |
184 | | - const decimal32_fast dec_result {mul_result.sig, mul_result.exp, mul_result.sign}; |
185 | | -
|
186 | | - const auto res_add {detail::check_non_finite(dec_result, z)}; |
187 | | - if (res_add != zero) |
188 | | - { |
189 | | - return res_add; |
190 | | - } |
191 | | -
|
192 | | - bool lhs_bigger {dec_result > z}; |
193 | | - if (dec_result.isneg() && z.isneg()) |
194 | | - { |
195 | | - lhs_bigger = !lhs_bigger; |
196 | | - } |
197 | | - bool abs_lhs_bigger {abs(dec_result) > abs(z)}; |
198 | | -
|
199 | | - // To avoid the rounding step we promote the constituent pieces to the next higher type |
200 | | - detail::decimal64_components promoted_mul_result {static_cast<std::uint64_t>(mul_result.sig), |
201 | | - mul_result.exp, mul_result.sign}; |
202 | | -
|
203 | | - detail::normalize<decimal64>(promoted_mul_result.sig, promoted_mul_result.exp); |
204 | | -
|
205 | | - auto sig_z {static_cast<std::uint64_t>(z.full_significand())}; |
206 | | - auto exp_z {z.biased_exponent()}; |
207 | | - detail::normalize<decimal64>(sig_z, exp_z); |
208 | | - detail::decimal64_components z_components {sig_z, exp_z, z.isneg()}; |
209 | | -
|
210 | | - if (!lhs_bigger) |
211 | | - { |
212 | | - detail::swap(promoted_mul_result, z_components); |
213 | | - abs_lhs_bigger = !abs_lhs_bigger; |
214 | | - } |
215 | | -
|
216 | | - detail::decimal64_components result {}; |
217 | | -
|
218 | | - if (!promoted_mul_result.sign && z_components.sign) |
219 | | - { |
220 | | - result = detail::d64_sub_impl<detail::decimal64_components>(promoted_mul_result.sig, promoted_mul_result.exp, promoted_mul_result.sign, |
221 | | - z_components.sig, z_components.exp, z_components.sign, |
222 | | - abs_lhs_bigger); |
223 | | - } |
224 | | - else |
225 | | - { |
226 | | - result = detail::d64_add_impl<detail::decimal64_components>(promoted_mul_result.sig, promoted_mul_result.exp, promoted_mul_result.sign, |
227 | | - z_components.sig, z_components.exp, z_components.sign); |
228 | | - } |
229 | | -
|
230 | | - return {result.sig, result.exp, result.sign}; |
231 | | -} |
232 | | -
|
233 | | -constexpr auto fmad64f(decimal64_fast x, decimal64_fast y, decimal64_fast z) noexcept -> decimal64_fast |
234 | | -{ |
235 | | - // First calculate x * y without rounding |
236 | | - constexpr decimal64_fast zero {0, 0}; |
237 | | -
|
238 | | - const auto res {detail::check_non_finite(x, y)}; |
239 | | - if (res != zero) |
240 | | - { |
241 | | - return res; |
242 | | - } |
243 | | -
|
244 | | - auto sig_lhs {x.full_significand()}; |
245 | | - auto exp_lhs {x.biased_exponent()}; |
246 | | - detail::normalize<decimal64>(sig_lhs, exp_lhs); |
247 | | -
|
248 | | - auto sig_rhs {y.full_significand()}; |
249 | | - auto exp_rhs {y.biased_exponent()}; |
250 | | - detail::normalize<decimal64>(sig_rhs, exp_rhs); |
251 | | -
|
252 | | - auto mul_result {detail::d64_mul_impl<detail::decimal64_fast_components>(sig_lhs, exp_lhs, x.isneg(), sig_rhs, exp_rhs, y.isneg())}; |
253 | | - const decimal64_fast dec_result {mul_result.sig, mul_result.exp, mul_result.sign}; |
254 | | -
|
255 | | - const auto res_add {detail::check_non_finite(dec_result, z)}; |
256 | | - if (res_add != zero) |
257 | | - { |
258 | | - return res_add; |
259 | | - } |
260 | | -
|
261 | | - bool lhs_bigger {dec_result > z}; |
262 | | - if (dec_result.isneg() && z.isneg()) |
263 | | - { |
264 | | - lhs_bigger = !lhs_bigger; |
265 | | - } |
266 | | - bool abs_lhs_bigger {abs(dec_result) > abs(z)}; |
267 | | -
|
268 | | - // To avoid the rounding step we promote the constituent pieces to the next higher type |
269 | | - detail::decimal128_components promoted_mul_result {static_cast<detail::uint128>(mul_result.sig), |
270 | | - mul_result.exp, mul_result.sign}; |
271 | | -
|
272 | | - detail::normalize<decimal128>(promoted_mul_result.sig, promoted_mul_result.exp); |
273 | | -
|
274 | | - auto sig_z {static_cast<detail::uint128>(z.full_significand())}; |
275 | | - auto exp_z {z.biased_exponent()}; |
276 | | - detail::normalize<decimal128>(sig_z, exp_z); |
277 | | - detail::decimal128_components z_components {sig_z, exp_z, z.isneg()}; |
278 | | -
|
279 | | - if (!lhs_bigger) |
280 | | - { |
281 | | - detail::swap(promoted_mul_result, z_components); |
282 | | - abs_lhs_bigger = !abs_lhs_bigger; |
283 | | - } |
284 | | -
|
285 | | - detail::decimal128_components result {}; |
286 | | -
|
287 | | - if (!promoted_mul_result.sign && z_components.sign) |
288 | | - { |
289 | | - result = detail::d128_sub_impl<detail::decimal128_components>( |
290 | | - promoted_mul_result.sig, promoted_mul_result.exp, promoted_mul_result.sign, |
291 | | - z_components.sig, z_components.exp, z_components.sign, |
292 | | - abs_lhs_bigger); |
293 | | - } |
294 | | - else |
295 | | - { |
296 | | - result = detail::d128_add_impl<detail::decimal128_components>( |
297 | | - promoted_mul_result.sig, promoted_mul_result.exp, promoted_mul_result.sign, |
298 | | - z_components.sig, z_components.exp, z_components.sign); |
299 | | - } |
300 | | -
|
301 | | - return {result.sig, result.exp, result.sign}; |
302 | | -} |
303 | | -
|
304 | | -constexpr auto fmad128f(decimal128_fast x, decimal128_fast y, decimal128_fast z) noexcept -> decimal128_fast |
305 | | -{ |
306 | | - return x * y + z; |
307 | | -} |
308 | | -*/ |
309 | | - |
310 | 17 | BOOST_DECIMAL_EXPORT constexpr auto fma(decimal32 x, decimal32 y, decimal32 z) noexcept -> decimal32 |
311 | 18 | { |
312 | 19 | return x * y + z; |
|
0 commit comments