Searching for three positive integers as a solution to an equation

Question

I recently ran into a question on M.SE asking,

For positive integers x, y, z find a solution s.t. \$ \frac{x}{y+z} + \frac{y}{x+z} + \frac{z}{x+y} = 4 \$

so rather than think, I made a quick brute force c++ program to check for solutions:

unsigned lim = 251;
for(double x=1; x<lim; x++) {
    for(double y=1; y<lim; y++) {
        for(double z=1; z<lim; z++) {
            printf("x=%.0f, y=%.0f, z=%.0f\n",x,y,z);
            if( std::abs((x/(y+z) + y/(x+z) + z/(y+x)) - 4) <= 1E-15) {
                printf("solution: x=%.0f, y=%.0f, z=%.0f\n",x,y,z);
            }
        }
    }
}

All of this being within main, of course. This is obviously slow so I was looking for optimizations. A couple of my thoughts were to remove the std::abs call because the inside expression must be positive, perhaps I could simply check == 4 rather than account for precision. I couldn't figure out a way to drop a for loop, because you can't isolate a variable. Otherwise I'm not sure what to do.

Questions:

• How can this code be sped up? My goal is lim=1E3 in under a minute. As of now it takes 31.0643 s for lim=251.
• Can the for-loops be reduced? As of now formatting is not really an issue for me, but can this be simplified down to a single expression and/or loop? Or maybe look a little cleaner, without sacrificing optimization?

Update: After some research I'm now aware that the smallest known solution to this has numbers with 81 digits. I'm not concerned with this, just the above questions.

Answer 1

Here are some simple optimizations you can make that can push your limit a bit higher:

1. Don't use floating point arithmetic. Floating point arithmetic is slow compared to integer arithmetic (this is not always true with the advent of things like FPUs, but it is a good rule of thumb). Even worse, floating point arithmetic is prone to things like rounding errors, which means you might find solutions which don't work and miss solutions that do work.

   One way to avoid using floating point arithmetic here is to expand out the equation into the form:

   \$x(x+y)(x+z) + y(y+x)(y+z) + z(z+x)(z+y) = 4(y+z)(x+z)(x+y)\$

   In this form you can check whether the equation is satisfied just using integer arithmetic (make sure you use a integer type large enough so that neither side overflows, however).

2. Make use of symmetry (thanks @Deduplicator). Another simple observation you can make is that this equation is symmetric in \$x\$, \$y\$, and \$z\$, i.e. if you have a solution \$(x, y, z)\$, then any permutation of this solution also works. One way to use this fact is to only loop over \$1 \leq y \leq x\$ and \$1 \leq z \leq y\$. This cuts down on your total number of iterations by a factor of \$3! = 6\$.

3. Reduce number of degrees of freedom. One interesting observation is that, once you fix \$x\$ and \$y\$, you don't have many choices left for \$z\$. So instead of looping over all \$z\$ in the range \$[1, L]\$, you can more efficiently compute the possible values of \$z\$. (This removes one of your for loops).

   How do you find the possible values of \$z\$? Well, if we know what \$x\$ and \$y\$ are, the above equation reduces to a cubic in \$z\$; i.e., something for the form \$c_3z^3 + c_2z^2 + c_1z + c_0 = 0\$. There are a couple ways to solve this for \$z\$; the easiest way I can think of (but definitely not the fastest) is to try all the factors of \$c_0\$ (by the Rational Root Theorem, any integer \$z\$ that satisfies this must divide \$c_0\$). This does involve factoring \$c_0\$, however, which can take a while depending on how you implement it.

   A faster method is to use binary search to find the roots, but you have to be a bit careful here, since the cubic might not be monotone increasing/decreasing over the entire interval. The correct way to do this is to first find the roots of the derivative of the cubic (the derivative of the cubic is a quadratic, so you can do this with just the quadratic formula); these are the critical points of the cubic, so the cubic is monotone increasing/decreasing on the intervals between these points, and you can binary search here.

4. Learn about elliptic curves. So it turns out that the smallest solutions to this Diophantine equation are immense, and it's pretty much impossible to find any via any method like this (I forget how large they are, but I believe each of x, y, and z is at least 10 digits long). The question then becomes, how did people find this answer in the first place?

   The trick here is that the cubic above (in fact, pretty much any homogeneous cubic in three variables) is an instance of a mathematical object known as an elliptic curve, and that these objects have a ton of nice properties. One of these nice properties is that there is a method whereby you can take any two integer solutions to this equation and "add" them to get a third integer solution. So one approach which works is to start with an integer solution (but not a positive integer solution) such as (11, 9, -5), and repeatedly "add" it to itself until you end up with a positive integer solution (as far as I know, there's no guarantee you ever will end up with a positive integer solution, but in this case you do after a couple of steps).

EDIT:

I ran some tests to compare the above optimizations (see here). On my machine for limit=1000:

• OP's code (removing the printf per iteration) takes ~14.5 seconds.
• OP's code with optimization 1 takes ~4.2 seconds.
• OP's code with optimizations 1 + 2 takes ~0.7 seconds.
• OP's code with optimizations 1 + 2 + 3 takes ~0.25 seconds.

For larger limits, these differences are more pronounced. For limit=5000:

• OP's code with optimizations 1 + 2 + 3 takes ~7.9 seconds.
• OP's code with optimizations 1 + 2 takes ~86 seconds.
• The other two cases each took at least 5 minutes.

Answer 2

Stop printing so much

I ran your program and found that it spent most of its time printing every single combination (even though none of them were solutions). After I removed the line that prints the current triplet, I was able to run your program in 0.31 seconds (limit = 251).

I then extended the limit to 1001, and the program ran in 19.91 seconds, but it didn't find any solution.