diff --git a/tests/test_Array.cpp b/tests/test_Array.cpp
index 74dd2c37a96c84a7f72dc3d19138dae07d5ac6a2..7efc346e29693e2ee6302dfd4498ac2a9cb51b04 100644
--- a/tests/test_Array.cpp
+++ b/tests/test_Array.cpp
@@ -325,6 +325,14 @@ TEST_CASE("Array", "[utils]")
 
   SECTION("reproducible floating point sum")
   {
+    auto direct_sum = [](auto array) {
+      auto sum = array[0];
+      for (size_t i = 1; i < array.size(); ++i) {
+        sum += array[i];
+      }
+      return sum;
+    };
+
     SECTION("reproducible double sum")
     {
       Array<double> array(10'000);
@@ -340,6 +348,8 @@ TEST_CASE("Array", "[utils]")
       const auto sum_after_shuffle = sum(array);
 
       REQUIRE(sum_before_shuffle == sum_after_shuffle);
+
+      REQUIRE(sum_before_shuffle == Catch::Approx(direct_sum(array)));
     }
 
     SECTION("reproducible float sum")
@@ -357,6 +367,8 @@ TEST_CASE("Array", "[utils]")
       const auto sum_after_shuffle = sum(array);
 
       REQUIRE(sum_before_shuffle == sum_after_shuffle);
+
+      REQUIRE(sum_before_shuffle == Catch::Approx(direct_sum(array)));
     }
 
     SECTION("reproducible TinyVector<3,double> sum")
@@ -379,6 +391,11 @@ TEST_CASE("Array", "[utils]")
       const auto sum_after_shuffle = sum(array);
 
       REQUIRE(sum_before_shuffle == sum_after_shuffle);
+
+      auto naive_sum = direct_sum(array);
+      REQUIRE(sum_before_shuffle[0] == Catch::Approx(naive_sum[0]));
+      REQUIRE(sum_before_shuffle[1] == Catch::Approx(naive_sum[1]));
+      REQUIRE(sum_before_shuffle[2] == Catch::Approx(naive_sum[2]));
     }
 
     SECTION("reproducible TinyVector<3,float> sum")
@@ -401,6 +418,11 @@ TEST_CASE("Array", "[utils]")
       const auto sum_after_shuffle = sum(array);
 
       REQUIRE(sum_before_shuffle == sum_after_shuffle);
+
+      auto naive_sum = direct_sum(array);
+      REQUIRE(sum_before_shuffle[0] == Catch::Approx(naive_sum[0]));
+      REQUIRE(sum_before_shuffle[1] == Catch::Approx(naive_sum[1]));
+      REQUIRE(sum_before_shuffle[2] == Catch::Approx(naive_sum[2]));
     }
 
     SECTION("reproducible TinyMatrix<2, 3> sum")
@@ -424,6 +446,14 @@ TEST_CASE("Array", "[utils]")
       const auto sum_after_shuffle = sum(array);
 
       REQUIRE(sum_before_shuffle == sum_after_shuffle);
+
+      auto naive_sum = direct_sum(array);
+      REQUIRE(sum_before_shuffle(0, 0) == Catch::Approx(naive_sum(0, 0)));
+      REQUIRE(sum_before_shuffle(0, 1) == Catch::Approx(naive_sum(0, 1)));
+      REQUIRE(sum_before_shuffle(0, 2) == Catch::Approx(naive_sum(0, 2)));
+      REQUIRE(sum_before_shuffle(1, 0) == Catch::Approx(naive_sum(1, 0)));
+      REQUIRE(sum_before_shuffle(1, 1) == Catch::Approx(naive_sum(1, 1)));
+      REQUIRE(sum_before_shuffle(1, 2) == Catch::Approx(naive_sum(1, 2)));
     }
 
     SECTION("reproducible TinyMatrix<2, 3, float> sum")
