Algorithm for Egyptian fractions, want to optimize the maximum denominator size

Any rational number can be expanded into a finite sum of unit fractions with distinct denominators, called Egyptian fractions.

There is no 'optimal' algorithm in terms of denominator size or number of fractions.

Some of the best known algorithms:

• 'Splitting' method, based on the relation: $\frac{1}{n}=\frac{1}{n+1}+\frac{1}{n(n+1)}$.

All three of these methods usually give very large maximum denominators. The best overview of different algorithms can be seen here.

I tried to make an algorithm which gives smaller maximum denominators. The algorithms can be called Splitting-Joining (SJ), because it employs a combination of:

• adding fractions to reduce the size and length - for example $\frac{1}{8}+\frac{1}{8}=\frac{1}{4}$;
• if the first step is not possible, then we split the repeating fractions according to the rule $\frac{1}{a b}=\frac{1}{a(a+b)}+\frac{1}{b(a+b)}$, where $a,b$ are two divisors of $n$ closest in size, but not equal. For example: $\frac{1}{21}+\frac{1}{21}=\frac{1}{21}+\frac{1}{30}+\frac{1}{70}$.

This is the idea anyway.

Here is my naive realization of the algorithm.

p=5;
q=37;
Q1=Table[q,{j,1,p}];
Q=Split[Sort[Q1]];
M=Length[Q];
P=Table[Length[Q[[j]]],{j,1,M}];
While[Total[P]>M,
Q=Split[Sort[Q1]];
M=Length[Q];
P=Table[Length[Q[[j]]],{j,1,M}];
G=Table[GCD[Q[[j,1]],P[[j]]],{j,1,M}];
G1=Table[GCD[Q[[j,1]],P[[j]]-1],{j,1,M}];
A=Table[Q[[j,1]],{j,1,M}];
B=Table[Q[[j,1]],{j,1,M}];
Do[
If[P[[j]]!=1,
If[G[[j]]==1,If[G1[[j]]==1,
C0=Reverse[Divisors[Q[[j,1]]]];
L0=Length[C0];
A[[j]]=C0[[Floor[L0/2]]];
B[[j]]=Q[[j,1]]/A[[j]];
If[A[[j]]==B[[j]],
A[[j]]=C0[[Floor[L0/2]-1]];
B[[j]]=Q[[j,1]]/A[[j]]];
Q[[j]]=Join[{Q[[j,1]]},Table[A[[j]] (A[[j]]+B[[j]]),{k,1,P[[j]]-1}],Table[B[[j]] (A[[j]]+B[[j]]),{k,1,P[[j]]-1}]],
Q[[j]]=Join[{Q[[j,1]]},Table[Q[[j,1]]/G1[[j]],{k,1,(P[[j]]-1)/G1[[j]]}]]],Q[[j]]=Table[Q[[j,1]]/G[[j]],{k,1,Floor[P[[j]]/G[[j]]]}]]],{j,1,M}];
Q1=Flatten[Q]];


Here I use very convenient Split procedure on the list of denominators, combining them into groups of identical values.

Then I look if I can add all the fractions in each group, or at least all but one to obtain a shorter list of smaller values. I do it using the GCD function on the value of the denominator and the size of the group.

If this fails, I leave one denominator the same, and split others according to the rules above. I use the Divisors procedure to get $a$ and $b$ values.

Then I flatten the list of all denominators and repeat everything untill there is no repetitions in the list at all.

I'm not very good at Mathematica so far, so I would appreciate some pointers (I'm not asking for anything more) about improving the algorithm.

How to make it more efficient? Shorter? Are there bugs in this version?

The other question (which might be off-topic for this community):

How to optimize the length of the expansion, without enlarging the maximum denominator? Sometimes this algorithm gives very large expansions, compared to others. Maybe I missed something.

Here are the results, compared to Greedy and Engel algorithms (the sequence of denominators of unit fractions is provided as the result):

$$\frac{3}{37}$$

SJ: $\{19, 37, 703\}$

Greedy: $\{13, 241, 115921\}$

Engel: $\{13, 247, 9139\}$

$$\frac{5}{121}$$

SJ: $\{61, 62, 121, 3782, 7381, 11102, 22022\}$

Greedy: $\{25, 757, 763309, 873960180913, 1527612795642093418846225\}$

Engel: $\{25, 775, 31775, 1938275, 234531275\}$

