Pythagorean Triple finder

Question

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&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 <array>
#include <iostream>

namespace pythag {

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

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

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

  friend std::ostream& operator<<(std::ostream& stream, const triple t) noexcept {
    return stream << t._a << ", " << t._b << ", " << 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<vec3, 3>;

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

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

  static constexpr auto U = std::array<trans, 3>{{
      // 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 <100ms on i7 2600
      // [](const auto& t) { std::cout << t << "\n"; });  // print option. slow obviously
      [&sum_peri, &count](const auto& t) { sum_peri += t.perimeter(); ++count; });

  std::cout << count << "\n";        // these 2 lines are just a way to prevent
  return (sum_peri ^ count) & 0xff;  // entire programme being optimised away
}

Answer 1

Use of the lambda to process each triple

Sure, why not. Nothing wrong with using a lambda in main() here instead of writing a normal function.

Way of constraining the template to prevent unsuitable lambda's being passed. (C++20 concepts needed to do it cleanly?)

If you could live with pythag::triple::find() not being constexpr, then you don't need C++20 for this at all. Just use std::function<> for the type of proc:

#include <functional>
...
void find(long max_perim, std::function<void(triple &)> proc) {
    ...
}

Use of initialisers for static data. Also is there a way to get away from the {{ .. {{ .. {},{},{} .. }} .. }} craziness? 4 extra opening and closing braces.

Yes, by not using std::array<>, but just declaring a multidimensional, "C-style" array:

  static constexpr vec3 U[3][3] {
      {{  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]
  };

You need to change transform() as well:

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

Use of specific type for vec3 and a using alias for trans

I think you could've made vec3 an array as well, either std::array<> (at the cost of even more braces) or just int[3]. Since you use trans in multiple places, it's good to have made an alias for it. But what's weird is that you write:

static constexpr auto U = std::array<trans, 3>{...};

The auto is totally unnecessary here, you could've written it simply as:

static constexpr std::array<trans, 3> U{...};

Body of the transform() function. Is there a map() ... splat way?

There's std::inner_product which could replace dot(), and std::transform() could have replaced the body of transform(), if your triples and vectors were arrays. It would have looked like:

triple transform(const trans& T) const {
    triple result;
    return std::transform(T.begin(), T.end(), result.values.begin(),
        [this](vec3 V){
            return std::inner_product(V.begin(), V.end(), this->values.begin(), 0);}
        }
    );
}

Instead of const long _a, _b, _c you'd have to have something like std::array<long, 3> values to make it work.

General function signatures / use of const, constexpr, attributes etc

Does it really need to be constexpr? It is restricting you (you can't use std::function<> for find() for example), and I don't see any reason why you would ever need to know the number of Pythogorean triples up to some value at compile-time.

Of course, marking things const, constexpr, [[nodiscard]] and so on where possible is a good thing. Keep doing that.

I would however change the type of the components of vec3 to long, to match that of the triple itself. You could even think of making an alias for the value type, or make the whole class templated, so you can decide which value type to use.

Use of "by value" vs "by ref" semantics for triples / tranforms etc. There are some "arguably redundant" const T&s in there which were just as fast as "by value" (probably copy elided).

With optimization enabled, any decent compiler will probably generate the same assembly code here, regardless of whether you pass by value or by const reference. You even missed the possibility to make operator<<() take a const reference.

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.

What algorithm is best depends on the application. Your recursive DFS algorithm is probably the most elegant one to implement, but apart from the weird output order, the main issue is that you can run out of stack space if the recursion goes too deep. So try to implement another algorithm that doesn't require recursive function calls.

Answer 2

Here are some comments on your revised code, many of which apply to your original code as well.

for (auto& 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&& formula : this->child_formulas) {
    this->transform(formula).find(max_perim, proc);
}

Notice the use of braces around the loop body; the use of auto&&, Which Always Works; the use of this-> 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-> says whose it is. Both of those things were in question, in the old code.)

---

class triple : public std::array<long, 3> {

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 <typename F>
void find(long max_perim, F&& 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<class F>
void find(int max_perim, const F& 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<F>(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<class F, class = std::enable_if_t<std::is_invocable_v<const F&, triple&>>>
void find(int max_perim, const F& proc) {

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

template<class F> requires std::is_invocable<const F&, triple&>
void find(int max_perim, const F& proc) {

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

void find(int max_perim, const std::invocable<triple&> auto& 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<class F, class = std::enable_if_t<std::is_invocable_v<const F&, triple&>>>
void find(int max_perim, const F& proc) {
    proc(*this);
}

this->find(42, 7);

versus (B)

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

this->find(42, 7);

versus (C)

template<class F>
void find(int max_perim, const F& proc) {
    proc(*this);
}

this->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 > max_depth) return;
    cnt_this_depth = cnt_next_depth;
    cnt_next_depth = 0;
  }
  for (auto& 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<triple> this_level = ...;
std::vector<triple> next_level;
for (int i=0; i < max_depth; ++i) {
    for (auto&& t : this_level) {
        proc(t);
        for (auto&& 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& operator<<(std::ostream& stream, const triple& t) noexcept {
  // Frustrating boiler plate. Any terse alternatives, that do it quickly and correctly?
  char comma[] = {'\0', ' ', '\0'}; // NOLINT
  for (auto d: t) {
    stream << comma << d;
    comma[0] = ',';
  }
  return stream;
}

No good answer. The idiomatic way would be

friend std::ostream& operator<<(std::ostream& stream, const triple& t) {
    bool first = true;
    for (auto&& d : t) {
        if (!first) stream << ", ";
        stream << 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&& here only because I saw your comment that you want to stick with "Almost Always Auto" style.

---

  triple transform(const trans& 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->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<trans, 3>{{
      // 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& 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() > 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);

?

Answer 3

An answer to document how to integrate some of the excellent feedback points from @G. Sliepen in his accepted answer above.

The following have been changed from the original question above:

• Remove constexpr since it is not required and opens up options for using the STL algorithms (except for the static U, where constexpr allows inline definition of this constant)
• Use std::array<long|int, 3> to store both the triple and the vec3 transform data. This allows more algorithmic use than ._a, ._b, ._c and ._x, ._y, ._z. (But see below for how I retained simple accessor syntax)
• Use std::inner_product and std::transform to do the matrix multiplication. std::accumulate for the perimeter.
• Addition of a "Breadth first search" (BFS) finder.

Things I decided not to adopt:

• std::function, which is known to be very slow, was deemed inappropriate for inner loop lambda. So the constraining of the lambda type remains a TODO.
• C-Arrays: Lose their size when you pass them around (since Kernighan and Ritchie) and decay to a pointer. So std::array<> adds value as we can easily use it in algorithms, because it knows its size and you can call std::begin() etc. We learn to live with the "many braces".

The one extra wrinkle was to make triple inherit from std::array<long, 3> and vec3 inherit from std::array<int, 3>. Then similarly for trans. This allows strong typing (ie not just a using alias which just decays away) but retains the easy and direct access via operator[] and usage in algorithms or ranged for loops without writing a single line of code. There is no vtable or similar overhead, because this is not virtual, polymorphic inheritance.

Performance is identical to the original above. Code size looks bigger (goldbolt link), but that is all due to the new BFS algo which uses std::queue, heap allocators etc. If we remove that for a fairer comparison the code size is identical to above, and speed on -O3 is identical. ie STL std::inner_product and std::transform and std::accumulate are truly "zero cost" here, and arguably "negative cost" given the improved readability and more general algorithm.

An unrelated annoyance is the boilerplate in operator<<. Still haven't found a terse way of expressing that.

#include <algorithm>
#include <array>
#include <cstdio>
#include <functional>
#include <iostream>
#include <numeric>
#include <queue>

namespace pythag {

class triple : public std::array<long, 3> {
public:
  template <typename F>
  void find(long max_perim, F&& proc) {
    // recursive DFS
    if (perimeter() > max_perim) return;
    proc(*this);
    for (auto& T: U) transform(T).find(max_perim, proc);
  }

  template <typename F>
  void find_by_level(int max_depth, F&& proc) {
    // iterative BFS with level tracking / limiting
    std::queue<triple> q;
    q.push(*this);
    int depth          = 0;
    int cnt_this_depth = 1;
    int cnt_next_depth = 0;
    while (!q.empty()) {
      auto t = q.front();
      q.pop();
      proc(t);
      cnt_next_depth += U.size(); // always 3
      if (--cnt_this_depth == 0) {
        if (++depth > max_depth) return;
        cnt_this_depth = cnt_next_depth;
        cnt_next_depth = 0;
      }
      for (auto& T: U) q.push(t.transform(T));
    }
  }

  [[nodiscard]] long perimeter() const noexcept {
    return std::accumulate(this->begin(), this->end(), 0L);
  }

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

private:
  struct vec3 : public std::array<int, 3> {};   // strong types by inheritance
  struct trans : public std::array<vec3, 3> {}; // to distinguish from triple

  [[nodiscard]] triple transform(const trans& 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->begin(), 0L);
    });
    return res;
  }

  static constexpr auto U = std::array<trans, 3>{{
      // 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]
  }};
};

} // namespace pythag

int main() {
  using pythag::triple;

  // basic usage demo
  auto print = [](const triple& t) { std::cout << t << "\n"; };
  std::cout << "Primitive triples up to a perimeter of 200\n";
  triple{{3, 4, 5}}.find(200, print);
  std::cout << "\nPrimitive triples up to 2 levels of transformation away from {3,4,5}\n";
  triple{{3, 4, 5}}.find_by_level(2, print);

  // // performance test
  // long sum_peri = 0;
  // long count    = 0;
  // pythag::triple{{3, 4, 5}}.find(
  //     100'000'000, // produces 7'023'027 triples in <100ms on i7 2600
  //     [&sum_peri, &count](const auto& t) {
  //       sum_peri += t.perimeter();
  //       ++count;
  //     });

  // // prevent entire programme being optimised away without <iostream>
  // return (sum_peri ^ count) & 0xff;
}