Euler's brick generator in C++

Question

I stumbled upon this unsolved problem of math named Euler's brick. I wrote a program in C++ to generate solutions for Euler's brick problem.

It searches from 1 to 10,000 in about 10 seconds. It took 10.977153 minutes to search from 1 to 50,000. How can I increase the performance further so I can search for larger limits? Are there any important optimizations that I am missing?

Here is my current code:

#include <iostream>
#include <vector>
#include <array>
#include <math.h>
#include <chrono>

inline bool is_whole(double x) {
    return x - int(x) == 0;
}

int main() {
    auto start = std::chrono::high_resolution_clock::now();
    std::vector<std::array<int, 3>> solutions = {};
    unsigned long int i = 1;
    unsigned long int limit = 10000;
    while (i < limit) {
        int a = i;
        for (int b = i; b < limit; b++) {
            if (is_whole(sqrt(double((a*a + b*b))))) {
                for (int c = b; c < limit; c++) {
                    if (
                        is_whole(sqrt(double((a*a + c*c)))) &&
                        is_whole(sqrt(double((b*b + c*c))))
                    ) {
                        std::array<int, 3> sol = {a, b, c};
                        solutions.push_back(sol);
                    }
                }
            }
        }
        ++i;
    }
    auto stop = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::microseconds>(stop - start);
    std::cout << solutions.size() << " solutions found in " << duration.count() << " microseconds..." << std::endl;
    for (std::array<int, 3> sol : solutions) {
        if (is_whole(sqrt(sol[0] * sol[0] + sol[1] * sol[1] + sol[2]))) {
            std::cout << "\n\nPerfect cuboid found!!" << "[" << sol[0] << ", " << sol[1] << ", " << sol[2] << "]\n" << std::endl; 
        }
        std::cout << "[" << sol[0] << ", " << sol[1] << ", " << sol[2] << "], ";
    }
    return 0;
}

Answer 1

I recommend including <cmath> rather than the deprecated <math.h>, which should only be used in bilingual (C/C++) contexts. The C++ version ensures that identifiers are properly namespaced (e.g. std::sqrt), helping make your code clearer.

I think it would be clearer to define a simple Box type instead of std::array<int,3>, preferably with unsigned dimensions rather than signed.

There's lots of conversion between signed and unsigned types (including a which is always equal to static_cast<int>(i)) that can be eliminated. Remove i and make a, b and c share a single type (name it, in case you need to change it later).

