While looking through some of my older programming projects, I found my old Egyptian Form project, which implements a greedy algorithm described by Solomon W. Golomb [1] for turning any fraction of positive integers into a sum of reciprocals of positive integers in C++. Since I want to stretch out my programming legs some more in preparation for (hopefully) entering Georgia Tech’s OMSCS program, I decided to see if there’s anything related to this that I could do some mathematics on while also programming a basic demonstration. That’s when I found the following conjecture:
Erdős–Straus Conjecture. For every integer $N\geq 2,$ there exist integers $0 \lt x\leq y \leq z$ such that $\frac{4}{N} = \frac{1}{x} + \frac{1}{y} + \frac{1}{z}.$
This conjecture is still open, but it does lend itself nicely to computational verification. And it also features some interesting optimizations that allow for huge performance gains, like taking something that originally didn’t finish after several minutes (I ended the process before it could finish) to taking just about 1 ms!
Note that this is a journey I went on on my own, so there may be better and more correct ways of computing these than what I’ll look at in this blog post. I haven’t spent much time looking at the literature for this problem, except for the identities in Mordell [2] mentioned on the Wikipedia page for the conjecture. This is just the result of a couple of days of experimentation, discovery, and programming.
A Working Algorithm
Let’s start by stating our goals.
- We want to find examples (or counterexamples) to this conjecture programmatically.
- We only care about finding if an example exists or not; we won’t try to enumerate them all.
- We want our code to be fast.
- We want to do a little math along the way to make things even faster!
I’ll be using the Number Theory Library (NTL) by Victor Shoup for the arbitrary integer size support. Below is the skeleton of our erdos_straus function.
void erdos_straus( NTL::ZZ const &N )
{ // Finding x,y,z such that 4/N = 1/x + 1/y + 1/z
if( N < 2 )
{
return;
}
// Code goes here!
std::cout << "Counterexample at N = " << N << "?\n";
}
To search for something, we need to know where to look (otherwise, we're outta luck!), so we have to nail down where our search space is. The naïve approach would be to just loop over all positive integers $x,y,z,$ but what's our stopping condition?
Proposition. Let $N\gt2$ be an integer. Suppose there exist integers $0\lt x\leq y\leq z$ satisfying the Erdős–Straus conjecture. Then,
- $\left\lceil \frac{N}{4} \right\rceil \leq x \leq \left\lfloor \frac{3N}{4} \right\rfloor$
- $x\leq y\leq \left\lfloor \frac{2Nx}{4x - N} \right\rfloor$
- $z = \frac{Nxy}{4xy - N\left(x + y\right)}$
Furthermore, if we find $x$ and $y$ in these intervals such that $z = \frac{Nxy}{4xy - N\left(x + y\right)}\geq y,$ then the Erdős–Straus conjecture is satisfied for $N.$ That is, we have $0\lt x\leq y\leq z$ such that $\frac{4}{N} = \frac{1}{x} + \frac{1}{y} + \frac{1}{z}.$
Proof.
1. The left side of the inequality is obvious, since $\frac{4}{N} = \frac{1}{x} + \frac{1}{y} + \frac{1}{z}\geq \frac{1}{x},$ and the ceiling comes for free.
For the right side, note that $x$ is maximized precisely when $y$ and $z$ are minimized, which would force $x = y = z.$ Therefore, $\frac{4}{N} \leq \frac{3}{x},$ which means $x\leq \frac{3N}{4}.$ Since $x$ is an integer, this gives us the floor for free, so $x\leq \left\lfloor \frac{3N}{4} \right\rfloor.$
2. The lower bound on $y$ is given by the conjecture's statement.
For the upper bound, maximizing $y$ minimizes $z,$ forcing $y = z,$ and so $\frac{4}{N} \leq \frac{1}{x} + \frac{2}{y}.$ Subtracting $\frac{1}{x}$ from both sides and combining fractions, we have $\frac{4x - N}{Nx} \leq \frac{2}{y}.$ Taking the reciprocal of both sides, multiplying by $2$ to isolate $y,$ and noting that $y$ being an integer gives us a free floor, we have $y\leq \left\lfloor \frac{2Nx}{4x - N} \right\rfloor,$ as required.
3. By basic algebraic manipulations, $$\frac{1}{z} = \frac{4}{N} - \frac{1}{x} - \frac{1}{y} = \frac{4x - N}{Nx} - \frac{1}{y} = \frac{4xy - Ny - Nx}{Nxy} = \frac{4xy - N\left(x + y\right)}{Nxy}.$$ Taking the reciprocal gives the desired result.
The final statement of the proposition follows from reversing the algebraic manipulation.
QED
Corollary. Given the same hypothesis as above, $\left(4xy - N\left(x + y\right)\right) \mid Nxy.$
Proof.
To satisfy the statement of Erdős–Straus, $z = \frac{Nxy}{4xy - N\left(x + y\right)}$ must be an integer.
QED
Not only do we now have a finite range of values to search, but we even have some conditions which will allow us to check that our values really do satisfy Erdős–Straus. We already have enough to write a function!
void erdos_straus( NTL::ZZ const &N )
{ // Finding x,y,z such that 4/N = 1/x + 1/y + 1/z
if( N < 2 )
{
return;
}
auto const x_min{ ( N + 3 ) / 4 };
// ceil(N/4) = floor((N+3) / 4)
auto const x_max{ ( 3 * N ) / 4 };
// floor(3N/4)
for( auto x{ x_min }; x <= x_max; ++x )
{ // ceil(N/4) <= x <= floor(3N/4)
auto const y_denominator{ 4 * x - N }; // 4x - N
if( y_denominator <= 0 )
{
continue;
}
auto const y_max{ ( 2*N*x ) / y_denominator };
// floor(2Nx / (4x-N))
for( auto y{ x }; y <= y_max; ++y )
{ // x <= y <= floor(2Nx / (4x-N))
auto const z_denominator{ 4 * x * y - N * ( x + y ) };
// 4xy - N(x+y)
if( z_denominator <= 0 )
{
continue;
}
auto const z_numerator{ N * x * y };
// Nxy
if( z_numerator % z_denominator != 0 )
{ // z must be an integer.
continue;
}
auto const z{ z_numerator / z_denominator };
// z = Nxy / (4xy - N(x+y))
if( z >= y )
{
std::cout << "4/" << N << " = 1/" << x
<< " + 1/" << y
<< " + 1/" << z << ".\n";
return;
}
}
}
std::cout << "Counterexample at N = " << N << "?\n";
}
Using std::chrono to time running erdos_straus( 10000 ) gives
4/10000 = 1/2501 + 1/6252501 + 1/39093762502500.
Finished in 1s305ms.
On one hand, this is a success! On the other hand, the 1.305 seconds on such a small input is a little concerning, since we really want to test up into the billions in a reasonable amount of time.
What's the runtime of our current algorithm? Since $x$ has to search a space of size $\left\lfloor \frac{3N}{4} \right\rfloor - \left\lceil \frac{N}{4} \right\rceil = \Theta\left(N\right)$ and $y$ also searches a space of size $\Theta\left(N\right)$, this algorithm takes $\Omega\left(N^{2}\right)$ time in the worst case (ignoring any time spent on arbitrary-size integer arithmetic, which is non-trivial). This is actually an exponential time, specifically pseudo-polynomial time!
We would really want to get a polynomial time, but if we can at least get a logarithmic factor like $\log N$ instead of one of the $N$'s, that would be nice. However, before that, we'll make a quick detour to OpenMP to make our current implementation parallel (this is what I did next while working on this problem).
First, let me describe the system I'm running on. I have a Ryzen 9 5950X with 32 GB DDR4-3600 @ 14-16-16-32. This gives me 16 cores, 32 threads, 64 MB of L3 cache, and a decent RAM-to-CPU latency, which will hopefully translate into a huge performance boost when we run this on multiple threads.
Since each $x$-value's loop works independently, we can trivially slice this in such a way that $k$ cores work on each $k$th value of $x$ until a solution is found.
void erdos_straus( NTL::ZZ const &N )
{ // Finding x,y,z such that 4/N = 1/x + 1/y + 1/z
if( N < 2 )
{
return;
}
auto const x_min{ ( N + 3 ) / 4 };
// ceil(N/4) = floor((N+3) / 4)
auto const x_max{ ( 3 * N ) / 4 };
// floor(3N/4)
bool found{ false };
auto thread_count = omp_get_max_threads();
#pragma omp parallel shared(found)
{
auto const thread_num = omp_get_thread_num();
for( auto x{ x_min + thread_num }; x <= x_max && !found; x += thread_count )
{ // ceil(N/4) <= x <= floor(3N/4)
auto const y_denominator{ 4 * x - N }; // 4x - N
if( y_denominator <= 0 )
{
continue;
}
auto const y_max{ ( 2*N*x ) / y_denominator };
// floor(2Nx / (4x-N))
for( auto y{ x }; y <= y_max && !found; ++y )
{ // x <= y <= floor(2Nx / (4x-N))
auto const z_denominator{ 4 * x * y - N * ( x + y ) };
// 4xy - N(x+y)
if( z_denominator <= 0 )
{
continue;
}
auto const z_numerator{ N * x * y };
// Nxy
if( z_numerator % z_denominator != 0 )
{ // z must be an integer.
continue;
}
auto const z{ z_numerator / z_denominator };
// z = Nxy / (4xy - N(x+y))
if( z >= y )
{
#pragma omp critical
if( !found )
{
found = true;
std::cout << "4/" << N << " = 1/" << x
<< " + 1/" << y
<< " + 1/" << z << ".\n";
}
break; // Can't return inside OpenMP block.
}
}
}
}
if( !found )
{
std::cout << "Counterexample at N = " << N << "?\n";
}
return;
}
Changes:
- I've added a boolean
foundwhich each loop will check for and will be changed totruewhen an example is found. - The variables
thread_count = omp_get_max_threads()andthread_num = omp_get_thread_num()do as they say; these help OpenMP cut the search space up for the cores. - The
#pragma omp parallel shared(found)line tells OpenMP to run the block in parallel. - The
#pragma omp criticaltells OpenMP that other cores shouldn't run the corresponding block while one is executing it -- this prevents a race condition on the variablefoundwhen changing it fromfalsetotrue. - Since we're no longer returning inside the loop, we wrap the "counterexample" print out in a conditional.
Running erdos_straus( 10000 ) now gives:
4/10000 = 1/2529 + 1/218020 + 1/17230393125.
Finished in 91ms.
A huge improvement already! Though a side effect is that sometimes we get different answers, like
4/10000 = 1/2532 + 1/197813 + 1/78259768125.
Finished in 84ms.
This okay, though, since this just means that a different slice of the search space got there first. Remember, we only care about existence, not exactly which solution we find. But this also means that we need a new benchmark.
What happens if we run erdos_straus( 150000 )?
4/150000 = 1/37532 + 1/43982813 + 1/3868975634803125.
Finished in 18s218ms.
Again, we're hitting a wall. This will be our new benchmark to beat.
What now?
Continued Fractions
Let's look at $\frac{4}{N} = \frac{1}{x} + \frac{1}{y} + \frac{1}{z}$ again. Since we're dealing with $x$ by slicing it and feeding it to multiple cores, let's see if we can optimize how we're dealing with $y$. This is where continued fractions come in.
Here, we'll take a continued fraction to be of the form $a_{0} + \cfrac{1}{a_{1} + \cfrac{1}{a_{2} +\cdots\vphantom{\cfrac{1}{a_{3}}}}}.$
We will write this more succinctly as $\left[ a_{0}; a_{1}, a_{2},\ldots\right],$ and we'll take as a known fact that any rational number can be written as a finite continued fraction.
Fix $x$ and let $S = \frac{4}{N} - \frac{1}{x} = \frac{4x - N}{Nx}.$ Let $\left[ a_{0}; a_{1}, a_{2},\ldots, a_{n}\right]$ be a finite continued fraction expansion for $S.$ Since $S\leq 1$ ($\frac{4}{N}\leq 2$ and equality forces $x = 1$), we can further ignore $a_{0},$ since it will always be either $0$ or $1.$ We'll be using convergents (and, further, continuants of these convergents) to search for $y$.
Definition. The convergents of a continued fraction are a sequence of elements defined by $C_{k} = \frac{A_{k}}{B_{k}},$ where the $A_{k}$ and $B_{k}$ are the (numerator and denominator, respectively) continuants of $C_{k},$ given by the recurrence $$A_{k+1} = a_{k+1}A_{k} + A_{k-1}$$ $$B_{k+1} = a_{k} B_{k} + B_{k-1}$$ with initial values $A_{-1} = 1, A_{0} = 0$ and $B_{-1} = 0, B_{0} = 1.$
N.B. The actual definition of a convergent is that the $k$th convergent is the result of cutting the continued fraction off at $a_{k}.$ I'm using this definition here because it's more convenient for us.
We will use the following two theorems, the first without proof.
Theorem. Let $S$ be any non-negative real number and let $C_{k} = \frac{A_{k}}{B_{k}}$ be the sequence of convergents for a continued fraction of $S.$ Then, $$\left|S - C_{k}\right| \leq \frac{1}{B_{k}B_{k+1}} \lt \frac{1}{B_{k}^{2}}.$$
Theorem. The continuants of a continued fraction for a non-negative real number grow exponentially in $k.$
Proof.
We'll prove this for the denominator continuants, which is relevant to us, as the numerator case is similar.
First, note that $B_{-1} = F_{0}$ and $B_{0} = F_{1},$ where $F_{k}$ is the Fibonacci sequence.
Now, suppose $B_{k}\geq F_{k+1}.$ Since $B_{k+1} = a_{k+1}B_{k} + B_{k-1}\geq B_{k} + B_{k-1} \geq F_{k+2},$ we've proven that $B_{k}\geq F_{k+1}$ by induction.
Finally, since $F_{k}{ \sim } \frac{\varphi^{k}}{\sqrt{5}},$ where $\varphi$ is the golden ratio, this means we've proven that $B_{k}$ grows exponentially, as required.
QED
Corollary. The convergents converge to $S.$
Now that we've looked at a few properties about continued fractions, let's go back to the actual problem we're trying to solve: finding $y$ efficiently.
Instead of searching through every potential $y$, why not restrict to just good candidates? And what better candidates than the denominator continuants, $B_{k}?$
When $S = \frac{1}{y} + \frac{1}{z},$ we have $\left|S - \frac{1}{y}\right| = \frac{1}{z},$ which, since $y\leq z,$ means that $S$ and $\frac{1}{y}$ should be "close" to each other. This suggests that maybe there's a $B_{k}$ that can act as our $y.$ And since we've established that the continuants of convergents grow exponentially, this substantially reduces our search space.
With that, we can write a (mostly) correct algorithm for finding solutions to Erdős–Straus questions! I'll write the algorithm out first, then we'll have to talk about something important. Don't mind the two auxiliary functions, continued_fraction and generate_continuants, that we'll be using -- they can be found in the GitHub repository for this post.
void erdos_straus( NTL::ZZ const &N )
{ // Finding x,y,z such that 4/N = 1/x + 1/y + 1/z
if( N < 2 )
{
return;
}
auto const x_min{ ( N + 3 ) / 4 };
// ceil(N/4) = floor((N+3) / 4)
auto const x_max{ ( 3 * N ) / 4 };
// floor(3N/4)
bool found{ false };
auto thread_count = omp_get_max_threads();
#pragma omp parallel shared(found)
{
auto const thread_num = omp_get_thread_num();
for( auto x{ x_min + thread_num }; x <= x_max && !found; x += thread_count )
{ // Considering 4/N = 1/x + 1/y + 1/z,
// we have 4/N - 1/x =
// (4x - N)/(Nx) = 1/y + 1/z
auto const numerator{ 4 * x - N };
auto const denominator{ N * x };
if( denominator <= 0 )
{
continue;
}
// Generate the denominator continuants for the convergents of the continued fraction
auto const cf{ continued_fraction( numerator, denominator ) };
auto const continuants{ generate_continuants( cf ) };
for( auto const &B : continuants )
{ // y ~= denominator of a convergent
// (1/y should be close to (4x - N)/(Nx))
auto const y{ B };
if( y < x || y >( 2 * N * x ) / numerator )
{ // If y is too small or big, skip it.
continue;
}
// 1/z = 4/N - 1/x - 1/y = numerator/denominator - 1/y
// = (y * numerator - denominator) / (y * denominator)
auto const z_numerator{ denominator * y };
auto const z_denominator{ numerator * y - denominator };
if( z_denominator <= 0 || z_numerator % z_denominator != 0 )
{
continue;
}
auto const z{ z_numerator / z_denominator };
if( z >= y )
{
#pragma omp critical
if( !found )
{
found = true;
std::cout << "4/" << N << " = 1/" << x
<< " + 1/" << y
<< " + 1/" << z << ".\n";
}
break; // Can't return inside OpenMP block.
}
}
}
}
if( !found )
{
std::cout << "Counterexample at N = " << N << "?\n";
}
return;
}
Let's try an example.
Recall that erdos_straus( 150000 ) took 18s218ms before. Now, with this optimization, we get
4/150000 = 1/37514 + 1/100483929 + 1/23559713203162500.
Finished in 1ms.
That's fast! In fact, we can solve something much larger now.
What about erdos_straus( 1000000000 )?
4/1000000000 = 1/283029552 + 1/2142244860 + 1/1475658593750250000000.
Finished in 16s261ms.
This is 2 seconds faster than the previous implementation, and that one was just on 150000. That's a pretty good optimization! However, there's an issue. I haven't stated a theorem to prove the correctness of this algorithm, and for good reason.
That's because this algorithm isn't correct. It fails for several values of $N,$ such as $2,3,5,6,7,8,9,13,$ and a few more. But the examples seem to get much less common as $N$ grows large, so I suspect that the following "theorem," while not perfect, works a "good amount" of the time.
Theorem. The erdos_straus algorithm above has a search space of size $\Theta\left(N \log N\right)$ when the fallback is unnecessary.
To fix this, whenever we get to the end of erdos_straus, we call a function called erdos_straus_fallback which is just the one we had before the continued fractions stuff. That means our worst-case analysis is the same, but it seems to be extremely uncommon. So uncommon, in fact, that I missed these bad inputs when I had some other miscellaneous optimizations, which I'll be covering in a bit. The code with the fallback implemented can be found in the GitHub repository for this post, but don't look yet because you'll spoil the last section of this post.
Miscellaneous Optimizations
While we've done this whole business with reducing the search space using continued fractions, you may have already noticed a very simple optimization.
Suppose $N$ satisfies Erdős–Straus with $0\lt x\leq y\leq z.$ Then, for any $k\geq 1,$ $kN$ satisfies Erdős–Straus with $0\lt kx\leq ky\leq kz.$ Just divide both sides in $\frac{4}{N} = \frac{1}{x} + \frac{1}{y} + \frac{1}{z}$ by $k.$ So maybe we can knock out a few small primes and use the pre-computed results to handle larger multiples.
Proposition. Suppose $N\geq 2$ is an integer.
- If $N$ is divisible by 2, then Erdős–Straus is satisfied by $$\left(x,y,z\right) = \left(\frac{N}{2}, N, N\right)$$
- If $N$ is divisible by 3, then Erdős–Straus is satisfied by $$\left(x,y,z\right) = \left(\frac{N}{3}, 2N, 2N\right)$$
- If $N$ is divisible by 5, then Erdős–Straus is satisfied by $$\left(x,y,z\right) = \left(\frac{2N}{5}, N, 2N\right)$$
- If $N$ is divisible by 7, then Erdős–Straus is satisfied by $$\left(x,y,z\right) = \left(\frac{2N}{7}, \frac{28N}{7}, \frac{28N}{7}\right)$$
Proof.
1. For $N = 2$, $$\frac{4}{N} = 2 = 1 + \frac{1}{2} + \frac{1}{2}.$$ When $2\mid N,$ $N = 2k$ for some positive integer $k,$ and so $$\frac{4}{N} = \frac{4}{2k} = \frac{1}{k} + \frac{1}{2k} = \frac{1}{2k} = \frac{1}{N/2} + \frac{1}{N} + \frac{1}{N}.$$
For the rest, I'll avoid the general case and just exhibit the special cases where $k = 1.$
2. $$\frac{4}{3} = 1 + \frac{1}{6} + \frac{1}{6}.$$
3. $$\frac{4}{5} = \frac{1}{2} + \frac{1}{5} + \frac{1}{10}.$$
4. $$\frac{4}{7} = \frac{1}{2} + \frac{1}{28} + \frac{1}{28}.$$
QED
We could continue with this, but each successive prime contributes less and less. Instead, we should expand the number of congruence classes we cover on smaller integers. For instance, covering multiples of 11 only helps with 1 out of every 11 integers, while covering two congruence classes modulo 3 (0 and 2) covers 2 out of every 3 integers.
In Mordell [2], calculations are performed which (when worked out) give solutions to several other congruence classes in terms of $N.$ In these next propositions, we will be looking at those calculations, adding a bit, and then turning that into the final piece of the C++ code!
Proposition. [2, Pg. 287, Equation (3)] Suppose $N\geq 2$ and there exist positive integers $a,b,c,d$ such that $aN + b + c = 4abcd,$ then $N\geq 2$ satisfies Erdős–Straus with $$x,y,z\in \left\lbrace bcd, Nabd, Nacd\right\rbrace,$$ in an order where $0\lt x\leq y\leq z$.
Proof.
By algebraic manipulation, $$\frac{1}{bcd} + \frac{1}{Nabd} + \frac{1}{Nacd} = \frac{Na + c + b}{Nabcd} = \frac{4abcd}{Nabcd} = \frac{4}{N}.$$
QED
This may look like it comes from nowhere, but this will actually enable us to find formulas for integers congruent to $2 \pmod{3},$ $2,3 \pmod{4},$ $1,2,3 \pmod{5},$ $3,5,6 \pmod{7},$ and $5\pmod{8}.$ You could probably get even more (I'll leave it as an exercise to the reader), but we'll stick with these. We'll actually cover a very large number of integers by just special cases, which I'll show with some numbers after we get the formulas.
I'll list these in the same order as they're listed in [2], since some of the arguments require that $N$ doesn't fall into previous cases. This will be a big wall of statements, but I think it'll be easy to understand once we get into the flow.
Proposition. [2, Pg. 287-288] Suppose $N\geq 2$ is an integer.
- If $N \equiv 3 \pmod{4},$ then Erdős–Straus is satisfied by $$\left(x,y,z\right) = \left(\frac{N + 1}{4}, \frac{N\left(N + 1\right)}{2}, \frac{N\left(N + 1\right)}{2}\right).$$
- If $N \equiv 5 \pmod{8},$ then Erdős–Straus is satisfied by $$\left(x,y,z\right) = \left(\frac{N + 3}{4}, \frac{N\left(N + 3\right)}{8}, \frac{N\left(N + 1\right)}{4}\right).$$
- If $N \equiv 2 \pmod{3},$ then Erdős–Straus is satisfied by $$\left(x,y,z\right) = \left(\frac{N + 1}{3}, N, \frac{N\left(N + 1\right)}{3}\right).$$
- If $N \equiv 3 \pmod{7},$ then Erdős–Straus is satisfied by $$\left(x,y,z\right) = \left(\frac{2N + 1}{7}, 2N, \frac{N\left(2N + 1\right)}{7} \right).$$
- If $N \equiv 5 \pmod{7},$ then Erdős–Straus is satisfied by $$\left(x,y,z\right) = \left(\frac{2\left(N + 2\right)}{7}, 2N, \frac{N\left(N + 2\right)}{7} \right).$$
- If $N \equiv 6 \pmod{7},$ then Erdős–Straus is satisfied by $$\left(x,y,z\right) = \left( \frac{2\left(N + 1\right)}{7}, 2N, \frac{2N\left(N + 1\right)}{7} \right).$$
- If $N \equiv 1 \pmod{5},$ then Erdős–Straus is satisfied by $$\left(x,y,z\right) = \left( \frac{4\left(N + 4\right)}{15}, 4N, \frac{N\left(N + 4\right)}{15} \right).$$
- If $N \equiv 2 \pmod{5},$ then Erdős–Straus is satisfied by $$\left(x,y,z\right) = \left( \frac{2\left(2N + 1\right)}{15}, 4N, \frac{4N\left(2N + 1\right)}{15} \right).$$
- If $N \equiv 3 \pmod{5},$ then Erdős–Straus is satisfied by $$\left(x,y,z\right) = \left( \frac{4\left(N + 2\right)}{15}, 4N, \frac{2N\left(N + 2\right)}{15} \right).$$
Proof.
We will be using the previous proposition for this, so remember $aN + b + c = 4abcd$ gives $x,y,z\in \left\lbrace bcd, Nabd, Nacd\right\rbrace.$
1. Setting $a = 2, b = 1, c = 1$ gives $2N + 2 = 8d,$ which can be rearranged into $N = 4d - 1\equiv -1 \pmod{4}\equiv 3 \pmod{4}$ and $d = \frac{N + 1}{4}.$ Therefore, we have $bcd = \frac{N+1}{4},$ $Nabd = \frac{2N\left(N+1\right)}{4} = \frac{N\left(N+1\right)}{2} = Nacd.$
Since $d$ is any positive integer, this covers all relevant cases where $N \equiv 3 \pmod {4}.$ We will omit lines similar to this one from all subsequent cases.
Since the $0 \pmod{2}$ case covers both $0,2 \pmod{4},$ we need only consider $N \equiv 1 \pmod{4}$ from now on.
2. Setting $a = 1, b = 1, c = 2$ gives $N + 3 = 8d,$ which can be rearranged into $N = 8d - 3\equiv -3 \pmod{8} \equiv 5 \pmod{8}$ and $d = \frac{N + 3}{8}.$ Therefore, we have $bcd = \frac{N+3}{4},$ $Nabd = \frac{N\left(N + 3\right)}{8},$ and $Nacd = \frac{2N\left(N + 3\right)}{8} = \frac{N\left(N + 3\right)}{4}.$
Since the results for the congruence classes modulo 2 and 4 cover everything except $1,5 \pmod{8},$ we need only consider $N \equiv 1 \pmod{8}$ from now on.
For the next cases, we rearrange the equation $aN + b + c = 4abcd$ to $aN + b = c\left(4abd - 1\right) = cq.$
3. Take $q = 3$ and $a = b = 1.$ This forces $d = 1$ and gives $N + 1 = 3c,$ which can be rearranged into $N = 3c - 1 \equiv -1 \pmod{3} \equiv 2 \pmod{3}$ and $c = \frac{N + 1}{3}.$ Therefore, $bcd = \frac{N + 1}{3},$ $Nabd = N,$ and $Nacd = \frac{N\left(N + 1\right)}{3}.$
Since we've covered $0,2 \pmod{3},$ we need only consider $N\equiv 1 \pmod{3}$ from now on.
4-6. For the $\pmod{7}$ cases, take $q = 7.$ Then, $4abd - 1 = 7,$ which means $abd \mid 2.$
4. Take $b = d = 1,$ then $a = 2,$ and so $aN + b = cq$ becomes $2N + 1 = 7c,$ which (using the fact that $2^{-1} \equiv 4 \pmod{7}$) can be rearranged into $N = 2^{-1}\left(7c - 1\right)\equiv -4 \pmod{7} \equiv 3 \pmod{7}$ and $c = \frac{2N+1}{7}.$ Therefore, $bcd = \frac{2N + 1}{7},$ $Nabd = 2N,$ and $Nacd = \frac{N\left(2N + 1\right)}{7}.$
5. Take $a = d = 1,$ then $b = 2,$ and so $aN + b = cq$ becomes $N + 2 = 7c,$ which can be rearranged into $N = 7c - 2\equiv -2 \pmod{7}\equiv 5\pmod{7}$ and $c = \frac{N + 2}{7}.$ Therefore, $bcd = \frac{2\left(N + 2\right)}{7},$ $Nabd = 2N,$ and $Nacd = \frac{2N\left(N+2\right)}{7}.$
6. Take $a = b = 1,$ then $d = 2,$ and so $aN + b = cq$ becomes $N + 1 = 7c,$ which can be rearranged into $N = 7c - 1\equiv -1 \pmod{7}\equiv 6\pmod{7}$ and $c = \frac{N + 1}{7}.$ Therefore, $bcd = \frac{2\left(N + 1\right)}{7},$ $Nabd = 2N,$ and $Nacd = \frac{2N\left(N+1\right)}{7}.$
Since we've covered $0,3,5,6 \pmod{7},$ we need only consider $N\equiv 1, 2, 4 \pmod{7}.$
7-9. For the $\pmod{5}$ parts, take $q = 15.$ Then, $4abd - 1 = 15,$ which means $abd = 4.$
We will be making use of the fact that $N \equiv 1 \pmod{3}$ to turn $N \pmod{15}$ statements into $N \pmod{5}$ statements.
7. Take $a = d = 1,$ then $b = 4,$ and so $aN + b = cq$ becomes $N + 4 = 15c,$ which can be rearranged into $N = 15c - 4 \equiv -4 \pmod{15} \equiv 11 \pmod{15},$ which gives $N \equiv 1 \pmod{5},$ and $c = \frac{N+4}{15}.$ Therefore, $bcd = \frac{4\left(N + 4\right)}{15},$ $Nabd = 4N,$ and $Nacd = \frac{N\left(N+4\right)}{15}.$
8. Take $a = d = 2,$ then $b = 1,$ and so $aN + b = cq$ becomes $2N + 1 = 15c,$ which (using the fact that $2^{-1} \equiv 8 \pmod{15}$) can be rearranged into $N = 2^{-1}\left(15c - 1\right) \equiv -8 \pmod{15} \equiv 7 \pmod{15},$ which gives $N \equiv 2 \pmod{5},$ and $c = \frac{2N + 1}{15}.$ Therefore, $bcd = \frac{2\left(2N + 1\right)}{15},$ $Nabd = 4N,$ and $Nacd = \frac{4N\left(2N+1\right)}{15}.$
9. Take $b = d = 2,$ then $a = 1,$ and so $aN + b = cq$ becomes $N + 2 = 15c,$ which can be rearranged into $N = 15c - 2 \equiv -2 \pmod{15} \equiv 13 \pmod{15},$ which gives $N \equiv 3 \pmod{5},$ and $c = \frac{N+2}{15}.$ Therefore, $bcd = \frac{4\left(N + 2\right)}{15},$ $Nabd = 4N,$ and $Nacd = \frac{2N\left(N+2\right)}{15}.$
I leave it as an exercise to the reader to show that other combinations of values for $a,b,d$ when $q = 15$ don't give any new identities which can be used.
QED
We program these optimizations into a function special_cases that gets called at the top of erdos_straus, only use the continued fractions stuff for things that these miss, and only when all else fails do we use the old and slow erdos_straus_fallback. See the GitHub repository for the final code.
With this, all that's left is to do some benchmarking. Instead of plugging a single number in, since we'll either hit a special case or it'll perform as before, let's see how this does over a range of numbers. We'll start by computing erdos_straus on all of the integers from 2 to 15 million.
Running erdos_straus on all integers from 2 to 15000000.
Finished in 1m14s184ms
In each case below, the numbers are those which were not caught in the cases above it.
Total number = 0 (mod 2): 7500000.
Total number = 3 (mod 4): 3750000.
Total number = 5 (mod 8): 1875000.
Total number = 0 (mod 3): 625000.
Total number = 2 (mod 3): 625000.
Total number = 0 (mod 7): 89286.
Total number = 3 (mod 7): 89286.
Total number = 5 (mod 7): 89285.
Total number = 6 (mod 7): 89286.
Total number = 0 (mod 5): 53572.
Total number = 1 (mod 5): 53571.
Total number = 2 (mod 5): 53571.
Total number = 3 (mod 5): 53571.
Total number of special cases: 14946428.
Total number of integers which had to run through the main algorithm: 53571.
Total number which used the fallback: 0.
I added some extra infrastructure for collecting data. Notice how, even going through the first 15 million integers, we only used the main algorithm 53572 times (that's ~0.36% of the time), and not once did we need to go into the fallback. Of course, we're making a lot of use of the special cases since we're looking at a bunch of non-primes, so here's a run on the first 1.5 million primes.
Building vector of primes.
Prime vector of size 1500000 took 4ms to build.
Running erdos_straus on each prime up to 23879519.
Finished in 1m37s727ms
In each case below, the numbers are those which were not caught in the cases above it.
Total number = 0 (mod 2): 1.
Total number = 3 (mod 4): 750237.
Total number = 5 (mod 8): 375008.
Total number = 0 (mod 3): 0.
Total number = 2 (mod 3): 187681.
Total number = 0 (mod 7): 0.
Total number = 3 (mod 7): 31275.
Total number = 5 (mod 7): 31230.
Total number = 6 (mod 7): 31179.
Total number = 0 (mod 5): 0.
Total number = 1 (mod 5): 23295.
Total number = 2 (mod 5): 23402.
Total number = 3 (mod 5): 23412.
Total number of special cases: 1476720.
Total number of primes which had to run through the main algorithm: 23279.
Total number which used the fallback: 0.
I have a file with the first 50 million primes, so I load the number I want of them into RAM inside of a std::vector and let it rip. Again, no fallback code usage, and this is up to 23879519. I want to make the following conjecture, though I haven't really tested enough to have any confidence in it.
Conjecture (Optimistic). There exists an integer $N\geq 2$ such that the continued fractions algorithm (without the fallback) correctly determines if Erdős–Straus is satisfiable for all integers $n\geq N.$
Colloquially, this would mean that the continued fractions algorithm is asymptotically correct, which would mean any failure cases could be handled separately. Of course, it could just be that I've covered enough special cases to make up for the failure of the continued fractions algorithm until we hit much larger numbers. For a more realistic conjecture, we have the following.
Conjecture (Probably More Realistic). The set of integers $N\geq 2$ such that the continued fractions algorithm (without the fallback) fails to determine if Erdős–Straus is satisfiable has natural density $0.$
In other words, as we get to larger integers, the probability that we find a "bad" $N$ gets smaller and smaller. Either way, I think this algorithm works pretty well for a couple of days staring at a Kindle Scribe with pages upon pages of calculations and derivations and typing up a bunch of code (then retyping it for this blog post).
Conclusion
Optimizing code can be a real headache sometimes. I spent longer than I care to admit trying to prove the correctness of the continued fractions algorithm before trying to plug 2 in without the special cases enabled -- a case of premature optimization being the root of all evil? Looking into whether the failures of the continued fraction algorithm stop at some point, or at least have low (0?) natural density, seems like an interesting topic for further research once I have a clearer mind and fresher eyes for this problem.
The main goal of this post was really to do some programming and then try to optimize the code a bit. Since I've gone from the naïve brute-force taking 18 seconds for N = 150 thousand to calculating for N = 1 billion in just 16 seconds, I think I've done my job. I will likely still work on the code occasionally when I feel like doing more, though GitHub is where that'll be happening.
Of course, as mentioned in the beginning, this work was done in the blind as far as algorithms for this conjecture are concerned. I would be surprised if there weren't faster algorithms available for this problem, written by people much more knowledgeable than I am and who have spent more than just a couple of leisurely days on it.
References
[1] "An Algebraic Algorithm for the Representation Problems of the Ahmes Papyrus," by Solomon W. Golomb, American Mathematical Monthly, Vol. 69, 1962
[2] Mordell, Louis J. (1967), Diophantine Equations, Academic Press, pp. 287–290