Academic Articles & Supplements

Ten Algorithms for Egyptian Fractions

David Eppstein

Author

David Eppstein

Title

Ten Algorithms for Egyptian Fractions

Description

The ancient Egyptians used a number system involving sums of distinct unit fractions. Mathematica is used to explore the many techniques that have been devised for translating rational numbers into these sums.

Ten Algorithms for Egyptian Fractions

by
David Eppstein

Mathematica in Education and Research
Volume 4 Number 2
Copyright 1995 TELOS/Springer-Verlag

When we use fractional numbers today, there are two ways we usually represent them: as fractions (ratios of integers) such as 5/6, and as decimal numbers such as 0.8333.
Computers typically use binary versions of either of these two representations. But these are not the only possibilities. The ancient Egyptians used a third method: instead of writing down a single fraction, they would write a sum of several distinct unit fractions, each having numerator one. For instance the Egyptians would have written 5/6 as 1/2 + 1/3 (of course, they would have used hieroglyphics instead of Arabic numerals). Today such sums are known as Egyptian fractions. (We will see another important modern representation, continued fractions, later.)

Any number has infinitely many Egyptian fraction representations, although there are only finitely many having a given number of terms [Ste92]. It is not known how the Egyptians found their representations, but today many algorithms are known for this problem, each behaving differently in terms of the number of unit fractions produced, the size of the denominators of the fractions, and the time taken to find the representations. For a good but brief introduction to Egyptian fraction algorithms and their implementation in Mathematica, see Wagon's book [Wag91]. Here we examine a number of algorithms in more detail, implement them, and analyze their performance. We also include some investigations into how many unit fractions are needed to represent rational numbers having small numerators.

We will represent Egyptian fractions in Mathematica simply as a list of unit fractions. The original rational number represented by such a list can be recovered by Apply[Plus,%]. Throughout we use q to denote the rational number we are trying to represent, or x/y when we want to talk about its numerator and denominator separately.

Methods Based on Approximation

The most natural and obvious method of finding an Egyptian fraction representation for a given number is to approximate the number as closely as possible by a single unit fraction, and then to use the same method to represent the remainder. For instance, the largest unit fraction less than 5/6 is 1/2, and removing 1/2 from 5/6 leaves 1/3, so this idea leads to the representation 1/2+1/3 mentioned above. There are several ways of translating this idea into a specific algorithm.

The Greedy Method

The greedy method produces an Egyptian fraction representation of a number q by letting the first unit fraction be the largest unit fraction less than q, and then continuing in the same manner to represent the remaining value. If q>1, we first separate out the integer part Floor[q] before representing the remaining fraction. Our implementation works by first computing a list of the fractions left after subtracting each successive term in the greedy representation, and then subtracting a shifted copy of this list from itself.

GreedyPart[q_Integer] := 0;GreedyPart[Rational[1,y_]] := 0;GreedyPart[q_Rational] := q - If[q < 0 || q > 1, Floor[q], Rational[1,1+Quotient[1,q]]]; SubtractShifted[l_] := Drop[l,-2] - Take[l,{2,-2}];EgyptGreedy[q_] := SubtractShifted[FixedPointList[GreedyPart,q]]

Let us now make sure that this routine correctly finds an Egyptian fraction representation, and analyze its performance. If we start with x/y, the remaining value after one step is
(y mod x) / y(Ceiling[y/x]), which has a smaller numerator. Similarly, the numerator decreases after each further step, so the algorithm always halts. The resulting sequence of fractions clearly adds to the original input, so the only way this method could go wrong would be if two of the fractions were equal (this is not allowed in Egyptian fractions). But this can't happen, since the denominators of the fractions must be strictly increasing: if we had two successive terms 1/a and 1/b with b † a, we could have chosen 1/(a-1) instead of 1/a.

Since the numerator decreases after each step, the number of terms in the representation of x/y is at most x. In many cases we could expect each successive numerator to be randomly distributed modulo the previous numerator; if this were really true we would instead only expect to see O(log x) terms. The denominator is at most squared each step, so the largest denominator is at most y^(2^x) or more generally y^(2^k) where k is the number of terms. The number of operations performed by the algorithm is proportional to k, but some of these operations might involve arithmetic on very large numbers. We demonstrate this method with a simple example.

EgyptGreedy[18/23]

That example was fairly well behaved; Wagon [Wag91] suggests trying this method on 31/311, which produces a representation with 10 terms, the maximum denominator having over 500 digits. (As we will see later, other methods produce much better representations for 31/311.)

We can easily modify our code to test the assertion that the numerators of the fractions remaining at each step do indeed decrease.

EgyptGreedyNumerators[q_] :=
Numerator[Drop[FixedPointList[GreedyPart,q],-2]]

EgyptGreedyNumerators[18/23]

The Harmonic Method

The greedy method treats the integer and fractional parts of a number differently. Instead, we can always remove the largest unit fraction that is smaller than both x/y and the previously removed unit fraction, even if x/y is larger than one. We treat this separately from the greedy method as it must be implemented somewhat differently — FixedPointList now needs two values, the remaining fraction and the bound on the denominator. Once we have found the sequence of remaining fractions, we remove the denominator bounds by Transpose (faster than applying First to each member of the list) and subtract the shifted list from itself as before. Our function takes two arguments, the first being the number we want to represent and the second being the largest denominator already included in the representation. The same method can also be used to generate Egyptian fractions in which the first term is arbitrarily small, simply by supplying a large integer in the second argument.