Overuse of std::endl - I don't think there's any need to flush the output stream anywhere before it's closed (after main() returns), so just replace all those with plain newlines ("\n"). And make sure we don't finish the program with a partially-written line (at best, that's annoying to users).

is_whole() appears subject to floating-point rounding errors. I'd replace it with an is_perfect_square() that accepts an integer type and can be exact.

This test looks incorrect:

        if (is_whole(sqrt(sol[0] * sol[0] + sol[1] * sol[1] + sol[2]))) {

I'm guessing you meant sol[2] * sol[2] as the last addend.

To make the search faster, we'll need to abandon brute-force search, and use one of the algorithms that generates Pythagorean triples, at least for the first two dimensions.

The Wikipedia entry for Perfect cuboid contains a lot of constraints that could prune the search.

Answer 2

Your algorithm seems a little naive (in the sense it doesn't make use of available a priori facts). My java implementation did LIMIT = 50k in ~14.5s.

A = 43344; B = 47300; C = 49665

A = 43440; B = 45612; C = 49775

14.448s

Process finished with exit code 0

• Precalc squares of [1 - limit] (for array lookup rather than multiplcation in loops) (makes a bit of difference)

• Get and filter the Pythagorean numbers (per https://en.wikipedia.org/wiki/Euler_brick, "At least two edges of an Euler brick are divisible by 4" and "At least two edges of an Euler brick are divisible by 3", so you can chop off a LOT of candidates early. It also reduces the number of sqrt()s. (about 2/3 faster than not filtering here -- 14.5s vs 42s)

In java:

    // Generate and filter Pythagorean Numbers
    // O(n^2)
    long[] square = IntStream.range(0, N + 1).asLongStream().map(x -> x * x).toArray();
    for (int a = 1; a < N; ++a) {
        Long aSquare = square[a];
        for (int b = a + 1; b < N; ++b) {
            // Because at least 2 of a, b, c need to be divisible by 3 or 4,
            // at least one of a or b must, too.
            if ((a % 4 == 0 || b % 4 == 0) && 
                    (a % 3 == 0 || b % 3 == 0)
            ) {
                Long bSquare = square[b];
                if (Math.sqrt(aSquare + bSquare) % 1 == 0) { // hacky java isInteger? check
                    Set<Integer> values = pythagoreanNumbers.computeIfAbsent(a, k -> new HashSet<>());
                    values.add(b);
                }
            }
        }
    }

• If you save the Pythagorean numbers with an implied c, you can store it as a map of (key = a an Integer, value = Set{b an Integer| a^2 + b^2 = c^2 and sqrt(c^2) an integer}, and you only need to check entries with ||set|| > 1 because EITHER it is not a candidate or a -> { b, c }, b -> { c } => (a,b), (b,c), (a,c) [aka, a Euler Brick] are in the map. That is a -> { b, c } is sufficient to verify it is a brick. However, there may be more than one brick for a given a: e.g, (748, 1989, 4080) and (748, 4080, 6336), so only remove ||B|| = 1 (not sure how much faster, as I didn't do it straight up...)

• Then check the triples. In java:

    Stream<Map.Entry<Integer, Set<Integer>>> candidates = pythagoreanNumbers.entrySet().stream().filter(e -> e.getValue().size() > 1).sorted(Map.Entry.comparingByKey());
  
    // check candidates
    // O(N^3)
    // ** Without a lot of checking, my feeling is this is really between O(lg n) and O(N^2) (because the outer N
    // (a's) tends to be MUCH smaller than N, and the innermost N^2 (b's and c's) is even smaller on average
    // (~ lg N? after filtering?)), but using O(N^3) based on # of nested loops
    candidates.forEach(e -> {
        Integer[] objects = e.getValue().stream().sorted().toArray(Integer[]::new);
        Integer a = e.getKey();
        for (int i = 0; i < objects.length - 1; ++i) {
            Integer b = objects[i];
            for (int j = i + 1; j < objects.length; ++j) {
                Integer c = objects[j];
                if ((a % 11 == 0 || b % 11 == 0 || c % 11 == 0) && // At least one edge of an Euler brick is divisible by 11 
                                                                   // Maybe not so useful, but put in for 'completeness' of using known facts and because I suspect % is enough more efficient than sqrt() when x gets big 
                        (Math.sqrt(square[b] + square[c]) % 1 == 0) // hacky java isInteger? check again
                ) {
                    System.out.println("A = " + a + "; B = " + b + "; C = " + c);
                }
            }
        }
    });
  

Answer 3

In the outer loop, for example, you check all pairs a > b. Let's say a = 10,000 then you check 10,000 values for b. Loop over the length of the diagonal instead. The diagonal goes from 10,000 to 14,142 or so, so there are only 4,142 different values. You instantly save a factor 2.5.

Now at each step you calculate a square root and check whether it is an integer. You don't need to do that. Every time you you change d by 1, the square root while change by some value which changes slowly. It will be a bit tricky to work out. But with the example a = 10,000, the change in b will be less than 0.5 at d = 10,020. So you do your calculation until the change in b changes by less than 0.5. At d = 10,018, b changes by 16.927. At d = 10,019, b changes by 16.462. So you calculate b rounded to an integer ones. You increase it by 17, compare d^2 - a^2 and b^2, if b^2 is larger than you decrease b by 1, and if both sides are the same then you found a solution. As d grows, you have to change the increment in b to 16, 15, 14 etc. So usually you will have only two checks involving integers only. And d^2 - b^2 changes by 2d+1 if you increase d by 1, so that is very simple.

If you get to much bigger numbers: You can use the gcd algorithm to do many steps in one. For example with a = 1,000,000, d = 1,200,000, b increases by 1.809 if you increase d by 1, and that increase changes very, very slowly. So you run until you have a square root only slightly larger than an integer (when b^2 is slightly less than d^2 - a^2), and then you can increase d by 5 and b increases by 9.045. You can do that several times, until the fractional part is more than 0.809, then you increase d by 4. The number of times you need to increase by 4 will change slowly as the difference between square roots changes. All this is quite tricky, but it can reduce the number of steps to something like a^(2/3) instead of a.