$$\frac{5}{122}$$

SJ: $\{61, 62, 122, 3782\}$

Greedy: $\{25, 1017, 3101850\}$

Engel: $\{25, 1025, 125050\}$

$$\frac{10}{39}$$

SJ: $\{7, 13, 39, 91\}$

Greedy: $\{4, 156\}$

Engel: $\{4, 156\}$

$$\frac{1023}{1024}$$

SJ: $\{2, 4, 8, 16, 32, 64, 128, 256, 512, 1024\}$

Greedy: $\{2, 3, 7, 44, 9462, 373029888\}$

Engel: $\{2, 4, 8, 16, 32, 64, 128, 256, 512, 1024\}$

$$\frac{41}{182}$$

SJ: $\{20, 28, 30, 35, 65, 70, 78, 91, 130, 140, 182, 260\}$

Greedy: $\{5, 40, 3640\}$

Engel: $\{5, 40, 3640\}$

$$\frac{18}{23}$$

SJ: $\{3, 4, 12, 23, 24, 69, 78, 552, 598\}$

Greedy: $\{2, 4, 31, 2852\}$

Engel: $\{2, 4, 32, 736\}$

So far I see that the size of the maximum denominator is the best for SJ, however the length of the expansion sometimes suffers. And moreover, the denominators are not the best possible. For example we can actually expand:

$$\frac{18}{23}=\frac{1}{2}+\frac{1}{4}+\frac{1}{46}+\frac{1}{92}$$

This one I've done using the same ideas, but doing it by hand, which allowed me to see the patterns my program missed. The same expansion is also obtained by the Binary Remainder Method, which seems to also give the better denominator sizes.

Edit

The Binary Remainder Method seems to be mostly superior to the current version of SJ. For the listed above numbers it gives the following results:

$\frac{3}{37}$ - $\{16, 74, 296, 592\}$ (SJ gives {19, 37, 703});

$\frac{5}{121}$ - $\{32, 121, 968, 1936, 3872\}$ (SJ gives {61, 62, 121, 3782, 7381, 11102, 22022});

$\frac{5}{122}$ - $\{32, 122, 976, 1952\}$ (SJ gives {61, 62, 122, 3782});

$\frac{10}{39}$ - $\{4, 156\}$ (SJ gives {7, 13, 39, 91});

$\frac{41}{182}$ - $\{8, 16, 32, 182, 1456, 2912\}$ (SJ gives {20, 28, 30, 35, 65, 70, 78, 91, 130, 140, 182, 260});

Some are worse, but most are better. Still, I think SJ method could use some improvement.

Edit 2:

I found out that splitting the whole group instead of leaving one element unchanged makes the algorithm work better in some cases (but not others):

If we change:

Q[[j]] = Join[{Q[[j, 1]]},   Table[A[[j]] (A[[j]] + B[[j]]), {k, 1, P[[j]] - 1}], Table[B[[j]] (A[[j]] + B[[j]]), {k, 1, P[[j]] - 1}]]


to:

Q[[j]] = Join[Table[A[[j]] (A[[j]] + B[[j]]), {k, 1, P[[j]]}], Table[B[[j]] (A[[j]] + B[[j]]), {k, 1, P[[j]]}]]


We get great improvement for $5/22$ ({6, 22, 66}), $5/122$ ({31, 122, 1891}), $18/23$ ({2, 4, 46, 92}) and some others. But in some cases the result is less optimal.

Another big problem - the current realization of the algorithm is very slow when it comes to large numerators and denominators.

• Some code for the "greedy algorithm" is in here. Would you be okay with an (improved) implementation of the "binary remainder" method? Feb 3 '17 at 20:15
• @J.M., the implementation for binary remainder is available at the page I linked. So I already have it. I'm interested in improving this new method, because sometimes it gives better results than BR. See the edit especially. But if you have improved implementation of BR, I'm interested too Feb 3 '17 at 20:17
• To be clear - both BR and SJ methods do not give the best results in many cases. For example, using another simple method I found for $5/122$: {44,88,242,484,968} and for $41/182$: {8,14,56,91}. But since this is Mathematica community, I'm asking specifically about how to improve SJ algorithm (especially the code itself) Feb 6 '17 at 8:23
• *In the comment above I meant $5/121$ Feb 6 '17 at 13:08
• What do you really want to optimize for with respect to denominator size? After reading through your question, I can tell you want to minimize the maximum denominator, but you're also concerned with expansion length. From a "mathematically aesthetic" viewpoint, it seems to me (a matter of opinion) that the most elegant expansion is the one in which the sum of the denominators is lowest, regardless of the number of unit fractions. This is an interesting problem. Feb 9 '17 at 7:49

