| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | /* | ||
| 2 | POLYFLOAT - Extended precision floating points | ||
| 3 | Copyright : POLYFLOAT Contributors & Maintainers | ||
| 4 | SPDX-License-Identifier: BSL-1.0 | ||
| 5 | */ | ||
| 6 | //====================================================================================================================== | ||
| 7 | #pragma once | ||
| 8 | |||
| 9 | #include <polyfloat/details/callable.hpp> | ||
| 10 | #include <polyfloat/types/concepts.hpp> | ||
| 11 | #include <polyfloat/types/traits.hpp> | ||
| 12 | #include <type_traits> | ||
| 13 | #include <polyfloat/module/core/is_normalized.hpp> | ||
| 14 | |||
| 15 | /////////////////////////////////////////////////////////////////////////////////////// | ||
| 16 | // https://hal.science/hal-04575249 | ||
| 17 | /////////////////////////////////////////////////////////////////////////////////////// | ||
| 18 | |||
| 19 | namespace plf::_ | ||
| 20 | { | ||
| 21 | |||
| 22 | 174750 | template<typename T> auto inner_is_normalized(T h, T m, T l = T(0)) noexcept | |
| 23 | { | ||
| 24 | auto eqnum = eve::is_equal[eve::numeric]; | ||
| 25 | 174750 | return (eqnum(h + m, h) || eve::is_infinite(h + m)) && (eqnum(m + l, m) || eve::is_infinite(l + m)); | |
| 26 | } | ||
| 27 | |||
| 28 | 268482 | template<typename Z> constexpr auto inner_is_normalized(Z const& z) noexcept | |
| 29 | { | ||
| 30 | 268482 | auto eqnum = eve::is_equal[eve::numeric]; | |
| 31 | if constexpr (plf::dimension_v<Z> == 1) | ||
| 32 | { | ||
| 33 | 268466 | return eve::true_(eve::as(z)); | |
| 34 | } | ||
| 35 | else if constexpr (plf::dimension_v<Z> == 2) | ||
| 36 | { | ||
| 37 | 16 | auto [h, l] = z; | |
| 38 | 16 | return eqnum(h + l, h); | |
| 39 | } | ||
| 40 | else if constexpr (plf::dimension_v<Z> == 3) | ||
| 41 | { | ||
| 42 | auto [h, m, l] = z; | ||
| 43 | return eqnum(h + m, h) && eqnum(m + l, m); | ||
| 44 | } | ||
| 45 | } | ||
| 46 | |||
| 47 | 268462 | template<typename T> constexpr auto is_n_1or3_t_p2(T x) noexcept // Graillat & Muller algorithm 5 | |
| 48 | { | ||
| 49 | using u_t = eve::underlying_type_t<T>; | ||
| 50 | 268462 | constexpr unsigned long Q = 1ul << (eve::nbmantissabits(eve::as(u_t())) - 2); | |
| 51 | 268462 | constexpr unsigned long P = Q + 1; | |
| 52 | 268462 | auto l = P * x; | |
| 53 | 268462 | auto r = Q * x; | |
| 54 | 268462 | return eve::is_not_equal(l - r, x); | |
| 55 | } | ||
| 56 | |||
| 57 | template<typename T> | ||
| 58 | 268462 | constexpr auto inner_cr_dw_fp_add_with_err(T xh, T xl, T c) noexcept // Graillat & Muller algorithm 7 | |
| 59 | // https://hal.science/hal-04575249 | ||
| 60 | { | ||
| 61 | using u_t = eve::underlying_type_t<T>; | ||
| 62 | 268462 | constexpr u_t c9_8 = u_t(9) / 8; | |
| 63 | 268462 | constexpr u_t c7_8 = u_t(7) / 8; | |
| 64 | 536924 | auto [sh, sl] = eve::two_add(xh, c); | |
| 65 | 536924 | auto [vh, vl] = eve::two_add(xl, sl); | |
| 66 | 536924 | auto [wh, wl] = eve::two_add( | |
| 67 | vh, sh); // quick seems to be incorrectly called here in the original article as vh and sh can be in wrong order | ||
| 68 | 268462 | auto in = is_n_1or3_t_p2(vh) || eve::is_eqz(vl); | |
| 69 | 268462 | auto samsgn = eve::sign(vl) == eve::sign(vh); | |
| 70 | 268462 | auto z = eve::if_else(in, wh, sh + eve::if_else(samsgn, c9_8, c7_8) * vh); | |
| 71 | 268462 | auto delta = eve::sub[in](wl, z - wh); | |
| 72 | 268462 | return eve::zip(z, delta, vl); | |
| 73 | // return eve::if_else(eve::is_nan(z), eve::zip(z, u_t(0),u_t(0)), eve::zip(z, delta, vl)); | ||
| 74 | } | ||
| 75 | |||
| 76 | 268462 | template<typename Z> constexpr auto inner_three_add(Z const& x, Z const& y, Z const& z) noexcept | |
| 77 | { | ||
| 78 | 536924 | auto [xh, xl] = eve::two_add[eve::pedantic](x, y); | |
| 79 | 268462 | auto r = inner_cr_dw_fp_add_with_err(xh, xl, z); | |
| 80 | 268462 | auto isnf = eve::is_not_finite(x); | |
| 81 |
2/2✓ Branch 1 taken 521 times.
✓ Branch 2 taken 135096 times.
|
268462 | if (eve::any(isnf)) |
| 82 | { | ||
| 83 | 1039 | auto h = eve::zip(x, eve::zero(as(x)), eve::zero(as(x))); | |
| 84 | 1039 | r = eve::if_else(eve::is_not_finite(x), h, r); | |
| 85 | } | ||
| 86 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 135617 times.
|
268462 | EVE_ASSERT(eve::all(inner_is_normalized(r)), "[plf::_::inner_three_add] output is not normalized"); |
| 87 | 268462 | return r; | |
| 88 | } | ||
| 89 | } | ||
| 90 |