#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>

#include <algebra/TinyMatrix.hpp>
#include <algebra/TinyVector.hpp>
#include <utils/Array.hpp>
#include <utils/PugsAssert.hpp>
#include <utils/Types.hpp>

#include <deque>
#include <list>
#include <random>
#include <set>
#include <unordered_set>
#include <valarray>
#include <vector>

#include <utils/ReproducibleSumUtils.hpp>
#include <utils/SmallArray.hpp>
#include <utils/Timer.hpp>

// Instantiate to ensure full coverage is performed
template class Array<int>;

// clazy:excludeall=non-pod-global-static

TEST_CASE("Array", "[utils]")
{
  Array<int> a(10);
  REQUIRE(a.size() == 10);

  for (size_t i = 0; i < a.size(); ++i) {
    a[i] = 2 * i;
  }

  REQUIRE(((a[0] == 0) and (a[1] == 2) and (a[2] == 4) and (a[3] == 6) and (a[4] == 8) and (a[5] == 10) and
           (a[6] == 12) and (a[7] == 14) and (a[8] == 16) and (a[9] == 18)));

  SECTION("checking for copies")
  {
    Array<const int> b{a};

    REQUIRE(((b[0] == 0) and (b[1] == 2) and (b[2] == 4) and (b[3] == 6) and (b[4] == 8) and (b[5] == 10) and
             (b[6] == 12) and (b[7] == 14) and (b[8] == 16) and (b[9] == 18)));

    Array<int> c{a};

    REQUIRE(((c[0] == 0) and (c[1] == 2) and (c[2] == 4) and (c[3] == 6) and (c[4] == 8) and (c[5] == 10) and
             (c[6] == 12) and (c[7] == 14) and (c[8] == 16) and (c[9] == 18)));

    Array<int> d = std::move(c);

    REQUIRE(((d[0] == 0) and (d[1] == 2) and (d[2] == 4) and (d[3] == 6) and (d[4] == 8) and (d[5] == 10) and
             (d[6] == 12) and (d[7] == 14) and (d[8] == 16) and (d[9] == 18)));
  }

  SECTION("checking for fill")
  {
    Array<int> b(10);
    b.fill(3);

    REQUIRE(((b[0] == 3) and (b[1] == 3) and (b[2] == 3) and (b[3] == 3) and (b[4] == 3) and (b[5] == 3) and
             (b[6] == 3) and (b[7] == 3) and (b[8] == 3) and (b[9] == 3)));
  }

  SECTION("checking for affectations (shallow copy)")
  {
    Array<const int> b;
    b = a;

    REQUIRE(((b[0] == 0) and (b[1] == 2) and (b[2] == 4) and (b[3] == 6) and (b[4] == 8) and (b[5] == 10) and
             (b[6] == 12) and (b[7] == 14) and (b[8] == 16) and (b[9] == 18)));

    Array<int> c;
    c = a;

    REQUIRE(((c[0] == 0) and (c[1] == 2) and (c[2] == 4) and (c[3] == 6) and (c[4] == 8) and (c[5] == 10) and
             (c[6] == 12) and (c[7] == 14) and (c[8] == 16) and (c[9] == 18)));

    Array<int> d;
    d = std::move(c);

    REQUIRE(((d[0] == 0) and (d[1] == 2) and (d[2] == 4) and (d[3] == 6) and (d[4] == 8) and (d[5] == 10) and
             (d[6] == 12) and (d[7] == 14) and (d[8] == 16) and (d[9] == 18)));
  }

  SECTION("checking for affectations (deep copy)")
  {
    Array<int> b(copy(a));

    REQUIRE(((b[0] == 0) and (b[1] == 2) and (b[2] == 4) and (b[3] == 6) and (b[4] == 8) and (b[5] == 10) and
             (b[6] == 12) and (b[7] == 14) and (b[8] == 16) and (b[9] == 18)));

    b.fill(2);

    REQUIRE(((a[0] == 0) and (a[1] == 2) and (a[2] == 4) and (a[3] == 6) and (a[4] == 8) and (a[5] == 10) and
             (a[6] == 12) and (a[7] == 14) and (a[8] == 16) and (a[9] == 18)));

    REQUIRE(((b[0] == 2) and (b[1] == 2) and (b[2] == 2) and (b[3] == 2) and (b[4] == 2) and (b[5] == 2) and
             (b[6] == 2) and (b[7] == 2) and (b[8] == 2) and (b[9] == 2)));

    Array<int> c;
    c = a;

    REQUIRE(((c[0] == 0) and (c[1] == 2) and (c[2] == 4) and (c[3] == 6) and (c[4] == 8) and (c[5] == 10) and
             (c[6] == 12) and (c[7] == 14) and (c[8] == 16) and (c[9] == 18)));

    c = copy(b);

    REQUIRE(((a[0] == 0) and (a[1] == 2) and (a[2] == 4) and (a[3] == 6) and (a[4] == 8) and (a[5] == 10) and
             (a[6] == 12) and (a[7] == 14) and (a[8] == 16) and (a[9] == 18)));

    REQUIRE(((c[0] == 2) and (c[1] == 2) and (c[2] == 2) and (c[3] == 2) and (c[4] == 2) and (c[5] == 2) and
             (c[6] == 2) and (c[7] == 2) and (c[8] == 2) and (c[9] == 2)));

    Array<int> d{a.size()};
    copy_to(a, d);

    REQUIRE(((d[0] == 0) and (d[1] == 2) and (d[2] == 4) and (d[3] == 6) and (d[4] == 8) and (d[5] == 10) and
             (d[6] == 12) and (d[7] == 14) and (d[8] == 16) and (d[9] == 18)));

    REQUIRE(((c[0] == 2) and (c[1] == 2) and (c[2] == 2) and (c[3] == 2) and (c[4] == 2) and (c[5] == 2) and
             (c[6] == 2) and (c[7] == 2) and (c[8] == 2) and (c[9] == 2)));

    copy_to(c, d);

    REQUIRE(((a[0] == 0) and (a[1] == 2) and (a[2] == 4) and (a[3] == 6) and (a[4] == 8) and (a[5] == 10) and
             (a[6] == 12) and (a[7] == 14) and (a[8] == 16) and (a[9] == 18)));

    REQUIRE(((d[0] == 2) and (d[1] == 2) and (d[2] == 2) and (d[3] == 2) and (d[4] == 2) and (d[5] == 2) and
             (d[6] == 2) and (d[7] == 2) and (d[8] == 2) and (d[9] == 2)));
  }

  SECTION("checking for std container conversion")
  {
    {
      std::vector<int> v{1, 2, 5, 3};
      {
        Array<int> v_array = convert_to_array(v);

        REQUIRE(v_array.size() == v.size());
        REQUIRE(((v_array[0] == 1) and (v_array[1] == 2) and (v_array[2] == 5) and (v_array[3] == 3)));
      }

      {
        Array<const int> v_array = convert_to_array(v);

        REQUIRE(v_array.size() == v.size());
        REQUIRE(((v_array[0] == 1) and (v_array[1] == 2) and (v_array[2] == 5) and (v_array[3] == 3)));
      }
    }

    {
      std::vector<int> w;
      {
        Array<int> w_array = convert_to_array(w);
        REQUIRE(w_array.size() == 0);
      }
      {
        Array<const int> w_array = convert_to_array(w);
        REQUIRE(w_array.size() == 0);
      }
    }

    {
      std::valarray<int> v{1, 2, 5, 3};
      Array<int> v_array = convert_to_array(v);

      REQUIRE(v_array.size() == v.size());
      REQUIRE(((v_array[0] == 1) and (v_array[1] == 2) and (v_array[2] == 5) and (v_array[3] == 3)));
    }

    {
      std::set<int> s{4, 2, 5, 3, 1, 3, 2};
      Array<int> s_array = convert_to_array(s);

      REQUIRE(s_array.size() == s.size());
      REQUIRE(
        ((s_array[0] == 1) and (s_array[1] == 2) and (s_array[2] == 3) and (s_array[3] == 4) and (s_array[4] == 5)));
    }

    {
      std::unordered_set<int> us{4, 2, 5, 3, 1, 3, 2};
      Array<int> us_array = convert_to_array(us);

      REQUIRE(us_array.size() == us.size());

      std::set<int> s;
      for (size_t i = 0; i < us_array.size(); ++i) {
        REQUIRE((us.find(us_array[i]) != us.end()));
        s.insert(us_array[i]);
      }
      REQUIRE(s.size() == us_array.size());
    }

    {
      std::multiset<int> ms{4, 2, 5, 3, 1, 3, 2};
      Array<int> ms_array = convert_to_array(ms);

      REQUIRE(ms_array.size() == ms.size());
      REQUIRE(((ms_array[0] == 1) and (ms_array[1] == 2) and (ms_array[2] == 2) and (ms_array[3] == 3) and
               (ms_array[4] == 3) and (ms_array[5] == 4) and (ms_array[6] == 5)));
    }

    {
      std::list<int> l{1, 3, 5, 6, 2};
      Array<int> l_array = convert_to_array(l);

      REQUIRE(l_array.size() == l.size());
      REQUIRE(
        ((l_array[0] == 1) and (l_array[1] == 3) and (l_array[2] == 5) and (l_array[3] == 6) and (l_array[4] == 2)));
    }

    {
      std::deque<int> q{1, 3, 5, 6, 2};
      q.push_front(2);
      Array<int> q_array = convert_to_array(q);

      REQUIRE(q_array.size() == q.size());
      REQUIRE(((q_array[0] == 2) and (q_array[1] == 1) and (q_array[2] == 3) and (q_array[3] == 5) and
               (q_array[4] == 6) and (q_array[5] == 2)));
    }
  }

  SECTION("checking for Kokkos::View encaspulation")
  {
    {
      Kokkos::View<double*> kokkos_view("anonymous", 10);
      for (size_t i = 0; i < kokkos_view.size(); ++i) {
        kokkos_view[i] = i;
      }

      Array array = encapsulate(kokkos_view);

      REQUIRE(array.size() == kokkos_view.size());
      for (size_t i = 0; i < array.size(); ++i) {
        REQUIRE(&array[i] == &kokkos_view[i]);
      }
    }
  }

  SECTION("output")
  {
    Array<int> x{5};
    x[0] = 2;
    x[1] = 6;
    x[2] = 2;
    x[3] = 3;
    x[4] = 7;

    std::ostringstream array_ost;
    array_ost << x;
    std::ostringstream ref_ost;
    ref_ost << 0 << ':' << x[0];
    for (size_t i = 1; i < x.size(); ++i) {
      ref_ost << ' ' << i << ':' << x[i];
    }
    REQUIRE(array_ost.str() == ref_ost.str());
  }

  SECTION("checking for Array reductions")
  {
    Array<int> b(10);
    b[0] = 13;
    b[1] = 1;
    b[2] = 8;
    b[3] = -3;
    b[4] = 23;
    b[5] = -1;
    b[6] = 13;
    b[7] = 0;
    b[8] = 12;
    b[9] = 9;

    SECTION("Min")
    {
      REQUIRE(min(b) == -3);
    }

    SECTION("Max")
    {
      REQUIRE(max(b) == 23);
    }

    SECTION("Sum")
    {
      REQUIRE((sum(b) == 75));
    }

    SECTION("TinyVector Sum")
    {
      using N2 = TinyVector<2, int>;
      Array<N2> c(10);
      c[0] = N2{13, 2};
      c[1] = N2{1, 3};
      c[2] = N2{8, -2};
      c[3] = N2{-3, 2};
      c[4] = N2{23, 4};
      c[5] = N2{-1, -3};
      c[6] = N2{13, 17};
      c[7] = N2{0, 9};
      c[8] = N2{12, 13};
      c[9] = N2{9, -17};

      REQUIRE((sum(c) == N2{75, 28}));
    }

    SECTION("TinyMatrix Sum")
    {
      using N22 = TinyMatrix<2, 2, int>;
      Array<N22> d(10);
      d[0] = N22{13, 2, 0, 1};
      d[1] = N22{1, 3, 6, 3};
      d[2] = N22{8, -2, -1, 21};
      d[3] = N22{-3, 2, 5, 12};
      d[4] = N22{23, 4, 7, 1};
      d[5] = N22{-1, -3, 33, 11};
      d[6] = N22{13, 17, 12, 13};
      d[7] = N22{0, 9, 1, 14};
      d[8] = N22{12, 13, -3, -71};
      d[9] = N22{9, -17, 0, 16};

      REQUIRE((sum(d) == N22{75, 28, 60, 21}));
    }
  }

  SECTION("checking for floating Array min/max")
  {
    SECTION("Min")
    {
      Array<double> b(10);
      b[0] = 13;
      b[1] = 1;
      b[2] = 8;
      b[3] = -3;
      b[4] = 23;
      b[5] = -1;
      b[6] = 13;
      b[7] = 0;
      b[8] = 12;
      b[9] = 9;

      REQUIRE(min(b) == -3);

      b.fill(0);
      REQUIRE(min(b) == 0);

      b[0] = -13;
      b[1] = -1;
      b[2] = -8;
      b[3] = -3;
      b[4] = -23;
      b[5] = -1;
      b[6] = -13;
      b[7] = 0;
      b[8] = -12;
      b[9] = -9;

      REQUIRE(min(b) == -23);

      b[0] = 13;
      b[1] = 1;
      b[2] = 8;
      b[3] = 3;
      b[4] = 23;
      b[5] = 1;
      b[6] = 13;
      b[7] = 0;
      b[8] = 12;
      b[9] = 9;

      REQUIRE(min(b) == 0);

      b[7] = 3;
      REQUIRE(min(b) == 1);
    }

    SECTION("Max")
    {
      Array<double> b(10);
      b[0] = 13;
      b[1] = 1;
      b[2] = 8;
      b[3] = -3;
      b[4] = 23;
      b[5] = -1;
      b[6] = 13;
      b[7] = 0;
      b[8] = 12;
      b[9] = 9;

      REQUIRE(max(b) == 23);

      b[0] = -13;
      b[1] = -12;
      b[2] = -8;
      b[3] = -3;
      b[4] = -23;
      b[5] = -1;
      b[6] = -13;
      b[7] = -10;
      b[8] = -12;
      b[9] = -9;

      REQUIRE(max(b) == -1);

      b.fill(-13);
      REQUIRE(max(b) == -13);

      b[0] = 13;
      b[1] = 12;
      b[2] = 8;
      b[3] = 3;
      b[4] = 23;
      b[5] = 1;
      b[6] = 13;
      b[7] = 10;
      b[8] = 12;
      b[9] = 9;
      REQUIRE(max(b) == 23);
    }
  }

  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);

      for (size_t i = 0; i < array.size(); ++i) {
        array[i] = ((i < (array.size() / 10)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
      }

      const auto sum_before_shuffle = sum(array);

      std::shuffle(&(array[0]), &(array[0]) + array.size(), std::mt19937{std::random_device{}()});

      const auto sum_after_shuffle = sum(array);

      REQUIRE(sum_before_shuffle == sum_after_shuffle);

      REQUIRE(sum_before_shuffle == Catch::Approx(direct_sum(array)));

      ReproducibleSumManager::setReproducibleSums(false);
      REQUIRE(sum(array) == Catch::Approx(direct_sum(array)));
      ReproducibleSumManager::setReproducibleSums(true);
    }

    SECTION("reproducible float sum")
    {
      Array<float> array(10'000);

      for (size_t i = 0; i < array.size(); ++i) {
        array[i] = ((i < (array.size() / 10)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
      }

      const auto sum_before_shuffle = sum(array);

      std::shuffle(&(array[0]), &(array[0]) + array.size(), std::mt19937{std::random_device{}()});

      const auto sum_after_shuffle = sum(array);

      REQUIRE(sum_before_shuffle == sum_after_shuffle);

      REQUIRE(sum_before_shuffle == Catch::Approx(direct_sum(array)).epsilon(1E-4));
    }

    SECTION("reproducible TinyVector<3,double> sum")
    {
      Array<TinyVector<3, double>> array(10'000);

      for (size_t i = 0; i < array.size(); ++i) {
        array[i] =
          TinyVector<3, double>(((i < (array.size() / 10)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
                                17 * i - 3,   //
                                1);
      }

      const auto sum_before_shuffle = sum(array);

      std::shuffle(&(array[0]), &(array[0]) + array.size(), std::mt19937{std::random_device{}()});

      ReproducibleTinyVectorSum s0(array);

      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")
    {
      Array<TinyVector<3, float>> array(10'000);

      for (size_t i = 0; i < array.size(); ++i) {
        array[i] =
          TinyVector<3, float>(((i < (array.size() / 10)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
                               17 * i - 3,   //
                               1);
      }

      const auto sum_before_shuffle = sum(array);

      std::shuffle(&(array[0]), &(array[0]) + array.size(), std::mt19937{std::random_device{}()});

      ReproducibleTinyVectorSum s0(array);

      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]).epsilon(1E-4));
      REQUIRE(sum_before_shuffle[1] == Catch::Approx(naive_sum[1]).epsilon(1E-4));
      REQUIRE(sum_before_shuffle[2] == Catch::Approx(naive_sum[2]).epsilon(1E-4));
    }

    SECTION("reproducible TinyMatrix<2, 3> sum")
    {
      Array<TinyMatrix<2, 3, double>> array(10'000);

      for (size_t i = 0; i < array.size(); ++i) {
        array[i] = TinyMatrix<2, 3, double>{((i < (array.size() / 1000)) * 1E25 + 1E10) * ((i + 1) % 1'000) *
                                              std::sin(3 * i + 1),
                                            17 * i - 3,            //
                                            1,                     //
                                            std::cos(7 * i) * i,   //
                                            (2 - 3 * i) * 1E-30,   //
                                            (i * 100 + 1000)};
      }

      const auto sum_before_shuffle = sum(array);

      std::shuffle(&(array[0]), &(array[0]) + array.size(), std::mt19937{std::random_device{}()});

      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")
    {
      Array<TinyMatrix<2, 3, float>> array(10'000);

      for (size_t i = 0; i < array.size(); ++i) {
        array[i] =
          TinyMatrix<2, 3, float>(((i < (array.size() / 1000)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1),
                                  17 * i - 3,            //
                                  1,                     //
                                  std::cos(7 * i) * i,   //
                                  (2 - 3 * i) * 1E-30,   //
                                  (i * 100 + 1000));
      }

      const auto sum_before_shuffle = sum(array);

      std::shuffle(&(array[0]), &(array[0]) + array.size(), std::mt19937{std::random_device{}()});

      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(5E-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));
        }
      }
    }

    SECTION("non reproducible double sum")
    {
      Array<double> array(10'000);

      for (size_t i = 0; i < array.size(); ++i) {
        array[i] = ((i < (array.size() / 10)) * 1E25 + 1E10) * ((i + 1) % 1'000) * std::sin(3 * i + 1);
      }

      ReproducibleSumManager::setReproducibleSums(false);
      REQUIRE(sum(array) == Catch::Approx(direct_sum(array)));
      ReproducibleSumManager::setReproducibleSums(true);
    }
  }

  SECTION("checking for subArrayView")
  {
    Array<int> array{10};
    for (size_t i = 0; i < array.size(); ++i) {
      array[i] = 2 * i + 1;
    }

    auto view = subArrayView(array, 3, 6);

    REQUIRE(view.size() == 6);

    bool use_same_memory = true;
    for (size_t i = 0; i < view.size(); ++i) {
      use_same_memory &= (&(view[i]) == &(array[i + 3]));
    }
    REQUIRE(use_same_memory);

    view.fill(-3);
    for (size_t i = 0; i < view.size(); ++i) {
      REQUIRE(view[i] == -3);
    }

    for (size_t i = 0; i < array.size(); ++i) {
      int int_i = i;
      if ((i >= 3) and (i < 9)) {
        REQUIRE(array[i] == -3);
      } else {
        REQUIRE(array[i] == 2 * int_i + 1);
      }
    }

    for (size_t i = 0; i < view.size(); ++i) {
      view[i] = 3 - 2 * i;
    }

    for (size_t i = 0; i < array.size(); ++i) {
      int int_i = i;
      if ((i >= 3) and (i < 9)) {
        REQUIRE(array[i] == 3 - 2 * (int_i - 3));
      } else {
        REQUIRE(array[i] == 2 * int_i + 1);
      }
    }

    std::ostringstream os;
    os << view;
    REQUIRE(os.str() == "0:3 1:1 2:-1 3:-3 4:-5 5:-7");
  }

#ifndef NDEBUG

  SECTION("output with signaling NaN")
  {
    Array<double> x{5};
    x[0] = 2;
    x[2] = 3;

    std::ostringstream array_ost;
    array_ost << x;
    std::ostringstream ref_ost;
    ref_ost << 0 << ':' << 2 << ' ';
    ref_ost << 1 << ":nan ";
    ref_ost << 2 << ':' << 3 << ' ';
    ref_ost << 3 << ":nan ";
    ref_ost << 4 << ":nan";
    REQUIRE(array_ost.str() == ref_ost.str());
  }

  SECTION("checking for bounds violation")
  {
    REQUIRE_THROWS_WITH(a[10], "invalid index");
  }

  SECTION("invalid copy_to")
  {
    Array<int> b{2 * a.size()};
    REQUIRE_THROWS_WITH(copy_to(a, b), "incompatible Array sizes");
  }

  SECTION("checking for nan initialization")
  {
    Array<double> array(10);

    for (size_t i = 0; i < array.size(); ++i) {
      REQUIRE(std::isnan(array[i]));
    }
  }

  SECTION("checking for bad initialization")
  {
    Array<int> array(10);

    for (size_t i = 0; i < array.size(); ++i) {
      REQUIRE(array[i] == std::numeric_limits<int>::max() / 2);
    }
  }

  SECTION("invalid UnsafeArrayView")
  {
    Array<int> array(10);

    REQUIRE_THROWS_WITH(subArrayView(array, 3, 8), "required view is not contained in the Array");
    REQUIRE_THROWS_WITH(subArrayView(array, 9, 2), "required view is not contained in the Array");
  }

  SECTION("check for bounds in UnsafeArrayView")
  {
    Array<int> array(10);
    auto view = subArrayView(array, 2, 8);

    REQUIRE_THROWS_WITH(view[8], "invalid index");
  }
#endif   // NDEBUG
}