Pythagorean Triple finder Would like some feedback on a code style exercise motivated by this PythagTriple algorihtm. I started with the Rosetta Code - "efficient C" solution  (not the naive one) and refactored it heavily into modern C++17. 

It doesn't quite match their "Task spec". That that was not the focus. However, it could very easily be adopted to do so - just a few changes in main().Performance appears very similar to the C code I started with and the generated ASM seems reasonable. 

But the main focus here is style, ease of readability, maintenance and correctness. 

Specific feedback points I am looking for:


Use of the lambda to process each triple
Way of constraining the template to prevent unsuitable lambda's being passed. (C++20 concepts needed to do it cleanly?)
Use of initialisers for static data. Also is there a way to get away from the {{ .. {{ .. {},{},{} .. }} .. }}  craziness? 4 extra opening and closing braces. 
Use of specific type for vec3 and a using alias for trans
Body of the transform() function. Is there a map() ... splat way? 
General function signatures / use of const, constexpr, attributes etc
Use of "by value" vs "by ref" semantics for triples / tranforms etc. There are some "arguably redundant" const T&amp;s in there which were just as fast as "by value" (probably copy elided). 
More algorithmic question: The algorithm as implemented is constrained by max_perimeter and due to recursive DFS flow and the matrix transforms the triples come out in a (subjectively) "weird order". Several other constraints and orderings are conceivably possible. One idea is to implement several alternative .find_XX() methods to do that. 


EDIT: There is updated code in a separate answer below, integrating the feedback. 

#include &lt;array&gt;
#include &lt;iostream&gt;

namespace pythag {

class triple {
 public:
  const long _a;
  const long _b;
  const long _c;

  template &lt;typename F&gt;
  constexpr void find(long max_perim, F&amp;&amp; proc) {
    if (perimeter() &gt; max_perim) return;
    proc(*this);
    for (auto&amp; T : U) transform(T).find(max_perim, proc);
  }

  [[nodiscard]] constexpr long perimeter() const noexcept { return _a + _b + _c; }

  friend std::ostream&amp; operator&lt;&lt;(std::ostream&amp; stream, const triple t) noexcept {
    return stream &lt;&lt; t._a &lt;&lt; ", " &lt;&lt; t._b &lt;&lt; ", " &lt;&lt; t._c;
  }

private:
  struct vec3 {  // a strong type to distinguish from triple
    const int _x;
    const int _y;
    const int _z;
  };

  using trans = std::array&lt;vec3, 3&gt;;

  [[nodiscard]] constexpr long dot(vec3 V) const noexcept {
    return V._x * _a + V._y * _b + V._z * _c;
  }

  [[nodiscard]] constexpr triple transform(const trans&amp; T) const noexcept {
    return triple{dot(T[0]), dot(T[1]), dot(T[2])};
  }

  static constexpr auto U = std::array&lt;trans, 3&gt;{{
      // https://en.wikipedia.org/wiki/Pythagorean_triple#Parent.2Fchild_relationships
      // clang-format off
      {{{  1, -2,  2},     // vec3 U[0][0]
        {  2, -1,  2},     // vec3 U[0][1]
        {  2, -2,  3}}},   // vec3 U[0][1]

      {{{  1,  2,  2},     // vec3 U[1][0]
        {  2,  1,  2},     // vec3 U[1][1]
        {  2,  2,  3}}},   // vec3 U[1][2]

      {{{ -1,  2,  2},     // vec3 U[2][0]
        { -2,  1,  2},     // vec3 U[2][1]
        { -2,  2,  3}}},   // vec3 U[2][2]
      // clang-format on
  }};
};

} // namespace pythag

int main() {
  // basic usage demo
  long sum_peri = 0;
  long count = 0;

  pythag::triple{3, 4, 5}.find(
      100'000'000, // produces 7'023'027 "primitive" triples in &lt;100ms on i7 2600
      // [](const auto&amp; t) { std::cout &lt;&lt; t &lt;&lt; "\n"; });  // print option. slow obviously
      [&amp;sum_peri, &amp;count](const auto&amp; t) { sum_peri += t.perimeter(); ++count; });

  std::cout &lt;&lt; count &lt;&lt; "\n";        // these 2 lines are just a way to prevent
  return (sum_peri ^ count) &amp; 0xff;  // entire programme being optimised away
}
Here are some comments on your revised code, many of which apply to your original code as well.

