9
\$\begingroup\$

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;
}
\$\endgroup\$
13
  • 3
    \$\begingroup\$ No perfect Euler brick has been found, and it appears that this is the unsolved part. \$\endgroup\$ Jan 17, 2023 at 15:41
  • 1
    \$\begingroup\$ This bit of math will help you quickly generate triples worth examining: en.wikipedia.org/wiki/Tree_of_primitive_Pythagorean_triples \$\endgroup\$
    – J_H
    Jan 17, 2023 at 17:46
  • \$\begingroup\$ Hold on. On the linked page is says no perfect Euler brick has been found but then list the smallest one: The smallest Euler brick has sides (a,b,c)=(240,117,44) and face polyhedron diagonals d_(ab)=267, d_(ac)=244, and d_(bc)=125. What am I missing? \$\endgroup\$ Jan 17, 2023 at 19:38
  • 1
    \$\begingroup\$ Rather than have O(N^3). Could you reduce this to O(N^2) + O(N). 1: Find all integer pairs that have a perfect diagonal. 2: Match pairs with other pairs with a common side. This is the CS solution. Best would be to come up with a Maths solution that has a non brute force technique. \$\endgroup\$ Jan 17, 2023 at 19:44
  • 1
    \$\begingroup\$ no comment on the code, but it was proven in 2019 that no perfect Euler brick exists \$\endgroup\$
    – HDawG
    Jan 18, 2023 at 13:14

3 Answers 3

11
\$\begingroup\$

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.

\$\endgroup\$
2
  • \$\begingroup\$ The timing seems to report the wrong number as well. \$\endgroup\$
    – pacmaninbw
    Jan 17, 2023 at 16:14
  • \$\begingroup\$ Thanks... I tried adding the Pythagorean triplet idea... it works... but does not generate all the solutions... It also goes way above the limit when generating solutions due to the 2 * m * n. I have added the new code... Kindly please guide me... \$\endgroup\$ Jan 18, 2023 at 1:13
2
\$\begingroup\$

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);
                  }
              }
          }
      });
    
\$\endgroup\$
1
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.