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

Question

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:

• Greedy algorithm.

• Engel expansion.

• '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.

Answer 1

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
  ]

I took your posted code and made the following changes:

• 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.