PytharogianRectangles.cpp

// The algorithm itself
xt::xarray<int> next_pytharogian_triangles(xt::xarray<int> const& previous_stage)
{
   static const xt::xarray<int> stacked_matrices = {
      { -1, 2, 2 },
      { -2, 1, 2 },
      { -2, 2, 3 },
      { 1, 2, 2 },
      { 2, 1, 2 },
      { 2, 2, 3 },
      { 1, -2, 2 },
      { 2, -1, 2 },
      { 2, -2, 3 }
   };

   auto shape = previous_stage.shape();
   xt::xarray<int> next_three = xt::transpose(xt::linalg::dot(stacked_matrices, xt::transpose(previous_stage)));
   next_three.reshape({ 3 * shape[0], shape[1] });
   return next_three;
}

// To emulate a co-routine
class PytharogianRectangles
{
public:
   PytharogianRectangles()
      : m_currents({{ 3, 4, 5 }})
   {}

   xt::xarray<int> const& currents() const {
      return m_currents;
   }

   void next() {
      m_currents = next_pytharogian_triangles(m_currents);
   }

private:
   xt::xarray<int> m_currents;
};

// A basic bench
int bench()
{
   PytharogianRectangles rectangles;
   int count = 0;
   for (int depth = 0; depth < 10; ++depth) {
      auto triples = rectangles.currents();
      count += triples.shape()[0];
      if (depth + 1 < 10)
         rectangles.next();
   }
   return count;
}