As to making the code slightly more efficient and shorter, I can point you in the right direction:

 FaroFrac2[p_, q_] :=
Module[{Q, M, P, G, G1, A, B, C0, L0, PQ, AB},
Q = Table[q, {p}];
M = 1;
P = {p};
While[Total[P] > M,
Q = Split[Sort[Q]];
M = Length[Q];
P = Length[#] & /@ Q;
PQ = Transpose[{P, Q}];
G = GCD[#[[2, 1]], #[[1]]] & /@ PQ;
G1 = GCD[#[[2, 1]], #[[1]] - 1] & /@ PQ;
A = #[[1]] & /@ Q;
B = A;
Do[
If[P[[i]] != 1,
If[G[[i]] == 1,
If[G1[[i]] == 1,
C0 = Reverse[Divisors[Q[[i, 1]]]];
L0 = Length[C0];
A[[i]] = C0[[Floor[L0/2]]];
B[[i]] = Q[[i, 1]]/A[[i]];
If[A[[i]] == B[[i]],
A[[i]] = C0[[Floor[L0/2] - 1]];
B[[i]] = Q[[i, 1]]/A[[i]];
];
AB = Table[A[[i]] + B[[i]], {j, P[[i]] - 1}];
Q[[i]] =
Join[{Q[[i, 1]]}, (A[[i]]*#) & /@ AB, (B[[i]]*#) & /@ AB],
Q[[i]] =
Join[{Q[[i, 1]]},
Table[Q[[i, 1]]/G1[[i]], {j, (P[[i]] - 1)/G1[[i]]}]];
],
Q[[i]] = Table[Q[[i, 1]]/G[[i]], {j, Floor[P[[i]]/G[[i]]]}];
]
],
{i, 1, M}
];
Q = Flatten[Q];
];
Q
]


• Removed Q1 because it was a superfluous variable
• Simplified the initial variable instantiation
• Changed Table function calls to Map (/@) wherever possible
• Initialized B as a copy of A instead of doing the same computation twice
• In the "Split" computation, I created a list AB to save the additions so they could be mapped across the multiplications without having to perform the same addition twice.

These are (in my view) about the only efficiency-improvements that can be made without starting to change the basic structure of your current algorithm, as you posted it. There very well may be vast improvements to be made by reimplementing your math with new code to better utilize the functional aspects of Mathematica.

These minor changes yielded an equally minor improvement to the run-time:

For testing purposes I wrapped your original code in a function called FaroFrac, while the new version is named FaroFrac2.

You also asked about bugs: well, after checking the output on millions of inputs, I'd say it's looking pretty bug free. Of course, if p > q, it fails badly, but I assume you are not considering that a valid input. Here's a little test code I used to double check output.

FaroSimp[D_List] := Total[(1/#) & /@ D]
testTable =
Table[{RandomInteger[1000],
1000 + RandomInteger[1000000]}, {1000}];
faroRes = Map[{FaroFrac2[#[[1]], #[[2]]], #} &, testTable];
deltas = FaroSimp[#[[1]]] - #[[2, 1]]/#[[2, 2]] & /@ faroRes;
Select[deltas, # != 0 &]


This should give you a good idea of the "surface" Mathematica side of things; but as to your other question, regarding improvements to the mathematical nature of the problem, that will take much more and deeper consideration.

• Thank you. You've gone above and beyond what I expected. Now that I think about it any mathematical improvements of the algorithm are off-topic for this community. One small comment: the code with Q[[i]] = Join[(A[[i]]*#) & /@ AB, (B[[i]]*#) & /@ AB] instead of Q[[i]] =Join[{Q[[i, 1]]}, (A[[i]]*#) & /@ AB, (B[[i]]*#) & /@ AB] still gives better results for most fractions (for some reason), i.e. shorter expansions with smaller maximum denominator - wait - how do I get the same result as in my Edit 2. This simple change doesn't give correct results Feb 11 '17 at 7:51
• Never mind, I just forgot to change the definition of $AB$. Anyway, both versions of the code are valid, sometimes one works better than the other, sometimes not Feb 11 '17 at 8:40