| 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_not_equal.hpp> | ||
| 14 | #include <polyfloat/module/core/is_eqz.hpp> | ||
| 15 | #include <polyfloat/module/core/sign.hpp> | ||
| 16 | #include <polyfloat/module/core/two_add.hpp> | ||
| 17 | #include <polyfloat/module/core/sub.hpp> | ||
| 18 | |||
| 19 | /////////////////////////////////////////////////////////////////////////////////////// | ||
| 20 | // https://hal.science/hal-04575249 | ||
| 21 | /////////////////////////////////////////////////////////////////////////////////////// | ||
| 22 | |||
| 23 | namespace plf::_ | ||
| 24 | { | ||
| 25 | 14001 | template<typename T> constexpr auto is_not_1or3_times_pow2(T x) noexcept // Graillat & Muller algorithm 5 | |
| 26 | { | ||
| 27 | using u_t = eve::underlying_type_t<T>; | ||
| 28 | 14001 | constexpr unsigned long Q = 1ul << (eve::nbmantissabits(eve::as(u_t())) - 2); | |
| 29 | 14001 | constexpr unsigned long P = Q + 1; | |
| 30 | 14001 | auto l = P * x; | |
| 31 | 14001 | auto r = Q * x; | |
| 32 | 14001 | return plf::is_not_equal(l - r, x); | |
| 33 | } | ||
| 34 | |||
| 35 | 20 | template<typename T> constexpr auto cr_dw_fp_add_with_err(T xh, T xl, T c) noexcept // Graillat & Muller algorithm 7 | |
| 36 | // https://hal.science/hal-04575249 | ||
| 37 | { | ||
| 38 | using u_t = eve::underlying_type_t<T>; | ||
| 39 | 20 | constexpr u_t c9_8 = u_t(9) / 8; | |
| 40 | 20 | constexpr u_t c7_8 = u_t(7) / 8; | |
| 41 | 40 | auto [sh, sl] = two_add(xh, c); | |
| 42 | 40 | auto [vh, vl] = two_add(xl, sl); | |
| 43 | 40 | auto [wh, wl] = two_add( | |
| 44 | vh, sh); // quick seems to be incorrectly called here in the original article as vh and sh can be in wrong order | ||
| 45 | 20 | auto in = is_not_1or3_times_pow2(vh) || is_eqz(vl); | |
| 46 | 20 | auto samsgn = sign(vl) == sign(vh); | |
| 47 | 20 | auto z = if_else(in, wh, sh + if_else(samsgn, c9_8, c7_8) * vh); | |
| 48 | 20 | auto delta = sub[in](wl, z - wh); | |
| 49 | 20 | return eve::zip(z, delta, vl); | |
| 50 | } | ||
| 51 | |||
| 52 | 13979 | template<typename T> constexpr auto cr_dw_fp_add(T xh, T xl, T c) noexcept // Graillat & Muller algorithm 6 | |
| 53 | { | ||
| 54 | using u_t = eve::underlying_type_t<T>; | ||
| 55 | 13979 | constexpr u_t c9_8 = u_t(9) / 8; | |
| 56 | 13979 | constexpr u_t c7_8 = u_t(7) / 8; | |
| 57 | 27430 | auto [sh, sl] = two_add(xh, c); | |
| 58 | 27430 | auto [vh, vl] = two_add(xl, sl); | |
| 59 | 13979 | auto in = is_not_1or3_times_pow2(vh) || is_eqz(vl); | |
| 60 | 13979 | auto samsgn = eve::sign(vl) == sign(vh); | |
| 61 | 13979 | return sh + if_else(in, vh, if_else(samsgn, c9_8, c7_8)); | |
| 62 | } | ||
| 63 | } | ||
| 64 |