Searching for three positive integers as a solution to an equation

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.

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.
• Thanks! I actually just found the first known solution has x, y, z of the order 10^81 :) Do you know of any implementations of the root binary search algorithm you mentioned that you could link to? Commented Jun 14, 2017 at 19:34
• There are a lot of references online which describe how to find the roots of arbitrary-degree polynomials (e.g. en.wikipedia.org/wiki/…) that use some similar ideas but are inherently more complicated than the cubic case. If I get some time later I'll edit my answer with an implementation of the cubic root-finding optimization for this problem. Commented Jun 14, 2017 at 20:15
• One additional point, iff you have a solution (x,y,z), than all permutations are also solutions. So the search-space can be reduced to x <= y <= z Commented Jun 14, 2017 at 20:31
• Good point @Deduplicator; I'll add that to my answer Commented Jun 14, 2017 at 20:59
• @cat: this is true, and I sort of alluded to this with what I said about FPUs, but in this case point 1 still holds, just because (integer) multiplication is a much faster operation than (float) division. Perhaps I should have said "Don't use division" instead of "Don't use floating point arithmetic", but there are other reasons for not using floating point arithmetic that I pointed out (i.e. precision issues). Commented Jun 15, 2017 at 1:47

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.

• yeah I figured as much. Usually I put the print statements in there, because I go crazy and think the program has hung if there is no output. Commented Jun 14, 2017 at 19:21
• @Dando18 A good solution for that is to print occasionally (e.g. if (z % 100 == 0) /* do the print message */). Even though the modulo costs a little bit, it's still pretty cheap and you won't have to wonder if the program has hung. What JS1 is pointing out is that I/O is one of the most expensive things you can do. Commented Jun 14, 2017 at 21:14
• @Justin thanks. in this case I moved the print statement to just the outer loop and it significantly helped. Commented Jun 14, 2017 at 21:16