# 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&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), dot(T), dot(T)};
}

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
{  2, -1,  2},     // vec3 U
{  2, -2,  3}}},   // vec3 U

{{{  1,  2,  2},     // vec3 U
{  2,  1,  2},     // vec3 U
{  2,  2,  3}}},   // vec3 U

{{{ -1,  2,  2},     // vec3 U
{ -2,  1,  2},     // vec3 U
{ -2,  2,  3}}},   // vec3 U
// 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
}



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 {
{{  1, -2,  2},     // vec3 U
{  2, -1,  2},     // vec3 U
{  2, -2,  3}},    // vec3 U

{{  1,  2,  2},     // vec3 U
{  2,  1,  2},     // vec3 U
{  2,  2,  3}},    // vec3 U

{{ -1,  2,  2},     // vec3 U
{ -2,  1,  2},     // vec3 U
{ -2,  2,  3}},    // vec3 U
};


You need to change transform() as well:

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


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. 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.

• Thanks for your comments. Yes. std::array<>. I can never decide, when I want a "small, non-dynamically allocated" set of pretty trivial types which container to use. std::vector is overkill. C-style arrays feel "dirty" (eventhough that's my history). So I end up with std::array.... well.. at least it remembers the size, and that was useful at least once here. So yes, changing (back) to C-style [] arrays, opens up a couple of options. Other tradeoff which you are touching on is "strong types" vs "primitives" ie struct { long,long,long}. vs std::array. There are costs to both. tbc Jan 1, 2020 at 1:01
• Re: don't need auto. I was trying to stick to Herb Sutter's "auto to stick" style here. I am still trying to decide whether I like it. Jan 1, 2020 at 1:04
• Re: "does it need to be constexpr"? No, ;-) it doesn't. But at least I am getting to the point where I am "overremembering" to put it in. ;-) And you are right, it resticts things slightly. I did almosty get it to the point of calling from main with constexpr, it can almost do it. Jan 1, 2020 at 1:07
• Re: std::function. I have previously benchmarked this, and it really is seriously slow. So I tend to stay away from it almost as a matter of policy. And that lamda is in the inner loop... Jan 1, 2020 at 1:12
• Re: operator<< and const triple& That's a good spot. My analyser didn't help me spot this omission, probably because it is a friend... Jan 1, 2020 at 1:17

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 = ',';
}
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
{{2, -1, 2}},    // vec3 U
{{2, -2, 3}}}}}, // vec3 U

{{{{{1, 2, 2}},    // vec3 U
{{2, 1, 2}},    // vec3 U
{{2, 2, 3}}}}}, // vec3 U

{{{{{-1, 2, 2}},    // vec3 U
{{-2, 1, 2}},    // vec3 U
{{-2, 2, 3}}}}}, // vec3 U
}};


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


?

• Thanks! Some really useful stuff there. Especially on the constraining of the lambda template. I admit I was spraying around the && because I know that "works", but was looking for input exactly in this area. I choose T & U because these are Matrixes and it is common in linear algebra to denote matrices with capitals. But I agree there is some "domain overlap" ther with template params. Not sure I agree with the "formula approach". Frankly 3x3 Matrix multiplication, which is all this is, should be a bulk standard operation (without using Eigen et al)... TBC Jan 1, 2020 at 22:09
• Your point on inheritance is well made. However I sort of feel a bit stuck with no where to go. a) use plain std::array<> for everything plus free functions perhaps some using aliases. Perhaps the best for what is will remain a minimalist academic article. b) My "cheaty middle way" where I get some typing by inheritance and the API I need. But tightly bound to STL type and intheriting stuff I don't need. c) The third way is fully stand alone types with their own iterators etc (unless we go back ot doing matrix multiply with a,b,c * ,x,y,z which works fine for the minimal example. Jan 1, 2020 at 22:21
• I feel the a,b,c with a page of manually typed formulae (BTW very very similar to the original C code from Rosetta code) is missing the point that these are matrices and the operations are standardised. At the end of the day, it's an architecture question which will depend on the real life use case. And this will never have one. So it's just speculative opinion really? That's perhaps why 2 people have recommended going opposite ways on it already? Jan 1, 2020 at 22:24
• I think "faceless data ranges" perhaps illustrates that you didn't see these as, or didn't realise they were matrices and standard linear algebra matrix operations? Jan 1, 2020 at 22:25
• On the BFS level limit algo. Sure. "Not my code". I already didnt't like it when I copy pasted that very piece ;-) In fact I think it's lame that we need a heap allocating queue or vector at all to implement this algo. Jan 1, 2020 at 22:28

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.

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 = ',';
}
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
{{2, -1, 2}},    // vec3 U
{{2, -2, 3}}}}}, // vec3 U

{{{{{1, 2, 2}},    // vec3 U
{{2, 1, 2}},    // vec3 U
{{2, 2, 3}}}}}, // vec3 U

{{{{{-1, 2, 2}},    // vec3 U
{{-2, 1, 2}},    // vec3 U
{{-2, 2, 3}}}}}, // vec3 U
}};
};

} // 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;
}


• Turning operator<< into a one-liner: return stream << std::accumulate(t.begin(), t.end(), std::string{}, [](const std::string &a, long b){return (a.empty() ? a : a + ", ") + std::to_string(b);}); Jan 1, 2020 at 16:01
• Yeah. The advantage of the char[] trick (which came from a cppreference example) is that it is branchless. The entire problem of join, as I would call it, is not that trivial and I would have hoped had been solved in the STL. the std::accumulate way also results in a ton of sting concatenations (with potential allocations) and then piping out. Anyway, probably a topic for another post. Jan 1, 2020 at 16:14
• @G.Sliepen this would be good: return os::str::join(stream, t); except it should say std::join instead of os::str::join. Internally that calls begin / end and uses them to feed the stream, inserting a "glue" and "term" which can be passed as optional params. There are 4 different overloads of it. I have a 300 line header file full of these little 1-10 line functions. 80% of them are string related. Am I missing a trick / lib or does everyone have one of those? Jan 1, 2020 at 17:20