@@ -447,6 +477,680 @@ TEST_CASE("Array", "[utils]")
       const auto sum_after_shuffle = sum(array);
 
       REQUIRE(sum_before_shuffle == sum_after_shuffle);
+
+      auto naive_sum = direct_sum(array);
+      REQUIRE(sum_before_shuffle(0, 0) == Catch::Approx(naive_sum(0, 0)));
+      REQUIRE(sum_before_shuffle(0, 1) == Catch::Approx(naive_sum(0, 1)));
+      REQUIRE(sum_before_shuffle(0, 2) == Catch::Approx(naive_sum(0, 2)));
+      REQUIRE(sum_before_shuffle(1, 0) == Catch::Approx(naive_sum(1, 0)).epsilon(1E-4));
+      REQUIRE(sum_before_shuffle(1, 1) == Catch::Approx(naive_sum(1, 1)).margin(1E-6));
+      REQUIRE(sum_before_shuffle(1, 2) == Catch::Approx(naive_sum(1, 2)));
+    }
+
+    SECTION("scalar bin summation")
+    {
+      SECTION("small arrays")
+      {
+        Array<double> array1(122);
+        for (size_t i = 0; i < array1.size(); ++i) {
+          array1[i] = ((i < (array1.size() / 10)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
+        }
+
+        Array<double> array2(357);
+        for (size_t i = 0; i < array2.size(); ++i) {
+          array2[i] = ((i < (array2.size() / 10)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
+        }
+
+        Array<double> array3(283);
+        for (size_t i = 0; i < array3.size(); ++i) {
+          array3[i] = ((i < (array3.size() / 10)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
+        }
+
+        Array<double> array(array1.size() + array2.size() + array3.size());
+        {
+          size_t I = 0;
+          for (size_t i = 0; i < array1.size(); ++i) {
+            array[I++] = array1[i];
+          }
+
+          for (size_t i = 0; i < array2.size(); ++i) {
+            array[I++] = array2[i];
+          }
+
+          for (size_t i = 0; i < array3.size(); ++i) {
+            array[I++] = array3[i];
+          }
+        }
+
+        ReproducibleScalarSum s1(array1);
+        ReproducibleScalarSum s2(array2);
+        ReproducibleScalarSum s3(array3);
+        using RSumT = decltype(s1);
+
+        RSumT::Bin bin1 = s1.getSummationBin();
+        RSumT::Bin bin2 = s2.getSummationBin();
+        RSumT::Bin bin3 = s3.getSummationBin();
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin1, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin3, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin3, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin1, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+      }
+
+      SECTION("small and large arrays")
+      {
+        Array<double> array1(122);
+        for (size_t i = 0; i < array1.size(); ++i) {
+          array1[i] = ((i < (array1.size() / 100)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
+        }
+
+        Array<double> array2(7357);
+        for (size_t i = 0; i < array2.size(); ++i) {
+          array2[i] = ((i < (array2.size() / 1000)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
+        }
+
+        Array<double> array3(283);
+        for (size_t i = 0; i < array3.size(); ++i) {
+          array3[i] = ((i < (array3.size() / 30)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
+        }
+
+        Array<double> array(array1.size() + array2.size() + array3.size());
+        {
+          size_t I = 0;
+          for (size_t i = 0; i < array1.size(); ++i) {
+            array[I++] = array1[i];
+          }
+
+          for (size_t i = 0; i < array2.size(); ++i) {
+            array[I++] = array2[i];
+          }
+
+          for (size_t i = 0; i < array3.size(); ++i) {
+            array[I++] = array3[i];
+          }
+        }
+
+        ReproducibleScalarSum s1(array1);
+        ReproducibleScalarSum s2(array2);
+        ReproducibleScalarSum s3(array3);
+        using RSumT = decltype(s1);
+
+        RSumT::Bin bin1 = s1.getSummationBin();
+        RSumT::Bin bin2 = s2.getSummationBin();
+        RSumT::Bin bin3 = s3.getSummationBin();
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin1, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin3, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin3, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin1, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+      }
+
+      SECTION("large arrays")
+      {
+        Array<double> array1(5122);
+        for (size_t i = 0; i < array1.size(); ++i) {
+          array1[i] = ((i < (array1.size() / 100)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
+        }
+
+        Array<double> array2(7357);
+        for (size_t i = 0; i < array2.size(); ++i) {
+          array2[i] = ((i < (array2.size() / 1000)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
+        }
+
+        Array<double> array3(6283);
+        for (size_t i = 0; i < array3.size(); ++i) {
+          array3[i] = ((i < (array3.size() / 300)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
+        }
+
+        Array<double> array(array1.size() + array2.size() + array3.size());
+        {
+          size_t I = 0;
+          for (size_t i = 0; i < array1.size(); ++i) {
+            array[I++] = array1[i];
+          }
+
+          for (size_t i = 0; i < array2.size(); ++i) {
+            array[I++] = array2[i];
+          }
+
+          for (size_t i = 0; i < array3.size(); ++i) {
+            array[I++] = array3[i];
+          }
+        }
+
+        ReproducibleScalarSum s1(array1);
+        ReproducibleScalarSum s2(array2);
+        ReproducibleScalarSum s3(array3);
+        using RSumT = decltype(s1);
+
+        RSumT::Bin bin1 = s1.getSummationBin();
+        RSumT::Bin bin2 = s2.getSummationBin();
+        RSumT::Bin bin3 = s3.getSummationBin();
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin1, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin3, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin3, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin1, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+      }
+    }
+
+    SECTION("TinyVector bin summation")
+    {
+      SECTION("small arrays")
+      {
+        Array<TinyVector<3, double>> array1(122);
+        for (size_t i = 0; i < array1.size(); ++i) {
+          array1[i] =
+            TinyVector<3, double>(((i < (array1.size() / 10)) * 1E25 + 1E20) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
+                                  ((i < (array1.size() / 50)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                  ((i < (array1.size() / 20)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1));
+        }
+
+        Array<TinyVector<3, double>> array2(357);
+        for (size_t i = 0; i < array2.size(); ++i) {
+          array2[i] =
+            TinyVector<3, double>(((i < (array2.size() / 10)) * 1E25 + 1E20) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
+                                  ((i < (array2.size() / 50)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                  ((i < (array2.size() / 20)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1));
+        }
+
+        Array<TinyVector<3, double>> array3(283);
+        for (size_t i = 0; i < array3.size(); ++i) {
+          array3[i] =
+            TinyVector<3, double>(((i < (array3.size() / 10)) * 1E25 + 1E20) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
+                                  ((i < (array3.size() / 50)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                  ((i < (array3.size() / 20)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1));
+        }
+
+        Array<TinyVector<3, double>> array(array1.size() + array2.size() + array3.size());
+        {
+          size_t I = 0;
+          for (size_t i = 0; i < array1.size(); ++i) {
+            array[I++] = array1[i];
+          }
+
+          for (size_t i = 0; i < array2.size(); ++i) {
+            array[I++] = array2[i];
+          }
+
+          for (size_t i = 0; i < array3.size(); ++i) {
+            array[I++] = array3[i];
+          }
+        }
+
+        ReproducibleTinyVectorSum s1(array1);
+        ReproducibleTinyVectorSum s2(array2);
+        ReproducibleTinyVectorSum s3(array3);
+        using RSumT = decltype(s1);
+
+        RSumT::Bin bin1 = s1.getSummationBin();
+        RSumT::Bin bin2 = s2.getSummationBin();
+        RSumT::Bin bin3 = s3.getSummationBin();
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin1, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin3, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin3, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin1, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+      }
+
+      SECTION("small and large arrays")
+      {
+        Array<TinyVector<3, double>> array1(122);
+        for (size_t i = 0; i < array1.size(); ++i) {
+          array1[i] =
+            TinyVector<3, double>(((i < (array1.size() / 10)) * 1E25 + 1E20) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
+                                  ((i < (array1.size() / 50)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                  ((i < (array1.size() / 20)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1));
+        }
+
+        Array<TinyVector<3, double>> array2(7357);
+        for (size_t i = 0; i < array2.size(); ++i) {
+          array2[i] =
+            TinyVector<3, double>(((i < (array2.size() / 10)) * 1E25 + 1E20) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
+                                  ((i < (array2.size() / 50)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                  ((i < (array2.size() / 20)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1));
+        }
+
+        Array<TinyVector<3, double>> array3(283);
+        for (size_t i = 0; i < array3.size(); ++i) {
+          array3[i] =
+            TinyVector<3, double>(((i < (array3.size() / 10)) * 1E25 + 1E20) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
+                                  ((i < (array3.size() / 50)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                  ((i < (array3.size() / 20)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1));
+        }
+
+        Array<TinyVector<3, double>> array(array1.size() + array2.size() + array3.size());
+        {
+          size_t I = 0;
+          for (size_t i = 0; i < array1.size(); ++i) {
+            array[I++] = array1[i];
+          }
+
+          for (size_t i = 0; i < array2.size(); ++i) {
+            array[I++] = array2[i];
+          }
+
+          for (size_t i = 0; i < array3.size(); ++i) {
+            array[I++] = array3[i];
+          }
+        }
+
+        ReproducibleTinyVectorSum s1(array1);
+        ReproducibleTinyVectorSum s2(array2);
+        ReproducibleTinyVectorSum s3(array3);
+        using RSumT = decltype(s1);
+
+        RSumT::Bin bin1 = s1.getSummationBin();
+        RSumT::Bin bin2 = s2.getSummationBin();
+        RSumT::Bin bin3 = s3.getSummationBin();
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin1, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin3, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin3, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin1, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+      }
+
+      SECTION("large arrays")
+      {
+        Array<TinyVector<3, double>> array1(5122);
+        for (size_t i = 0; i < array1.size(); ++i) {
+          array1[i] =
+            TinyVector<3, double>(((i < (array1.size() / 10)) * 1E25 + 1E20) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
+                                  ((i < (array1.size() / 50)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                  ((i < (array1.size() / 20)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1));
+        }
+
+        Array<TinyVector<3, double>> array2(7357);
+        for (size_t i = 0; i < array2.size(); ++i) {
+          array2[i] =
+            TinyVector<3, double>(((i < (array2.size() / 10)) * 1E25 + 1E20) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
+                                  ((i < (array2.size() / 50)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                  ((i < (array2.size() / 20)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1));
+        }
+
+        Array<TinyVector<3, double>> array3(4283);
+        for (size_t i = 0; i < array3.size(); ++i) {
+          array3[i] =
+            TinyVector<3, double>(((i < (array3.size() / 10)) * 1E25 + 1E20) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
+                                  ((i < (array3.size() / 50)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                  ((i < (array3.size() / 20)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1));
+        }
+
+        Array<TinyVector<3, double>> array(array1.size() + array2.size() + array3.size());
+        {
+          size_t I = 0;
+          for (size_t i = 0; i < array1.size(); ++i) {
+            array[I++] = array1[i];
+          }
+
+          for (size_t i = 0; i < array2.size(); ++i) {
+            array[I++] = array2[i];
+          }
+
+          for (size_t i = 0; i < array3.size(); ++i) {
+            array[I++] = array3[i];
+          }
+        }
+
+        ReproducibleTinyVectorSum s1(array1);
+        ReproducibleTinyVectorSum s2(array2);
+        ReproducibleTinyVectorSum s3(array3);
+        using RSumT = decltype(s1);
+
+        RSumT::Bin bin1 = s1.getSummationBin();
+        RSumT::Bin bin2 = s2.getSummationBin();
+        RSumT::Bin bin3 = s3.getSummationBin();
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin1, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin3, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin3, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin1, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+      }
+    }
+
+    SECTION("TinyMatrix bin summation")
+    {
+      SECTION("small arrays")
+      {
+        Array<TinyMatrix<2, 3, double>> array1(122);
+        for (size_t i = 0; i < array1.size(); ++i) {
+          array1[i] =
+            TinyMatrix<2, 3, double>(((i < (array1.size() / 10)) * 1E25 + 1E20) * ((i + 1) % 1'000) *
+                                       std::sin(3 * i + 1),
+                                     ((i < (array1.size() / 50)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                     ((i < (array1.size() / 20)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1),
+                                     ((i < (array1.size() / 30)) * 1E6 + 1E11) * ((i + 1) % 100) * std::cos(4 * i + 1),
+                                     ((i < (array1.size() / 70)) * 1E8 + 1E13) * ((i + 1) % 100) * std::sin(4 * i + 1),
+                                     ((i < (array1.size() / 11)) * 1E8 + 1E9) * ((i + 1) % 100) * std::cos(4 * i + 1));
+        }
+
+        Array<TinyMatrix<2, 3, double>> array2(357);
+        for (size_t i = 0; i < array2.size(); ++i) {
+          array2[i] =
+            TinyMatrix<2, 3, double>(((i < (array2.size() / 20)) * 1E25 + 1E20) * ((i + 1) % 1'000) *
+                                       std::sin(3 * i + 1),
+                                     ((i < (array2.size() / 70)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                     ((i < (array2.size() / 32)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1),
+                                     ((i < (array2.size() / 45)) * 1E6 + 1E11) * ((i + 1) % 100) * std::cos(4 * i + 1),
+                                     ((i < (array2.size() / 66)) * 1E8 + 1E13) * ((i + 1) % 100) * std::sin(4 * i + 1),
+                                     ((i < (array2.size() / 13)) * 1E8 + 1E9) * ((i + 1) % 100) * std::cos(4 * i + 1));
+        }
+
+        Array<TinyMatrix<2, 3, double>> array3(283);
+        for (size_t i = 0; i < array3.size(); ++i) {
+          array3[i] =
+            TinyMatrix<2, 3, double>(((i < (array3.size() / 23)) * 1E15 + 1E17) * ((i + 1) % 1'000) *
+                                       std::sin(3 * i + 1),
+                                     ((i < (array3.size() / 67)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                     ((i < (array3.size() / 62)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1),
+                                     ((i < (array3.size() / 35)) * 1E8 + 1E11) * ((i + 1) % 100) * std::cos(4 * i + 1),
+                                     ((i < (array3.size() / 63)) * 1E7 + 1E12) * ((i + 1) % 100) * std::sin(4 * i + 1),
+                                     ((i < (array3.size() / 17)) * 1E9 + 1E8) * ((i + 1) % 100) * std::cos(4 * i + 1));
+        }
+
+        Array<TinyMatrix<2, 3, double>> array(array1.size() + array2.size() + array3.size());
+        {
+          size_t I = 0;
+          for (size_t i = 0; i < array1.size(); ++i) {
+            array[I++] = array1[i];
+          }
+
+          for (size_t i = 0; i < array2.size(); ++i) {
+            array[I++] = array2[i];
+          }
+
+          for (size_t i = 0; i < array3.size(); ++i) {
+            array[I++] = array3[i];
+          }
+        }
+
+        ReproducibleTinyMatrixSum s1(array1);
+        ReproducibleTinyMatrixSum s2(array2);
+        ReproducibleTinyMatrixSum s3(array3);
+        using RSumT = decltype(s1);
+
+        RSumT::Bin bin1 = s1.getSummationBin();
+        RSumT::Bin bin2 = s2.getSummationBin();
+        RSumT::Bin bin3 = s3.getSummationBin();
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin1, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin3, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin3, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin1, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+      }
+
+      SECTION("small and large arrays")
+      {
+        Array<TinyMatrix<2, 3, double>> array1(122);
+        for (size_t i = 0; i < array1.size(); ++i) {
+          array1[i] =
+            TinyMatrix<2, 3, double>(((i < (array1.size() / 10)) * 1E25 + 1E20) * ((i + 1) % 1'000) *
+                                       std::sin(3 * i + 1),
+                                     ((i < (array1.size() / 50)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                     ((i < (array1.size() / 20)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1),
+                                     ((i < (array1.size() / 30)) * 1E6 + 1E11) * ((i + 1) % 100) * std::cos(4 * i + 1),
+                                     ((i < (array1.size() / 70)) * 1E8 + 1E13) * ((i + 1) % 100) * std::sin(4 * i + 1),
+                                     ((i < (array1.size() / 11)) * 1E8 + 1E9) * ((i + 1) % 100) * std::cos(4 * i + 1));
+        }
+
+        Array<TinyMatrix<2, 3, double>> array2(7357);
+        for (size_t i = 0; i < array2.size(); ++i) {
+          array2[i] =
+            TinyMatrix<2, 3, double>(((i < (array2.size() / 20)) * 1E25 + 1E20) * ((i + 1) % 1'000) *
+                                       std::sin(3 * i + 1),
+                                     ((i < (array2.size() / 70)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                     ((i < (array2.size() / 32)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1),
+                                     ((i < (array2.size() / 45)) * 1E6 + 1E11) * ((i + 1) % 100) * std::cos(4 * i + 1),
+                                     ((i < (array2.size() / 66)) * 1E8 + 1E13) * ((i + 1) % 100) * std::sin(4 * i + 1),
+                                     ((i < (array2.size() / 13)) * 1E8 + 1E9) * ((i + 1) % 100) * std::cos(4 * i + 1));
+        }
+
+        Array<TinyMatrix<2, 3, double>> array3(283);
+        for (size_t i = 0; i < array3.size(); ++i) {
+          array3[i] =
+            TinyMatrix<2, 3, double>(((i < (array3.size() / 23)) * 1E15 + 1E17) * ((i + 1) % 1'000) *
+                                       std::sin(3 * i + 1),
+                                     ((i < (array3.size() / 67)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                     ((i < (array3.size() / 62)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1),
+                                     ((i < (array3.size() / 35)) * 1E8 + 1E11) * ((i + 1) % 100) * std::cos(4 * i + 1),
+                                     ((i < (array3.size() / 63)) * 1E7 + 1E12) * ((i + 1) % 100) * std::sin(4 * i + 1),
+                                     ((i < (array3.size() / 17)) * 1E9 + 1E8) * ((i + 1) % 100) * std::cos(4 * i + 1));
+        }
+
+        Array<TinyMatrix<2, 3, double>> array(array1.size() + array2.size() + array3.size());
+        {
+          size_t I = 0;
+          for (size_t i = 0; i < array1.size(); ++i) {
+            array[I++] = array1[i];
+          }
+
+          for (size_t i = 0; i < array2.size(); ++i) {
+            array[I++] = array2[i];
+          }
+
+          for (size_t i = 0; i < array3.size(); ++i) {
+            array[I++] = array3[i];
+          }
+        }
+
+        ReproducibleTinyMatrixSum s1(array1);
+        ReproducibleTinyMatrixSum s2(array2);
+        ReproducibleTinyMatrixSum s3(array3);
+        using RSumT = decltype(s1);
+
+        RSumT::Bin bin1 = s1.getSummationBin();
+        RSumT::Bin bin2 = s2.getSummationBin();
+        RSumT::Bin bin3 = s3.getSummationBin();
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin1, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin3, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin3, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin1, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+      }
+
+      SECTION("large arrays")
+      {
+        Array<TinyMatrix<2, 3, double>> array1(3762);
+        for (size_t i = 0; i < array1.size(); ++i) {
+          array1[i] =
+            TinyMatrix<2, 3, double>(((i < (array1.size() / 10)) * 1E25 + 1E20) * ((i + 1) % 1'000) *
+                                       std::sin(3 * i + 1),
+                                     ((i < (array1.size() / 50)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                     ((i < (array1.size() / 20)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1),
+                                     ((i < (array1.size() / 30)) * 1E6 + 1E11) * ((i + 1) % 100) * std::cos(4 * i + 1),
+                                     ((i < (array1.size() / 70)) * 1E8 + 1E13) * ((i + 1) % 100) * std::sin(4 * i + 1),
+                                     ((i < (array1.size() / 11)) * 1E8 + 1E9) * ((i + 1) % 100) * std::cos(4 * i + 1));
+        }
+
+        Array<TinyMatrix<2, 3, double>> array2(7357);
+        for (size_t i = 0; i < array2.size(); ++i) {
+          array2[i] =
+            TinyMatrix<2, 3, double>(((i < (array2.size() / 20)) * 1E25 + 1E20) * ((i + 1) % 1'000) *
+                                       std::sin(3 * i + 1),
+                                     ((i < (array2.size() / 70)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                     ((i < (array2.size() / 32)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1),
+                                     ((i < (array2.size() / 45)) * 1E6 + 1E11) * ((i + 1) % 100) * std::cos(4 * i + 1),
+                                     ((i < (array2.size() / 66)) * 1E8 + 1E13) * ((i + 1) % 100) * std::sin(4 * i + 1),
+                                     ((i < (array2.size() / 13)) * 1E8 + 1E9) * ((i + 1) % 100) * std::cos(4 * i + 1));
+        }
+
+        Array<TinyMatrix<2, 3, double>> array3(6283);
+        for (size_t i = 0; i < array3.size(); ++i) {
+          array3[i] =
+            TinyMatrix<2, 3, double>(((i < (array3.size() / 23)) * 1E15 + 1E17) * ((i + 1) % 1'000) *
+                                       std::sin(3 * i + 1),
+                                     ((i < (array3.size() / 67)) * 1E15 + 1E10) * ((i + 1) % 100) * std::cos(3 * i + 1),
+                                     ((i < (array3.size() / 62)) * 1E7 + 1E10) * ((i + 1) % 10) * std::tan(3 * i + 1),
+                                     ((i < (array3.size() / 35)) * 1E8 + 1E11) * ((i + 1) % 100) * std::cos(4 * i + 1),
+                                     ((i < (array3.size() / 63)) * 1E7 + 1E12) * ((i + 1) % 100) * std::sin(4 * i + 1),
+                                     ((i < (array3.size() / 17)) * 1E9 + 1E8) * ((i + 1) % 100) * std::cos(4 * i + 1));
+        }
+
+        Array<TinyMatrix<2, 3, double>> array(array1.size() + array2.size() + array3.size());
+        {
+          size_t I = 0;
+          for (size_t i = 0; i < array1.size(); ++i) {
+            array[I++] = array1[i];
+          }
+
+          for (size_t i = 0; i < array2.size(); ++i) {
+            array[I++] = array2[i];
+          }
+
+          for (size_t i = 0; i < array3.size(); ++i) {
+            array[I++] = array3[i];
+          }
+        }
+
+        ReproducibleTinyMatrixSum s1(array1);
+        ReproducibleTinyMatrixSum s2(array2);
+        ReproducibleTinyMatrixSum s3(array3);
+        using RSumT = decltype(s1);
+
+        RSumT::Bin bin1 = s1.getSummationBin();
+        RSumT::Bin bin2 = s2.getSummationBin();
+        RSumT::Bin bin3 = s3.getSummationBin();
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin1, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin3, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+
+        {
+          RSumT::Bin bin_sum = zero;
+
+          RSumT::addBinTo(bin3, bin_sum);
+          RSumT::addBinTo(bin2, bin_sum);
+          RSumT::addBinTo(bin1, bin_sum);
+
+          REQUIRE(sum(array) == RSumT::getValue(bin_sum));
+        }
+      }
     }
   }