HarmonicPart[{0,d_}] := {0,d};
HarmonicPart[{Rational[1,y_],d_}] := {0,d};
HarmonicPart[{q_,d_}] :=
Max[d,1+Quotient[1,q]] //
{q-1/#,#+1} &;

EgyptHarmonic[q_,d_] :=
Transpose[FixedPointList[HarmonicPart,{q,d}]][[1]] //
SubtractShifted

The algorithm constructs a fragment of the harmonic series 1/2+1/3+1/4+1/5+... until this would produce a result larger than the original input, at which point the algorithm switches to the Greedy method for the remainder. This switch must happen after at most Exp[O(x/y+Log d)] terms, because the Harmonic series diverges (the sum up to 1/k is roughly log k). Therefore the correctness of the algorithm follows from the same analysis we saw before. However it may produce many more terms, with larger denominators, than the greedy method. Each step at most squares the denominator, so when the switch happens, the denominator of the remaining fraction can be at most doubly exponentially small with respect to x/y, and the eventual number of terms is doubly exponential in x/y (singly exponential in d). By the same analysis as the greedy algorithm, the largest denominator of the eventual representation is then at most quadruply exponential in x/y and triply exponential in d. As for the greedy method, this is only the worst case, and we can expect in practice to see one fewer level of exponentials in both the number of terms and the largest denominator. Even so, this algorithm tends to produce large representations.

EgyptHarmonic[18/23,5]

The Odd Greedy Method

Each x/y with y odd is known to have an Egyptian fraction representation in which each denominator is odd [Bre54,Ste54]. Conversely, if y is even, at least one of the terms in its representation must also be even. The most straightforward method of finding an odd-denominator representation seems to be to modify the greedy method to only use odd denominators, but it is not known whether this really works.

OddGreedyPart[{0,d_}] := {0,d};
OddGreedyPart[{Rational[1,y_],d_}] := {0,d};
OddGreedyPart[{q_,d_}] :=
Max[d,1+Quotient[1,q]] //
If[OddQ[#],#,#+1] & //
{q-1/#,#+1} &;

EgyptOddGreedy[q_,d_:3] :=
Transpose[FixedPointList[OddGreedyPart,{q,d}]][[1]] //
SubtractShifted

Unlike the greedy method, the numerators of the remaining fractions do not decrease at each step.
There are two reasons: first, like the harmonic method, we use a parameter d to make sure that the fractions we generate are distinct; d is used until the series 1/3+1/5+1/7+1/9+... becomes larger than q, at which point it becomes unimportant, but in this stage the numerators will in general increase. Second, whenever parity forces us to use a larger denominator than the greedy method, the denominator will again increase. We now give an example with q<1/3 to demonstrate the second phenomenon.

EgyptOddGreedyNumerators[q_,d_:3] :=
Transpose[FixedPointList[OddGreedyPart,{q,d}]][[1]] //
Numerator[Drop[#,-2]] &

EgyptOddGreedy[10/39]

EgyptOddGreedyNumerators[10/39]

Proving whether this method always halts remains an important open problem in the theory of Egyptian fractions [Guy81,KW91]. A heuristic argument shows that the answer is likely to be positive. After enough fractions have been generated for d to become unimportant, each step reduces the remaining fraction from some value x/y to a smaller fraction in which the numerator is between 1 and 2x-1. If each successive numerator were randomly distributed in this range, we would expect to see the numerators decrease by a factor of Exp[Integrate[Log[x],{x,0,2}]/2] ¯ 0.73576 per step. Therefore we should expect the algorithm to produce roughly (Log n)/(1-Log 2) ¯ 3.26 Log n unit fractions before halting, where n is the numerator of the remaining fraction at the point that d becomes unimportant. Of course nothing here is actually random, which is why this argument is not rigorous. (Also it ignores the possibility that the numerator and denominator of the fraction remaining after some steps may have a common factor, but that only serves to reduce the number of terms.) To test this argument, we compare it with the actual performance of our algorithm.

TestOddGreedy[q_] := EgyptOddGreedyNumerators[q] // ListPlot[Log[#]/(1-Log[2]), PlotJoined->True, AxesOrigin->{0,0}]&

TestOddGreedy[1999999991/123412340001]

Our code normalizes the vertical axis to match our heuristic prediction of the number of steps remaining. At least for this example, the numerators seem to decrease much more quickly than our prediction, so the number of terms generated (15) is considerably smaller than the 70 we would expect. It is also interesting to note the large drops made by the numerator in the third, seventh, and final steps. A closer inspection reveals that these phenomena are due to cancellation between common factors of the numerators and denominators of intermediate terms: these three steps involve common factors of 63, 45, and 5739, and two other steps involve factors of three.
It is not clear to me why such large cancellations should occur.

Conflict Resolution Methods

We next examine two methods for Egyptian fraction representation that employ the following simple idea: from a fraction x/y we can form a representation in unit fractions by making x copies of 1/y. This is not an Egyptian fraction since the unit fractions are not distinct. However we can now search for conflicting pairs (two copies of the same fraction) and resolve the conflict by replacing the pair with some other fractions adding to the same value. The methods differ in the way they choose the replacement fractions. It is trivial to prove that such methods give correct representations, but it may be harder to prove that they always halt or to analyze how well they perform.

The Pairing Method

This method uses the conflict resolution idea above. Whenever we have a conflicting pair (two copies of some fraction 1/y), we replace them either by a single fraction 2/y if y is even, or by 2/(y+1)+2/(y(y+1)) if y is odd. (Note that in all cases, the fractions simplify to have unit numerators.) The order in which this is done does not matter. Note that this process may combine pairs of fractions to form integers; e.g. this happens with sufficiently many copies of 1/7. If this happens, we allow integers to be combined to make larger integers.
This type of method is a natural fit to the pattern-matching capabilities of Mathematica, so our implementation defines a function DoPairing in such a way that Mathematica repeatedly transforms its argument list using the replacement defined above.

DoPairing[p___,q:Rational[1,y_],q_,r___] := If[OddQ[y], DoPairing[p,2/(y+1), 2/(y(y+1)),r], DoPairing[p,2/y,r]];DoPairing[p___,q_Integer,r_Integer,s___] := DoPairing[p,q+r,s]; SetAttributes[DoPairing, Orderless];EgyptPairList[l_] := Reverse[List @@ DoPairing @@ l];EgyptPairing[Rational[x_,y_]] := EgyptPairList[Table[1/y, {x}]]

Each replacement of 1/y+1/y by 2/y reduces the number of terms, initially x, by one, which can happen at most x times. Each other replacement leaves the number of terms the same but reduces the list of terms in lexicographic order; one can only perform such reductions a finite number of times. Therefore the algorithm eventually halts, with a representation having at most x terms.

Next let us determine the largest denominator that can arise. One of the fractions must be at least 1/y, and in general if the remainder after the first few terms is a/b, the next largest fraction in the representation must be at least a/xb. So if we remove the fractions from the final representation in order by size, then at each step the denominator is at most increased to its square times x, and the largest denominator is at most (xy)^(2^x). But this seems somewhat pessimistic ± with the heuristic assumption that equal fractions are not usually generated from different starting pairs, we get at most x replacements and in this case the largest denominator is roughly y^x (or even fewer if some denominators of intermediate terms are divisible by two).

EgyptPairing[18/23]

Perhaps more important than the direct use of this method for finding Egyptian fractions is the following fact, which shows that if we want to find a representation with few terms, it suffices to represent the given number as a sum of unit fractions without worrying about distinctness.

Theorem: Let q be represented as a sum of t unit fractions, not necessarily distinct.
Then q has a t-term Egyptian fraction representation.

Proof: apply the function EgyptPairList defined above to the given representation.
Each step leaves the sum of the fractions unchanged, and either shrinks the list by one
fraction or leaves its length unchanged. The fact that this halts can be shown by the same
argument given for the termination of EgyptPairing.

It would be of interest to bound the number of replacement steps performed by EgyptPairList and EgyptPairing.

The Splitting Method

The next method we describe is similar to the pairing method, but less clever: we keep a list of unit fractions as before, and resolve conflicts by replacing fractions with smaller fractions adding to the same quantity. However, instead of replacing 2/y with 2/(y+1) + 2/(y(y+1)), we replace it with 1/y + 1/(y+1) + 1/(y(y+1)). In other words, when two fractions conflict, we leave one of them in place and split the other one, creating a list with one more fraction than before.

DoSplitting[p___,q:Rational[1,y_],q_,r___] := DoSplitting[p,q,1/(y+1),1/(y(y+1)),r]; SetAttributes[DoSplitting, Orderless];EgyptSplitting[Rational[x_,y_]] := Reverse[List @@ DoSplitting @@ Table[1/y, {x}]]

It is not obvious that this method halts, but this has been proven by Graham and Jewett [Wag91]; see also Beeckmans [Bee93]. If no fraction arises in two different ways (once as 1/(y+1) and once as 1/(y(y+1)), we could analyze the algorithm on input x/y as having x-1 levels of splitting, each of which essentially doubles the number of terms in the representation. The total number of terms produced would then be O(2^x), and the largest denominator would be O(y^(2^x)). This is a best-case analysis; in practice the results will be even worse.

EgyptSplitting[5/6]

Methods Based on the Binary Number System

The Binary Method

An Egyptian fraction representation can be formed from the binary representation of a number; e.g. 27/22 = 1.00111010001. (The underscored digits repeat in blocks of 10 after these initial digits.) The initial part before the underscores gives fractions 1/2^a, and the repeating part gives fractions 1/(2^a (2^b - 1)), where b is the length of the repetition. We must take some care that the a in the second type of fraction is nonnegative, so we modify the representation above so that there are as many nonrepeating terms as repeating terms:
27/22 = 1.0011101000101110100.

A similar technique works for some other bases than binary. For instance the only digit that causes any trouble in a base 6 representation is 5, but 5/6 = 3/6+2/6 so we can still use this method with base 6. On the other hand this method does not work well with decimal notation as we can not represent 4, 7, 8, or 9 as sums of distinct divisors of 10, so numbers with those digits in their decimal representation would cause a problem for a decimal version of this method.

To implement the binary method, we first define a function to find the binary (or other base) representation of q, returned as two lists of digits. The first member of the first list is the integer part of q, the rest of the first list is the nonrepeating part of the representation, and the second list is the repeating part. It turns out to be easier to find, instead of the digits themselves, certain values mod y from which the digits can be computed. This makes it easier to detect repeating blocks of digits, since they occur exactly when the same value mod y arises twice.

RationalDigits[q_Integer, base_] := {{q},{0}};RationalDigits[Rational[a_,b_], base_Integer] := Module[{nextunit,addunit,units, reppos,breakpt,finddigit,digitize}, nextunit = (Mod[base Last[#], b]&); addunit = (Module[{c=nextunit[#]}, If[MemberQ[#,c], #, Append[#,c]]]&); units = FixedPoint[addunit, {Mod[a,b]}]; reppos = Position[units, nextunit[units]]; breakpt = reppos[[1]][[1]]-1; finddigit = (Floor[base # / b]&); digitize = (finddigit /@ # &); {Prepend[digitize[Take[units,breakpt]],Floor[a/b]], digitize[Drop[units,breakpt]]}];

Once we have found the repeating binary representation of a fraction, it is simple to turn the nonzero digits of the representation into terms in an Egyptian fraction representation. Most of the complication in our implementation is due to the point noted earlier, that we should have at least as many nonrepeating digits as repeating digits.

EgyptBinary[q_Integer] := {q};EgyptBinary[q_Rational] := Module[{l = RationalDigits[q,2], tpow = ({2 #1[[1]], #2}&), invprod = (#[[2]]/#[[1]]&), tplist,invlist, firstlen,firstlist,firstpart, mul,seclist,secpart,full}, tplist = (FoldList[tpow, {1,#[[1]]}, Drop[#,1]]&); invlist = (invprod /@ tplist[#])&; firstlen = Max[Length[l[[1]]],Length[l[[2]]]]; firstlist = Take[Apply[Join,l],firstlen]; firstpart = invlist[firstlist]; mul = 2^Length[l[[2]]]-1; seclist = RotateRight[l[[2]], Length[l[[1]]]]; secpart = (#/mul& /@ invlist[seclist]); full = (If[#==0,{ }, #]& /@ Join[firstpart,secpart]); Flatten[full]];

The correctness of this algorithm is straightforward. The fact that it halts for input q=x/y follows from the fact that the list units computed in RationalDigits is a list of distinct values modulo y, so can have length at most y. The lists of binary digits for x/y have between them at most y elements, and the padding to make the repetition start far enough along at most doubles this. At most y/2 elements on the list can correspond to binary ones, so the eventual Egyption fraction has at most y terms. All denominators are at most 2^(2y).

EgyptBinary[27/22]

The Binary Remainder Method

Let p be a power of two larger than the denominator of x/y. By dividing xp by y we find r and s satisfying x p = s y + r. Then we can represent r/p and s/p as sums of inverse powers of two; but x/y = s/p + r/(p y) so by concatenating the representation of s/p with 1/y times the representation of r/p we get a representation of x/y. The division by y ensures that no overlap occurs between the fractions from the two parts of the representation. For convenience of implementation we call EgyptBinary to find the binary representations of r/p and s/p.

BiggerPower[bound_, base_] := base^Length[IntegerDigits[bound,base]];EgyptBinRem[Rational[x_,y_]] := Module[{p, r, s}, p = BiggerPower[y,2]; r = Mod[x p, y]; s = Quotient[x p, y]; Join[If[s==0,{},EgyptBinary[s/p]], If[r==0,{},(#/y&) /@ EgyptBinary[r/p]]]]

The method produces at most Log x + Log y terms; in practice it will typically produce half that many. Each denominator is at most 2y^2.

EgyptBinRem[18/23]

The binary remainder method appears in a paper of Stewart [Ste54], where he uses it to find representations with all denominators even. Similar methods that replace the term p=2^k by some other value have proven useful in many recent results about Egyptian fractions. Breusch [Bre54] and Stewart [Ste54] set p to small multiples of 3^k, to show that every odd-denominator fraction has a representation with all denominators odd. Tenenbaum and Yokota [TY90] use factorial values of p to find representations with (1+o(1))(Log y) / (Log Log y) terms having all denominators bounded by O(y (Log y)^2 / (Log Log y)). Vose [Vos85] uses an even more complicated value of p to show that any x/y has a representation with only O(Sqrt Log[y]) terms.

Continued Fraction Methods

The Continued Fraction Method

One can derive a good Egyptian fraction algorithm from continued fractions: the algorithm is quick, generates reasonably few terms, and uses fractions with very small denominators [Ble72].

Any real number q can be represented as a continued fraction:

in which all the values a[i] are integers. This terminates in a finite sequence if and only if q is rational.

The convergents of q are formed by truncating the sequence; they are alternately above and below q, and are useful for finding good rational approximations to the original number.
(For instance the famous approximation 355/113 ¯ „ can be found as a convergent in this way.) Successive convergents have differences that are unit fractions. The sequence of these differences gives something like an Egyptian fraction representation of q, but unfortunately every other fraction in the sequence is negative.

If h[i]/k[i] denotes the ith convergent, we can define a sequence of secondary convergents:

As j ranges from 0 to a[n+1] the secondary convergents give an increasing sequence ranging from the (i-1)st convergent to the (i+1)st convergent [NZ80]. As with the primary convergents, successive secondary convergents differ by a unit fraction. If we interleave the sequence of every other primary convergent, connected by the appropriate sequences of secondary convergents, the differences of this interleaved sequence give an Egyptian fraction representation of q.

We first find the continued fraction representation of q=x/y. (Mathematica provides a package for continued fractions, but one must supply a bound on the number of terms to compute. We don't need or want such a bound.) In order to use this method, the continued fraction must have an odd number of terms, so if necessary we replace the last term a[i] with two terms a[i]-1 and 1.

CFNextTerm[q_Integer] := 0;
CFNextTerm[q_Rational] := 1/(q-Floor[q])

ContinuedFractionList[q_] :=
Floor /@ Drop[FixedPointList[CFNextTerm, q],-2];

CFMakeOdd[l_] :=
If[OddQ[Length[l]],l,Join[Drop[l,-1],{Last[l]-1,1}]]

We next find the primary and secondary sequences of unit fractions from these continued fraction representations.

CFPSAux[{a_,b_},c_] := {b + a c, a};CFPrimarySeq[l_] := Transpose[Drop[FoldList[CFPSAux,{0,1},l],1]][[1]]; CFSecondarySeq[l_] := If[Length[l] < 3, l, Table[l[[1]] + i l[[2]], {i, 0, (l[[3]]-l[[1]])/l[[2]]-1}] ~Join~ CFSecondarySeq[Drop[l,2]]]

As described above, our final representation is formed by hooking together secondary sequences. We first separate out the integer part of the input, which we leave as is. The remaining fractions are formed by multiplying pairs of values in the secondary sequence.

EgyptContinuedFraction[q_] := CFSecondarySeq[CFPrimarySeq[CFMakeOdd[ ContinuedFractionList[q]]]] // 1/(Drop[#,1] Drop[#,-1])& // If[Floor[q]==0, #, Prepend[#, Floor[q]]]&

Termination of the algorithm follows from the termination of the continued fraction representation algorithm, which is essentially the same as Euclid's algorithm for integer GCD's. It is clear from the construction of the secondary sequence, and from the fact that the final result has denominators that are products of pairs of numbers in the secondary sequence, that all fractions are distinct. The fact that the sum of the fractions is the original input number
is a straightforward but tedious exercise in algebraic manipulation. The number of terms in the Egyptian fraction representation of x/y is the sum of the odd terms after the first in the continued fraction list, which is at most x. Each fraction is a difference between two secondary convergents with denominator at most y, so each fraction has denominator at most y^2.

EgyptContinuedFraction[18/23]

The Grouped Continued Fraction Method

The worst case for the continued fraction method above occurs when the continued fraction representation has only three terms producing a long secondary sequence. In this case the Egyptian fraction representation will involve long sequences of fractions of the form
1/(a+b i)(a+b(i+1)). If we add k consecutive values in such a sequence, we get
k/(a+b i)(a + b(i + k)); it may happen that this can be simplified to a unit fraction again. By performing several simplifications, we both reduce the number of terms in the overall representation and also reduce some denominators. For instance, the continued fraction method for 7/15 gives

But 1/15 + 1/35 + 1/63 = 1/9, and 1/99 + 1/143 + 1/195 = 1/45, so we can replace these triples and find the shorter representation

This phenomenon is not unusual, and Bleicher [Ble72] showed how to take advantage of it to dramatically reduce the number of terms produced by the continued fraction method. Some care is required: if in the above list we instead group the last five terms, we get

which is not an Egyptian fraction representation.

Our implementation finds all shortest representations rather than a single representation, so if they had distinct fractions we would return both representations above. We partition the secondary sequence into blocks of arithmetic progressions and find groupings separately within each progression; this is safe as the sum of all fractions from one progression is smaller than half of any fraction in a previous progression. Within a progression, we determine which groups of terms can be combined to form a unit fraction, and represent each group as an edge in a graph, labelled with the corresponding unit fraction. For the example above, the graph has eight vertices and ten edges, as follows:

Each edge is directed from left to right. The horizontal edges represent the original terms produced by the continued fraction method, while the longer edges represent the groupings that result in unit fractions. Our task then becomes one of finding the shortest path through this graph, with the restriction that we cannot use two edges with the same label.

Unfortunately finding paths without repeated labels is NP-complete, so an efficient algorithm for this subproblem is unlikely to exist. Fortunately most of the time our graphs have few repeated labels and the problem is not as hard as its worst case. We use the following heuristic: for increasing values of k, find all paths of k or fewer edges, and filter out the paths with repeated labels; if not all paths are filtered out, return the remaining list of paths. The theoretically fastest algorithm for listing all short paths takes constant time per path, after preprocessing time proportional to the time to find a single shortest path [Epp94], however for ease of implementation we use a simpler method invented by Byers and Waterman [BW84]. (The motivation of both papers was not Egyptian fractions, but rather comparison of DNA and protein sequences; this also turns out to be equivalent to a certain shortest path problem.)

First we include code to make an adjacency matrix for a graph, containing in each entry either the fraction corresponding to an edge in the graph, or the empty set if no such edge exists (i.e. if the corresponding sum of terms does not reduce to a unit fraction). The input to this routine is the secondary sequence of the continued fraction.

ECFMakeGraph[l_] := Table[If[i<j, (If[Numerator[#]===1,#,{}]&) [Rational[j-i, l[[i]] l[[j]] ]], {}], {i,1,Length[l]}, {j,1,Length[l]}]

Next we include a shortest path algorithm, which takes as input the adjacency matrix above and produces a vector of distances from vertices to the last vertex. This vector is needed for our bounded length path search.

ECFPathLengths[g_] := ECFPathLengths[g,Length[g]-1,{0}];ECFPathLengths[g_,i_,vec_] := Prepend[vec, Min@@Table[If[g[[i,j]]==={}, Infinity, vec[[j-i]]+1], {j,i+1,i+Length[vec]}]] // (If[i===1,#,ECFPathLengths[g,i-1,#]]&);

We now implement Byers and Waterman's algorithm for finding all paths that contain at most b more edges than are in the shortest path itself. We will call this algorithm repeatedly, using larger and larger values of b, until we find a path without repeated labels. Our implementation takes as input the graph, the value of b, the vertex to start at, the number of vertices, and the vector of distances produced above, but all but the first two can be omitted (in which case we supply appropriate values automatically).

The technique is simply to build the path one edge at a time. At each step we compute a value d measuring the amount by which the path length would increase if we followed the given edge instead of keeping to the shortest path (d=0 for shortest path edges). We subtract d from b and continue recursively as long as the result is nonnegative.

ECFBoundedPaths[g_,b_] := ECFBoundedPaths[g,b,1,Length[g],ECFPathLengths[g]];ECFBoundedPaths[g_,b_,i_,l_,v_] := If[i===l,{{}}, Join@@Table[If[g[[i,j]]==={},Infinity, 1+v[[j]]-v[[i]]] // (If[#>b,{},(Prepend[#,g[[i,j]]]&) /@ ECFBoundedPaths[g,b-#,j,l,v]]&), {j,i+1,l}]]

We next include code for removing from the list those paths that contain a duplicated fraction.
It is not clear that the paths will have the fractions listed in sorted order, so we sort them first.

ECFContainsDupl[{___,q_,q_,___}] := True;ECFContainsDupl[l_] := False; ECFFilterDuplSub[x_] := If[ECFContainsDupl[x],{},{x}];

ECFFilterDupls[l_] := Join @@ (ECFFilterDuplSub[Reverse[Sort[#]]]&) /@ l; ECFShortFilter[g_] := ECFShortFilter[g,ECFPathLengths[g],0];ECFShortFilter[g_,v_,b_] := ECFFilterDupls[ECFBoundedPaths[g,b,1,Length[g],v]] // (If[#==={},ECFShortFilter[g,v,b+1],#]&);

The next function applies all of the above steps for three-term continued fractions. The final algorithm applies this to several three-term subsequences of the whole continued fraction.

ECFArithSeq[a_,b_,c_]:=ECFShortFilter[ ECFMakeGraph[CFSecondarySeq[{a,b,c}]]]

The next function takes two lists of lists, and forms all pairwise concatenations of one item from the first list and one from the second. The obvious approach of using Outer[Join,...] doesn't work since Outer interprets lists of lists as tensors, so we use an alternate method based on Distribute.

OuterJoin[ll_,mm_] := Distribute[{ll,mm},List,List,List,Join];

We are finally ready to define the overall modified continued fraction method, which breaks the primary sequence into subsequences and calls ECFArithSeq on each one.

ECFSecondaryPaths[l_] := If[Length[l]<3,{{}}, OuterJoin[ECFArithSeq[l[[1]],l[[2]],l[[3]]], ECFSecondaryPaths[Drop[l,2]]]]

EgyptGroupedCF[q_] := ECFSecondaryPaths[CFPrimarySeq[ CFMakeOdd[ContinuedFractionList[q]]]]

Every step involves a fixed number of nested loops with indices bounded by the length of the secondary sequence, so (with the possible exception of finding a short repetition-free path) the overall time is polynomial in the numerator of the original rational number given as input.
It is not hard to see that the algorithm produces sequences of fractions formed by grouping the results of the continued fraction method, so the sum of the sequence is correct. It remains to verify that no fraction is duplicated. This is checked explicitly within each subsequence, and the entire sum of any subsequence is less than half any single fraction in previous subsequences, so no two separate subsequences can produce duplications.

As in the continued fraction method, the largest denominator in the representation of x/y is O[y^2]. The number of terms is still O[x] but it can also be analyzed in terms of y.
Bleicher [Ble72] shows that by choosing a prime p with gcd(a,p)=1 and p=O(log a),
and using groups with sizes equal to powers of p, one can find a representation with
O(p Log[b]/Log[p]) = O(Log[x]Log[y]/Log Log[y]) terms.
Since the actual representation is chosen to have minimum length, it can be no longer than this.
It remains unclear whether the implementation above really takes polynomial time, or whether there can be sufficiently many repeated labels that the algorithm for listing short paths has to list a large number of paths and slows down to exponential. However in practice this method seems to work well. (Bleicher's method of grouping can apparently be done in polynomial time.)

EgyptGroupedCF[31/311]

The graph constructed for 31/311 is too complicated to depict here. It has two paths of length five; however one of the paths is eliminated because it has two copies of the label 1/231.

A Hybrid Pairing / Continued Fraction Method

We can use potentially even fewer terms than the grouped continued fraction method, at the expense of possibly increasing the maximum denominator in the representation. We simply find shortest paths in the same graph constructed by that method, ignoring the possibility of repeated labels, and then make the unit fractions in the resulting representation distinct by applying EgyptPairList.

EHArithSeq[a_,b_,c_] := ECFBoundedPaths[ ECFMakeGraph[CFSecondarySeq[{a,b,c}]],0]EHSecondaryPaths[l_] := If[Length[l]<3,{{}}, OuterJoin[EHArithSeq[l[[1]],l[[2]],l[[3]]], EHSecondaryPaths[Drop[l,2]]]]EgyptHybrid[q_] := EgyptPairList /@ EHSecondaryPaths[CFPrimarySeq[ CFMakeOdd[ContinuedFractionList[q]]]];

This method uses O(Log[x]Log[y]/Log Log[y]) terms to represent any number x/y. In the following example, we see representations corresponding to both shortest paths in the graph constructed for 31/311.

EgyptHybrid[31/311]

Small Numerators

The algorithms described above work for any input. We now discuss techniques limited to specific numerators. The typical question here is how many terms are needed to represent fractions with a given numerator. For fractions 2/y the answer is clearly 2. Some fractions 3/y require 3 terms, as we see below. It is not known whether any fraction 4/y requires 4 terms.

More generally, good bounds are known on the number of terms needed to represent x/y measured as a function of y [Vos85], but there seems to be less work on measuring this minimum number of terms as a function only of x. As we note in the section on 4/y, a solution to this specific case would have implications for the general problem.

Numerator 3

The basic result for fractions of the form 3/y is that there is a two-term expansion if and only if y has a factor congruent to 2 mod 3. Klee and Wagon [KW91] credit this result to N.˚Nakayama; however they supply no citation, so we repeat the proof below.

Theorem: 3/y has a two-term expansion if and only if y has a factor congruent to 2 mod 3.

Proof: In one direction, the representation 3/(3n+2)=1/(n+1)+1/(n+1)(3n+2) is found by both the greedy and continued fraction methods. This idea can easily be extended to 3/y where y is a multiple of 3n+2.

In the other direction, suppose y=3n+1 and 3/y=1/a+1/b=(a+b)/ab. First note that a and b must be divisible by the same power of 3, since if a were divisible by 3^i and b by 3^j, with j>i, then a+b would not divisible by 3^j and the powers of 3 wouldn't cancel from the denominator. Let g=gcd(a,b), u=a/g, v=b/g, so 3/y=(u+v)/guv and 3 divides u+v; let u+v=3z. Then 1/a+1/b=3z/guv and g must factor as zw since gcd(uv,u+v)=1. So y=uvw. For u+v=3z, one of u and v (say u) must be 2 mod 3, giving the factor of y we seek.

Unfortunately this seems to imply that finding short representations, even in this special case, is computationally difficult: at least as difficult as factoring integers.

Some examples of two-term representations that would not be found by our general algorithms: 3/25=1/10+1/50; 3/55=1/22+1/110=1/20+1/220; 3/121=1/44+1/484

Numerator 4

The question of whether all fractions 4/y have 3-term representations is discussed by Mordell [Mor69], who attributes it to Erdos and Straus. Guy [Guy81] cites several other authors as having worked on the problem: Bernstein, OblÇth, Rosati, Shapiro, Straus, Yamamoto, and Franceschine. Others have worked on more general versions of this problem including Schinzel, Sierpinski, SedlÇcek, PalamÈ, Stewart, Webb, Breusch, Graham, and Vaughan.
A positive solution to this question would have more general implications: we could use such a solution as the basis for a conflict resolution method that, given a number x/y, would find an Egyptian fraction representation with x^(Log[3]/Log[4]) ¯ x^0.7925 terms.

Modular Conditions

Mordell shows that in any example 4/y requiring 4 terms in an Egyptian fraction representation, y must be 1 mod 24, –1 mod 5, and one of three values mod 7 (giving a total of 6 possible values mod 840, all squares of small numbers). If y is a minimal counterexample, it must be prime (since if y=ab we could divide all terms in a representation of 4/a by b).

If y is 2 or 3 mod 4, the greedy algorithm gives a 2 or 3 term representation. If y is 3 mod 6, we have the representation 3/y+1/y. And if y is 5 mod 6, we have the representation 1/ceil[y/6] + 3/(y ceil[y/6]) where the last term has a 2-term expansion by our previous analysis. So if 4/y is to fail to have a 3 term representation, y must be of the form 24n + 1. Several methods extend this analysis by representing 4/y when y (equivalently n) has certain values modulo small primes.

The representations 1/(6n+1) + 3/(24n+1)(6n+1) and 1/(18n+1)(24n+1) + 3/(18n+1) work if one of 6n+1, 18n+1, or 24n+1 is divisible by a prime p congruent to 5 mod 6. Thus for any of these primes one can derive rules for finding three-term representations of 4/y, that work whenever y has certain values mod p. We can use this technique to find representations when n is congruent to 4, 3, or 1 mod 5 (and so, rule out counterexamples for y congruent to anything but –1 mod 5).

The representation 1/(6n+k) + (4k-1)/(6n+k)(24n+1) works via a greedy method if a factor of the second denominator is (4k-2) mod (4k-1), or more generally if the factor is (4k-1-i) mod (4k-1) and i divides the denominator. In particular these work with k=2 when n is 2, 3, 4, or 6 mod 7 (with the corresponding values of i being 0, 1, 1, and 2). Therefore in any counterexample 4/y, y must be a quadratic residue mod 7.

Yet another type of rule is possible: consider the decomposition
1/(6n+k) +a/(6n+k)(24n+1) + b/(6n+k)(24n+1), where a+b = 4k-1. This is only possible when k is even, since otherwise one of a or b would be even and not divide the denominator.
For instance 4/(24n+1)=1/(6n+10) + 26/(6n+10)(24n+1) + 13/(6n+10)(24n+1)
where the last two simplify to unit fractions if n is 7 mod 13.

Particular Values

As noted above, the numbers y for which 4/y might possibly require four terms fall into six classes modulo 840: 1, 121, 169, 289, 361, and 529. We only need to consider prime n since if mn is a counterexample, so must be both m and n. Following are representations for all such cases through 12500. Most use rules like the ones described above that depend only on the values of y mod 11, 13, and 19, but 4/3361 uses a method that depends on y mod 29 and 4/8089 uses a method that depends on y mod 17.

4/1801 = 1/451 + 1/295364 + 1/3249004

4/2521 = 1/636 + 1/69748 + 1/131876031

4/2689 = 1/676 + 1/139828 + 1/908882

4/3049 = 1/772 + 1/60980 + 1/5884570

4/3361 = 1/841 + 1/974690 + 1/28266010

4/3529 = 1/892 + 1/80726 + 1/569764108

4/3889 = 1/975 + 1/345150 + 1/268457670

4/4201 = 1/1096 + 1/25208 + 1/13237351

4/4561 = 1/1244 + 1/13684 + 1/15603181

4/4729 = 1/1185 + 1/510732 + 1/201739140

4/5209 = 1/1308 + 1/296262 + 1/3086457516

4/5569 = 1/1402 + 1/200484 + 1/140539284

4/5881 = 1/1604 + 1/17644 + 1/25941091

4/6841 = 1/1713 + 1/1065486 + 1/7288989726

4/7681 = 1/1924 + 1/1136788 + 1/7389122

4/8089 = 1/2023 + 1/5775546 + 1/98184282

4/8521 = 1/2324 + 1/25564 + 1/54457711

4/8689 = 1/2175 + 1/1718250 + 1/14929874250

4/8761 = 1/2196 + 1/836676 + 1/3665059218

4/8929 = 1/2233 + 1/7250348 + 1/79753828

4/9241 = 1/2314 + 1/1644898 + 1/10691837

4/9601 = 1/2406 + 1/1008105 + 1/269500070

4/9769 = 1/2452 + 1/614226 + 1/12000747588

4/10369 = 1/2828 + 1/31108 + 1/80639713

4/12049 = 1/3016 + 1/2795368 + 1/18169892

4/12289 = 1/3078 + 1/1644678 + 1/30317171913

According to Guy, N. Franceschine has performed similar calculations for y<10^8.

References

[Bee93]
L. Beeckmans. The splitting algorithm for Egyptian fractions. J. Number Th. 43, 1993, pp. 173--185. This paper contains a proof that the splitting method terminates; Wagon [Wag91] credits the same result to Graham and Jewett.

[Ble72]
M. N. Bleicher. A new algorithm for the expansion of continued fractions. J. Number Th. 4, 1972, pp. 342--382. Defines two methods for Egyptian fraction representation: what he calls the Farey sequence method, which is equivalent to the continued fraction method described here, and what he calls the continued fraction method, which is a variant of what we call the grouped continued fraction method.

[Bre54]
R. Breusch. A special case of Egyptian fractions, solution to advanced problem 4512. Amer. Math. Monthly 61, 1954, pp. 200--201. Shows that every x/y with y odd has an Egyptian fraction representation with all denominators odd, by using a method similar to the binary remainder method but using the base 5(3^k) instead of 2^k.

[BW84]
T. H. Byers and M. S. Waterman. Determining all optimal and near-optimal solutions when solving shortest path problems by dynamic programming. Oper. Res. 32,1984, pp. 1381--1384. Byers and Waterman describe a simple algorithm for finding short paths in a graph, used in our implementation of the grouped continued fraction method.

[Epp94]
D. Eppstein. Finding the k shortest paths. 35th IEEE Symp. Foundations of Computer Science, 1994, pp. 154--165. I describe what is theoretically the fastest known algorithm for the shortest path problem used in the grouped continued fraction method, however my technique is rather more complex than that of [BW84] and has not been implemented.

[Guy81]
R.K. Guy. Unsolved Problems in Number Theory. Springer-Verlag, 1981, pp. 87—93. Guy lists a number of questions about Egyptian fractions, including the following:
Does 4/y have a 3-term representation for all y?

Does 5/y? Does x/y for x sufficiently large relative to y?

Does the odd greedy method terminate?

What different denominators are possible in a t-term representation of one?

Do all positive-density sets of integers have a subset forming the denominators of a representation of one?

If we assign to each of the integers one of a finite set of colors, can we pick a single color that can be used as the denominators for representations of all rationals? Guy also

[KW91]
V. Klee and S. Wagon. Old and New Unsolved Problems in Plane Geometry and Number Theory. Math. Assoc. of America, 1991, pp. 175--177 and 206—208. Concentrates primarily on the question of whether the odd greedy method halts; notes that a similar method for finding representations with even denominators does halt; contains a number of further references.

[Mor69]
L. J. Mordell. Diophantine Equations. Academic Press, 1969, pp. 287--290. Discusses the question of whether 4/y has a three-term representation; describes methods of finding such representations depending on the value of n modulo various small primes.

[NZ80]
I. Niven and H.S. Zuckerman. An Introduction to the Theory of Numbers, 4th ed., Wiley, 1980, p. 200. They describe in an exercise here the secondary sequences used in the continued fraction methods of Egyptian fraction representation.

[Ste54]
B. M. Stewart. Sums of distinct divisors. Amer. J. Math. 76, 1954, pp. 779--785. Uses the binary remainder method to find representations with all denominators even. Also, like [Bre54], uses a method similar to the binary remainder method (using base 35(3^k)) to find odd representations of fractions with odd denominators.

[Ste92]
I. Stewart. The riddle of the vanishing camel. Scientific American, June 1992, pp. 122--124. Stewart shows that any rational has only a finite number of t-term Egyptian fraction representations, for any fixed number t, and describes a brute force method for finding all such representations.

[TY90]
G. Tenenbaum and H. Yokota. Length and denominators of Egyptian fractions. J. Number Th. 35, 1990, pp. 150--156. This paper provides good simultanous bounds on the number of terms and the denominators of an Egyptian fraction representation, as discussed in the section on the binary remainder method.

[Vos85]
M. Vose. Egyptian fractions. Bull. Lond. Math. Soc. 17, 1985, p. 21. Shows that every x/y has a t-term representation where t=O(Sqrt Log[y]). The technique is similar to the binary remainder method.

[Wag91]
S. Wagon. Mathematica in Action. W.H. Freeman, 1991, pp. 271--277. Wagon implements the greedy and odd greedy methods, and describes the splitting method. He also mentions the open problem of whether the odd greedy method always terminates for the special case of fractions with numerator 2.

About the Author

David Eppstein received his Ph.D. in computer science from Columbia University in 1989.
After a postdoctorate at the Xerox Palo Alto Research Center, he took a position at the University of California, Irvine, where he is now an associate professor. Prof. Eppstein's primary research areas are graph algorithms and computational geometry; his work is supported by an NSF National Young Investigator award and by matching funds from Xerox Corp.

David Eppstein
Department of Information and Computer Science
University of California, Irvine CA 92717
eppstein@ics.uci.edu
http://www.ics.uci.edu/~eppstein/

Cite this as: David Eppstein, "Ten Algorithms for Egyptian Fractions" from the Notebook Archive (2002), https://notebookarchive.org/2018-10-10pt4o2