diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index ba59069..7a56760 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -131,7 +131,9 @@ jobs:
 
       - name: GCOV
         if: runner.os == 'Linux'
-        run: gcovr --xml-pretty --exclude-unreachable-branches --exclude-throw-branches -r src -s build -o build/coverage.xml
+        run: |
+          ./build/src/cpp/pffdtd-engine test
+          gcovr --xml-pretty --exclude-unreachable-branches --exclude-throw-branches -r src -s build -o build/coverage.xml
 
       - name: Upload coverage report
         uses: codecov/codecov-action@v4
diff --git a/src/cpp/main.cpp b/src/cpp/main.cpp
index 53029b7..dc5a7a0 100644
--- a/src/cpp/main.cpp
+++ b/src/cpp/main.cpp
@@ -185,15 +185,34 @@ auto main(int argc, char** argv) -> int {
   if (*test) {
     // NOLINTBEGIN
     using pffdtd::Double;
+    PFFDTD_ASSERT(static_cast<float>(Double{42.0F} + 2.0F) == 44.0F);
     PFFDTD_ASSERT(static_cast<float>(Double{42.0F} + Double{2.0F}) == 44.0F);
     PFFDTD_ASSERT(static_cast<float>(Double{42.0F} - Double{2.0F}) == 40.0F);
     PFFDTD_ASSERT(static_cast<float>(Double{42.0F} * Double{2.0F}) == 84.0F);
     PFFDTD_ASSERT(static_cast<float>(Double{42.0F} / Double{2.0F}) == 21.0F);
 
+    PFFDTD_ASSERT(static_cast<double>(Double{42.0} + 2.0) == 44.0);
     PFFDTD_ASSERT(static_cast<double>(Double{42.0} + Double{2.0}) == 44.0);
     PFFDTD_ASSERT(static_cast<double>(Double{42.0} - Double{2.0}) == 40.0);
     PFFDTD_ASSERT(static_cast<double>(Double{42.0} * Double{2.0}) == 84.0);
     PFFDTD_ASSERT(static_cast<double>(Double{42.0} / Double{2.0}) == 21.0);
+
+    auto a = Double{42.0};
+    PFFDTD_ASSERT(a == a);
+    PFFDTD_ASSERT(a != Double<double>{});
+    PFFDTD_ASSERT(Double<double>{} != a);
+    PFFDTD_ASSERT(static_cast<double>(+a) == +42.0);
+    PFFDTD_ASSERT(static_cast<double>(-a) == -42.0);
+
+    a += Double{2.0};
+    PFFDTD_ASSERT(static_cast<double>(a) == 44.0);
+    a -= Double{2.0};
+    PFFDTD_ASSERT(static_cast<double>(a) == 42.0);
+    a *= Double{2.0};
+    PFFDTD_ASSERT(static_cast<double>(a) == 84.0);
+    a /= Double{2.0};
+    PFFDTD_ASSERT(static_cast<double>(a) == 42.0);
+
     // NOLINTEND
   }
 
diff --git a/src/cpp/pffdtd/double.hpp b/src/cpp/pffdtd/double.hpp
index e024618..6e668bd 100644
--- a/src/cpp/pffdtd/double.hpp
+++ b/src/cpp/pffdtd/double.hpp
@@ -17,6 +17,7 @@ namespace pffdtd {
 
 /// https://hal.science/hal-01351529v3/document
 /// https://inria.hal.science/inria-00070314/document
+/// https://people.eecs.berkeley.edu/~jrs/papers/robustr.pdf
 /// https://github.com/sukop/doubledouble
 /// https://github.com/JuliaMath/DoubleFloats.jl
 /// https://github.com/FlorisSteenkamp/double-double
@@ -53,29 +54,34 @@ struct Double {
 
   friend constexpr auto operator-(Double x) noexcept -> Double { return {-x.high(), -x.low()}; }
 
+  friend constexpr auto operator+(Double lhs, Real rhs) noexcept -> Double {
+    auto [sl, sh] = twoSum(lhs.high(), rhs);
+    return twoSumFast(sh, lhs.low() + sl);
+  }
+
   friend constexpr auto operator+(Double lhs, Double rhs) noexcept -> Double {
     auto [r, e] = twoSum(lhs.high(), rhs.high());
     e += lhs.low() + rhs.low();
-    return twoSumQuick(r, e);
+    return twoSumFast(r, e);
   }
 
   friend constexpr auto operator-(Double lhs, Double rhs) noexcept -> Double {
     auto [r, e] = twoDifference(lhs.high(), rhs.high());
     e += lhs.low() - rhs.low();
-    return twoSumQuick(r, e);
+    return twoSumFast(r, e);
   }
 
   friend constexpr auto operator*(Double lhs, Double rhs) noexcept -> Double {
     auto [r, e] = twoProduct(lhs.high(), rhs.high());
     e += lhs.high() * rhs.low() + lhs.low() * rhs.high();
-    return twoSumQuick(r, e);
+    return twoSumFast(r, e);
   }
 
   friend constexpr auto operator/(Double lhs, Double rhs) noexcept -> Double {
     auto r      = lhs.high() / rhs.high();
     auto [s, f] = twoProduct(r, rhs.high());
     auto e      = (lhs.high() - s - f + lhs.low() - r * rhs.low()) / rhs.high();
-    return twoSumQuick(r, e);
+    return twoSumFast(r, e);
   }
 
   friend constexpr auto operator+=(Double& lhs, Double rhs) noexcept -> Double& {
@@ -98,6 +104,8 @@ struct Double {
     return lhs;
   }
 
+  friend constexpr auto operator==(Double lhs, Double rhs) noexcept -> bool = default;
+
   private:
   [[nodiscard]] static constexpr auto split(Real a) noexcept -> Double {
     constexpr auto const digits = FloatTraits<Real>::digits;
@@ -121,7 +129,7 @@ struct Double {
     return Double{r, e};
   }
 
-  [[nodiscard]] static constexpr auto twoSumQuick(Real x, Real y) noexcept -> Double {
+  [[nodiscard]] static constexpr auto twoSumFast(Real x, Real y) noexcept -> Double {
     auto r = x + y;
     auto e = y - (r - x);
     return Double{r, e};
diff --git a/src/cpp/pffdtd/sycl.hpp b/src/cpp/pffdtd/sycl.hpp
index aecf749..da3c0e0 100644
--- a/src/cpp/pffdtd/sycl.hpp
+++ b/src/cpp/pffdtd/sycl.hpp
@@ -3,6 +3,7 @@
 
 #pragma once
 
+#include "pffdtd/float.hpp"
 #include "pffdtd/time.hpp"
 
 #if not defined(PFFDTD_HAS_SYCL)
@@ -27,4 +28,11 @@ template<typename Accessor>
 auto toString(sycl::info::device_type type) -> std::string;
 auto summary(sycl::device const& dev) -> void;
 
+template<>
+struct FloatTraits<sycl::half> {
+  static constexpr auto digits      = 11;
+  static constexpr auto minExponent = -13;
+  static constexpr auto maxExponent = 16;
+};
+
 } // namespace pffdtd