for (auto&amp; T: U) transform(T).find(max_perim, proc);


You're using the names T and U to refer to things that are not template type parameters. In fact, U is a reference to a static data member of the enclosing class, whose only declaration appears many lines below this use. This is extremely surprising code. I might instead write

for (auto&amp;&amp; formula : this-&gt;child_formulas) {
    this-&gt;transform(formula).find(max_perim, proc);
}


Notice the use of braces around the loop body; the use of auto&amp;&amp;, Which Always Works; the use of this-&gt; to clarify for the human reader what scope we're expecting to find these names in; and especially the descriptive name child_formulas for U. (child_formulas says what U is; this-&gt; says whose it is. Both of those things were in question, in the old code.)



class triple : public std::array&lt;long, 3&gt; {


Don't inherit from standard types. It's legal C++, but it's bad practice.


  If you (or your project, or your company) didn’t write class Foo, then class Foo should not be granted control over the API of your own class. And that’s what you’re doing when you inherit from a base class: you’re granting that class control over your API.


I like your original x, y, z data members much better.



template &lt;typename F&gt;
void find(long max_perim, F&amp;&amp; proc) {


You're passing proc by perfect-forwarding, but then when you use it, you're using it as a plain old lvalue. This is fine (in the wild-west anything-goes world of templates), but it doesn't express the meaning of proc quite as well as I'd like. proc is supposed to be a callback that gets called for each triple, right? Calling a callback doesn't modify the callback. So, we can and should require that proc be const-callable, and then we just write

template&lt;class F&gt;
void find(int max_perim, const F&amp; proc) {


(Drive-by irrelevant style adjustments to save typing. long varies in size and is the same size as int on many platforms, so let's just use int until we're ready to step all the way up to a well-defined int64_t or __int128_t.)

If you really want to support non-const and/or rvalue-callable Fs, you can perfect-forward proc all over the place, but trust me, it's not worth it. (The std::forward&lt;F&gt;(proc)s will clutter your code, and all you're doing is enabling your users to write confusing and unintuitive client code.)

Vice versa, re your question about concepts: You can constrain this template to SFINAE away in C++17 (or C++11) like this:

template&lt;class F, class = std::enable_if_t&lt;std::is_invocable_v&lt;const F&amp;, triple&amp;&gt;&gt;&gt;
void find(int max_perim, const F&amp; proc) {


This is not much more boilerplate than the C++2a concepts version:

template&lt;class F&gt; requires std::is_invocable&lt;const F&amp;, triple&amp;&gt;
void find(int max_perim, const F&amp; proc) {


Notice that in C++2a you will also be able to write "simply"

void find(int max_perim, const std::invocable&lt;triple&amp;&gt; auto&amp; proc) {


but (A) this doesn't express exactly the same SFINAE-constraint as the other two, and (B) IMHO it looks more confusing and scary.

But should you constrain this function to SFINAE away when it's given a "bad" lambda type? IMHO, no, you shouldn't. There's no reason to SFINAE here. SFINAE is for when you need this version of find to drop out of the overload set so that some other more general version can pick up the call. That's not the situation that you have, here.

Compare the error messages you get from (A) (or the C++2a Concepts version (A2))

template&lt;class F, class = std::enable_if_t&lt;std::is_invocable_v&lt;const F&amp;, triple&amp;&gt;&gt;&gt;
void find(int max_perim, const F&amp; proc) {
    proc(*this);
}

this-&gt;find(42, 7);


versus (B)

template&lt;class F&gt;
void find(int max_perim, const F&amp; proc) {
    static_assert(std::is_invocable_v&lt;const F&amp;, triple&amp;&gt;);
    proc(*this);
}

this-&gt;find(42, 7);


versus (C)

template&lt;class F&gt;
void find(int max_perim, const F&amp; proc) {
    proc(*this);
}

this-&gt;find(42, 7);


and see which version feels most "user-friendly" for the client programmer. (Remember that someone using find correctly will never see any of these error messages, so you should optimize to help the guy who doesn't know how to use it.)



  cnt_next_depth += U.size(); // always 3
  if (--cnt_this_depth == 0) {
    if (++depth &gt; max_depth) return;
    cnt_this_depth = cnt_next_depth;
    cnt_next_depth = 0;
  }
  for (auto&amp; T: U) q.push(t.transform(T));


Code comments are usually a red flag, at least on CodeReview. ;) It seems that this += corresponds to the three calls to q.push below; i.e., your algorithm requires that cnt_next_depth exactly track q.size(). But you can't use q.size() because you are reusing q to store both the elements at the current level and the elements at the next level.

It would make more sense to use two different queues:

std::vector&lt;triple&gt; this_level = ...;
std::vector&lt;triple&gt; next_level;
for (int i=0; i &lt; max_depth; ++i) {
    for (auto&amp;&amp; t : this_level) {
        proc(t);
        for (auto&amp;&amp; formula : child_formulas) {
            next_level.push_back(t.transform(formula));
        }
    }
    this_level = std::move(next_level);
    next_level.clear();
}


As a bonus, it turns out that we don't even need them to be queues anymore; they can be plain old vectors, and we save a bunch of heap traffic.



friend std::ostream&amp; operator&lt;&lt;(std::ostream&amp; stream, const triple&amp; t) noexcept {
  // Frustrating boiler plate. Any terse alternatives, that do it quickly and correctly?
  char comma[] = {'\0', ' ', '\0'}; // NOLINT
  for (auto d: t) {
    stream &lt;&lt; comma &lt;&lt; d;
    comma[0] = ',';
  }
  return stream;
}


No good answer. The idiomatic way would be

friend std::ostream&amp; operator&lt;&lt;(std::ostream&amp; stream, const triple&amp; t) {
    bool first = true;
    for (auto&amp;&amp; d : t) {
        if (!first) stream &lt;&lt; ", ";
        stream &lt;&lt; d;
        first = false;
    }
    return stream;
}


Notice that this function is not noexcept, because any of these stream operations might throw. (For example, if you're writing to a stringstream and it throws bad_alloc.)

Also, in real life I'd write for (int d : t) to indicate exactly what type of "elements" I expect to be getting out of a triple. I'm using auto&amp;&amp; here only because I saw your comment that you want to stick with "Almost Always Auto" style.



  triple transform(const trans&amp; T) const noexcept {
    auto res = triple{};
    std::transform(T.begin(), T.end(), res.begin(), [this](vec3 V) {
      return std::inner_product(V.begin(), V.end(), this-&gt;begin(), 0L);
    });
    return res;
  }


This is interesting use of STL algorithms, but algorithms are really meant for operating on big anonymous ranges of data elements, not for constant-size situations like we have here. In order to use STL algorithms, you've been forced to anonymize your formulas into faceless data ranges:

  static constexpr auto U = std::array&lt;trans, 3&gt;{{
      // https://en.wikipedia.org/wiki/Pythagorean_triple#Parent.2Fchild_relationships
      {{{{{1, -2, 2}},    // vec3 U[0][0]
         {{2, -1, 2}},    // vec3 U[0][1]
         {{2, -2, 3}}}}}, // vec3 U[0][1]

      {{{{{1, 2, 2}},    // vec3 U[1][0]
         {{2, 1, 2}},    // vec3 U[1][1]
         {{2, 2, 3}}}}}, // vec3 U[1][2]

      {{{{{-1, 2, 2}},    // vec3 U[2][0]
         {{-2, 1, 2}},    // vec3 U[2][1]
         {{-2, 2, 3}}}}}, // vec3 U[2][2]
  }};


Compare that code to a more "code-driven" version:

triple transform(const Formula&amp; f) const noexcept {
    return f(x, y, z);
}

auto formulas[] = {
    +[](int x, int y, int z){ return triple{ x - 2*y + 2*z, 2*x - y + 2*z, 2*x - 2*y + 3*z}; },
    +[](int x, int y, int z){ return triple{ x + 2*y + 2*z, 2*x + y + 2*z, 2*x + 2*y + 3*z}; },
    +[](int x, int y, int z){ return triple{ -x + 2*y + 2*z, -2*x + y + 2*z, -2*x + 2*y + 3*z}; },
};


In this version, trans is gone already, and t.transform(f) is hardly pulling its weight.



if (perimeter() &gt; max_perim) return;


Consider what you'd do if you wanted to find all triples up to some maximum side length, or up to some maximum hypotenuse length. Does this approach generalize to

triple.find(callback_on_each_triple, stopping_condition);


?