| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | //====================================================================================================================== | ||
| 2 | /* | ||
| 3 | POLYFLOAT - Extended precision floating points | ||
| 4 | Copyright : POLYFLOAT Contributors & Maintainers | ||
| 5 | SPDX-License-Identifier: BSL-1.0 | ||
| 6 | */ | ||
| 7 | //====================================================================================================================== | ||
| 8 | #pragma once | ||
| 9 | #include <eve/eve.hpp> | ||
| 10 | #include <polyfloat/polyfloat.hpp> | ||
| 11 | #include <polyfloat/module/core/convert.hpp> | ||
| 12 | #include <polyfloat/module/core/if_else.hpp> | ||
| 13 | #include <polyfloat/module/core/oneminus.hpp> | ||
| 14 | #include <polyfloat/module/core/is_normalized.hpp> | ||
| 15 | #include <polyfloat/types/triple_real.hpp> | ||
| 16 | #include <polyfloat/types/double_real.hpp> | ||
| 17 | |||
| 18 | namespace plf | ||
| 19 | { | ||
| 20 | |||
| 21 | //==================================================================================================================== | ||
| 22 | //! @name Compound Assignment Operators | ||
| 23 | //! @related polyfloat | ||
| 24 | //! @{ | ||
| 25 | //==================================================================================================================== | ||
| 26 | |||
| 27 | /// Adds the polyflat value `other` to `self` and returns the new value of `self`. | ||
| 28 | template<concepts::polyfloat T1, concepts::polyfloat_like T2> | ||
| 29 | 8617945 | constexpr auto& operator+=(T1& self, T2 other) noexcept | |
| 30 | requires(dimension_v<T1> >= dimension_v<T2>) | ||
| 31 | { | ||
| 32 | 8617945 | T1 oth{plf::convert(other, eve::as<eve::element_type_t<T1>>())}; | |
| 33 | if constexpr (dimension_v<T1> == 2u) | ||
| 34 | { | ||
| 35 | 3211150 | auto [xhi, xlo] = self; | |
| 36 | 3211150 | auto [yhi, ylo] = oth; | |
| 37 | 6406834 | auto [a, b] = eve::two_add(xhi, yhi); | |
| 38 | 6406834 | auto [c, d] = eve::two_add(xlo, ylo); | |
| 39 | 6406834 | auto [aa, cc] = eve::two_add[eve::raw](a, c); | |
| 40 | 3211150 | b += d; | |
| 41 | 3211150 | b += cc; | |
| 42 | 6406834 | auto [e, f] = eve::two_add(aa, b); | |
| 43 | 3211150 | auto test = eve::is_not_finite(e + f); | |
| 44 | 3211150 | self = plf::if_else(test, T1(a), plf::_::from_pair(e, f)); | |
| 45 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 1605577 times.
|
3211150 | EVE_ASSERT(eve::all(plf::is_normalized(self)), "output of += is not normalized"); |
| 46 | 3211150 | return self; | |
| 47 | } | ||
| 48 | else if constexpr (dimension_v<T1> == 3u) | ||
| 49 | { | ||
| 50 | 5406795 | auto [xhi, xmd, xlo] = self; | |
| 51 | 5406795 | auto [yhi, ymd, ylo] = oth; | |
| 52 | 10779414 | auto [a1, b1] = eve::two_add(xhi, yhi); | |
| 53 | 10779414 | auto [c1, d1] = eve::two_add(xmd, ymd); | |
| 54 | 10779414 | auto [e1, f1] = eve::two_add(xlo, ylo); | |
| 55 | 10779414 | auto [a2, c2] = eve::two_add[raw](a1, c1); | |
| 56 | 5406795 | b1 += f1; | |
| 57 | 10779414 | auto [d2, e2] = eve::two_add(d1, e1); | |
| 58 | 10779414 | auto [a3, d3] = eve::two_add[raw](a2, d2); | |
| 59 | 10779414 | auto [b2, c3] = eve::two_add(b1, c2); | |
| 60 | 5406795 | c3 += e2; | |
| 61 | 10779414 | auto [c4, d4] = eve::two_add(c3, d3); | |
| 62 | 10779414 | auto [b3, c5] = eve::two_add(b2, c4); | |
| 63 | 10779414 | auto [a4, b4] = eve::two_add[raw](a3, b3); | |
| 64 | 5406795 | c4 += d4; | |
| 65 | 10779414 | auto [b5, c6] = eve::two_add[raw](b4, c5); | |
| 66 | 10779414 | auto [a5, b6] = eve::two_add[raw](a4, b5); | |
| 67 | 10779414 | auto [b7, c7] = eve::two_add[raw](b6, c6); | |
| 68 | 5406795 | auto test = eve::is_not_finite(a5 + b7 + c7); | |
| 69 | 5406795 | self = plf::if_else(test, T1(a1), plf::_::from_triple(a5, b7, c7)); | |
| 70 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 2703498 times.
|
5406795 | EVE_ASSERT(eve::all(plf::is_normalized(self)), "output of += is not normalized"); |
| 71 | 5406795 | return self; | |
| 72 | } | ||
| 73 | } | ||
| 74 | |||
| 75 | /// Substracts the polyfloat value `other` to `self` and returns the new value of `self`. | ||
| 76 | template<concepts::polyfloat T1, concepts::polyfloat_like T2> | ||
| 77 | 4922686 | constexpr auto& operator-=(T1& self, T2 other) noexcept | |
| 78 | requires(dimension_v<T1> >= dimension_v<T2>) | ||
| 79 | { | ||
| 80 | 4922686 | return self += -other; | |
| 81 | } | ||
| 82 | |||
| 83 | /// Multiply the polyfloat value `other` to `self` and returns the new value of `self`. | ||
| 84 | template<concepts::polyfloat T1, concepts::polyfloat_like T2> | ||
| 85 | 8450605 | constexpr auto& operator*=(T1& self, T2 other) noexcept | |
| 86 | requires(dimension_v<T1> >= dimension_v<T2>) | ||
| 87 | { | ||
| 88 | 8450605 | T1 oth{plf::convert(other, eve::as<eve::element_type_t<T1>>())}; | |
| 89 | if constexpr (dimension_v<T1> == 2u) | ||
| 90 | { | ||
| 91 | 3226682 | auto [xhi, xlo] = self; | |
| 92 | 3226682 | auto [yhi, ylo] = oth; | |
| 93 | 6437040 | auto [p00, e00] = eve::two_prod(xhi, yhi); | |
| 94 | 3226682 | auto p01 = xhi * ylo; | |
| 95 | 3226682 | auto p10 = xlo * yhi; | |
| 96 | 3226682 | p01 += p10; | |
| 97 | 3226682 | e00 += p01; | |
| 98 | 6437040 | auto [p, e] = eve::two_add[raw](p00, e00); | |
| 99 | 3226682 | auto test = eve::is_not_finite(p + e); | |
| 100 | 3226682 | self = plf::if_else(test, T1(p00), plf::_::from_pair(p, e)); | |
| 101 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 1613341 times.
|
3226682 | EVE_ASSERT(eve::all(plf::is_normalized(self)), "output of += is not normalized"); |
| 102 | 3226682 | return self; | |
| 103 | } | ||
| 104 | else if constexpr (dimension_v<T1> == 3u) | ||
| 105 | { | ||
| 106 | 5223923 | auto [xh, xm, xl] = self; | |
| 107 | 5223923 | auto [yh, ym, yl] = oth; | |
| 108 | 10419522 | auto [p00, e00] = eve::two_prod(xh, yh); | |
| 109 | 5223923 | auto p = p00; | |
| 110 | 10419522 | auto [p01, e01] = eve::two_prod(xh, ym); | |
| 111 | 10419522 | auto [p10, e10] = eve::two_prod(xm, yh); | |
| 112 | 5223923 | auto p02 = eve::mul(xh, yl); | |
| 113 | 5223923 | auto p11 = eve::mul(xm, ym); | |
| 114 | 5223923 | auto p20 = eve::mul(xl, yh); | |
| 115 | 10419522 | auto [a1, b1] = eve::two_add(p01, p10); | |
| 116 | 5223923 | p01 = a1; | |
| 117 | 5223923 | p10 = b1; | |
| 118 | 5223923 | e01 += e10; | |
| 119 | 5223923 | p02 += p20; | |
| 120 | 10419522 | auto [a2, b2] = eve::two_add(e00, p01); | |
| 121 | 5223923 | e00 = a2; | |
| 122 | 5223923 | p01 = b2; | |
| 123 | 5223923 | p02 += p11; | |
| 124 | 10419522 | auto [a3, b3] = eve::two_add[eve::raw](p00, e00); | |
| 125 | 5223923 | p00 = a3; | |
| 126 | 5223923 | e00 = b3; | |
| 127 | 5223923 | p01 += p10; | |
| 128 | 5223923 | e01 += p02; | |
| 129 | 5223923 | p01 += e01; | |
| 130 | 10419522 | auto [a4, b4] = eve::two_add(e00, p01); | |
| 131 | 5223923 | e00 = a4; | |
| 132 | 5223923 | p01 = b4; | |
| 133 | 10419522 | auto [a5, b5] = eve::two_add[eve::raw](p00, e00); | |
| 134 | 5223923 | p00 = a5; | |
| 135 | 5223923 | e00 = b5; | |
| 136 | 10419522 | auto [a6, b6] = eve::two_add[eve::raw](e00, p01); | |
| 137 | 5223923 | e00 = a6; | |
| 138 | 5223923 | p01 = b6; | |
| 139 | 10419522 | auto [a7, b7] = eve::two_add[eve::raw](p00, e00); | |
| 140 | 5223923 | p00 = a7; | |
| 141 | 5223923 | e00 = b7; | |
| 142 | 5223923 | auto test = eve::is_not_finite(p00 + e00 + p01); | |
| 143 | 5223923 | self = plf::if_else(test, T1(p), plf::_::from_triple(p00, e00, p01)); | |
| 144 |
1/2✗ Branch 2 not taken.
✓ Branch 3 taken 2611995 times.
|
5223923 | EVE_ASSERT(eve::all(plf::is_normalized(self)), "output of *= is not normalized"); |
| 145 | 5223923 | return self; | |
| 146 | } | ||
| 147 | } | ||
| 148 | } | ||
| 149 | |||
| 150 | #include <polyfloat/types/ops1.hpp> //define +, -, *, / * is used by plf::rec and /= | ||
| 151 | |||
| 152 | namespace plf | ||
| 153 | { | ||
| 154 | |||
| 155 | namespace _ | ||
| 156 | { | ||
| 157 | 2987193 | template<concepts::polyfloat_like T> auto rec(T a) noexcept | |
| 158 | { | ||
| 159 | if constexpr (dimension_v<T> == 1 && std::is_floating_point_v<T>) return eve::rec[pedantic](a); | ||
| 160 | else | ||
| 161 | { | ||
| 162 | 2987193 | a = eve::if_else(eve::is_denormal(hi(a)), T(hi(a) * eve::zero(as(hi(a)))), a); | |
| 163 | 4480791 | auto correct = [&a](auto x) { | |
| 164 | 1493598 | auto a0 = plf::hi(a); | |
| 165 |
16/16✓ Branch 2 taken 467755 times.
✓ Branch 5 taken 467755 times.
✓ Branch 12 taken 471429 times.
✓ Branch 15 taken 471429 times.
✓ Branch 22 taken 280856 times.
✓ Branch 25 taken 280856 times.
✓ Branch 32 taken 271068 times.
✓ Branch 35 taken 271068 times.
✓ Branch 42 taken 424 times.
✓ Branch 45 taken 424 times.
✓ Branch 52 taken 613 times.
✓ Branch 55 taken 613 times.
✓ Branch 62 taken 558 times.
✓ Branch 65 taken 558 times.
✓ Branch 72 taken 895 times.
✓ Branch 75 taken 895 times.
|
1493598 | x = plf::if_else(eve::is_infinite(a0), T(eve::copysign(eve::zero(eve::as(a0)), a0)), x); |
| 166 |
16/16✓ Branch 2 taken 467755 times.
✓ Branch 5 taken 467755 times.
✓ Branch 12 taken 471429 times.
✓ Branch 15 taken 471429 times.
✓ Branch 22 taken 280856 times.
✓ Branch 25 taken 280856 times.
✓ Branch 32 taken 271068 times.
✓ Branch 35 taken 271068 times.
✓ Branch 42 taken 424 times.
✓ Branch 45 taken 424 times.
✓ Branch 52 taken 613 times.
✓ Branch 55 taken 613 times.
✓ Branch 62 taken 558 times.
✓ Branch 65 taken 558 times.
✓ Branch 72 taken 895 times.
✓ Branch 75 taken 895 times.
|
1493598 | return plf::if_else(eve::is_eqz(a0), T(eve::copysign(eve::inf(eve::as(a0)), a0)), x); |
| 167 | }; | ||
| 168 | if constexpr (dimension_v<T> == 2) | ||
| 169 | { | ||
| 170 | 1523419 | auto a0 = plf::hi(a); | |
| 171 | 1523419 | auto x0 = eve::rec[pedantic](a0); | |
| 172 | 1523419 | auto x1 = x0 + x0 * plf::oneminus(a * x0); | |
| 173 | 1523419 | return correct(x1); | |
| 174 | } | ||
| 175 | else if constexpr (dimension_v<T> == 3) | ||
| 176 | { | ||
| 177 | 1463774 | auto a0 = plf::_::from_pair(plf::hi(a), plf::md(a)); | |
| 178 | 1463774 | auto x0 = plf::_::rec(a0); | |
| 179 | 1463774 | auto x1 = x0 + x0 * plf::oneminus(a * x0); | |
| 180 | 1463774 | return correct(x1); | |
| 181 | } | ||
| 182 | } | ||
| 183 | // ok aussi mais est-ce plus rapide ? | ||
| 184 | // auto a0 = plf::hi(a); | ||
| 185 | // auto correct = [&a0](auto x){ | ||
| 186 | // x = plf::if_else(eve::is_infinite(a0), T(eve::copysign(eve::zero(eve::as(a0)), a0)), x); | ||
| 187 | // return plf::if_else(eve::is_eqz(a0), T(eve::copysign(eve::inf(eve::as(a0)), a0)), x); | ||
| 188 | // }; | ||
| 189 | // auto x0 = eve::rec[pedantic](a0); | ||
| 190 | // auto x1 = x0 + x0* plf::oneminus(a*x0); | ||
| 191 | // if constexpr (dimension_v<T> == 2) return correct(x1); | ||
| 192 | // x1 = x1 + x1* plf::oneminus(a*x1); | ||
| 193 | // x1 = x1 + x1* plf::oneminus(a*x1); | ||
| 194 | // return correct(x1); | ||
| 195 | // } | ||
| 196 | |||
| 197 | // original | ||
| 198 | // else if constexpr (dimension_v<T> == 2) | ||
| 199 | // { | ||
| 200 | // auto [a0, b0] = a; | ||
| 201 | // auto x0 = eve::rec[pedantic](a0); | ||
| 202 | // auto x1 = x0 + x0 * (T(1) - a * x0); | ||
| 203 | // x1 = plf::if_else(eve::is_infinite(hi(a)), T(eve::copysign(eve::zero(eve::as(hi(a))), plf::hi(a))), x1); | ||
| 204 | // return if_else(eve::is_eqz(hi(a)), T(eve::copysign(eve::inf(eve::as(hi(a))), plf::hi(a))), x1); | ||
| 205 | // } | ||
| 206 | // else if constexpr (dimension_v<T> == 3) | ||
| 207 | // { | ||
| 208 | // auto [a0, b0, c0] = a; | ||
| 209 | // auto x0 = eve::rec[pedantic](a0); | ||
| 210 | // auto x1 = x0 + x0 * (T(1) - a * x0); | ||
| 211 | // auto x2 = x1 + x1 * (T(1) - a * x1); | ||
| 212 | // auto x3 = x2 + x2 * (T(1) - a * x2); | ||
| 213 | // x3 = plf::if_else(eve::is_infinite(hi(a)), T(eve::copysign(eve::zero(eve::as(hi(a))), plf::hi(a))), x3); | ||
| 214 | // return if_else(eve::is_eqz(hi(a)), T(eve::copysign(eve::inf(eve::as(hi(a))), plf::hi(a))), x3); | ||
| 215 | // } | ||
| 216 | } | ||
| 217 | } | ||
| 218 | |||
| 219 | /// Divide `self` by the value `other` and returns the new value of `self`. | ||
| 220 | 1391499 | template<concepts::polyfloat T1, concepts::polyfloat_like T2> constexpr auto& operator/=(T1& self, T2 other) noexcept | |
| 221 | { | ||
| 222 | 1391499 | T1 oth{plf::convert(other, eve::as<eve::element_type_t<T1>>())}; | |
| 223 | 1391499 | auto deno = eve::is_denormal(hi(oth)); | |
| 224 |
2/2✓ Branch 1 taken 695772 times.
✓ Branch 2 taken 4 times.
|
1391499 | if (eve::none(deno)) return self *= _::rec(oth); |
| 225 | else | ||
| 226 | { | ||
| 227 | 8 | auto a = self * _::rec(oth); | |
| 228 | 8 | auto b = T1(hi(self) / hi(oth)); | |
| 229 | 8 | auto test = deno && eve::is_infinite(hi(a)); | |
| 230 | 8 | return self = plf::if_else(test, b, a); | |
| 231 | } | ||
| 232 | EVE_ASSERT(eve::all(plf::is_normalized(self)), "output of /= is not normalized"); | ||
| 233 | return self; | ||
| 234 | } | ||
| 235 | } | ||
| 236 | |||
| 237 | #include <polyfloat/types/ops2.hpp> | ||
| 238 |